Модель потенциала Ленарда Джонса, метод молекулярной динамики

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
 
namespace Molecules // Имя проекта (решения) - Molecules
{
 
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }
 
        List<Atom> atoms_list = new List<Atom>(); // список (коллекция) экземпляров (атомов)
        int razmer, kolvo; // размер и кол-во атомов
        Random rand = new Random(DateTime.Now.Millisecond); //экземпляр класса Random        
 
        public class Atom  //класс описывающий отдельный атом
        {  
            public double velosity { get { return Math.Sqrt(dX * dX + dY * dY); } }
            public double X, Y; // Координаты центра атомов
            public double dX, dY; // Путь за 1 tick таймера
            private double rad, mass;
            public static double koeffColiss = 1; //Коэффициент коллизии (столкновения), 1 или -1, направление вперед или в противоположную сторону
            static Random rnd = new Random(DateTime.Now.Millisecond);
            static PictureBox Field; //поле PictireBox 
 
            public double radius
            {
                get { return rad; }
                set 
                {
                    rad = value;
                    mass = rad * rad;
                }
            }
 
            public double Mass
            {
                get { return mass; }
                set { rad = Math.Sqrt(mass = value); }
            }
 
            public double Velosity { get { return Math.Sqrt(dX * dX + dY * dY); } }
 
            public Atom(int Size, PictureBox field) //метод описывающий атом
            {
                Field = field; //задание поля, в котором рисовать
                radius = Size * 0.5; //расчет радиуса, в зависимости от размера, заданного пользователем
                dX = (rnd.NextDouble() - 0.5) * 5; // Путь за один tick таймера 
                dY = (rnd.NextDouble() - 0.5) * 5;
            }
 
 
            // Функция рисования
            public void drawBall(PaintEventArgs e)
            {
                e.Graphics.DrawEllipse(new Pen(Color.Blue), (float)(X - radius), (float)(Y - radius), 2 * (float)radius, 2 * (float)radius); //рисование круга
            }
 
            //рассчет траектории (отскакивание от стенок)
            public void beginAtom()
            { 
                if (X <= radius && dX < 0) dX *= -koeffColiss; // Если координата Х <= 10 то меняем направление движения по оси Х на противоположное
                else if (X >= Field.Width - radius && dX > 0) dX *= -koeffColiss;
                if (Y <= radius && dY < 0) dY *= -koeffColiss; // Если координата У <= 10 то меняем направление движения по оси У на противоположное
                else if (Y >= Field.Height - radius && dY > 0) dY *= -koeffColiss;
 
                X += dX;
                Y += dY;
            }
 
            public bool BallOwerlap(Atom ball1)
            {
                return (X - ball1.X) * (X - ball1.X) + (Y - ball1.Y) * (Y - ball1.Y) <=
                     (radius + ball1.radius) * (radius + ball1.radius); // условие столкновения
            }
 
            public void collisAtoms(Atom ball1) // Здесь рассчитывается взаимодействие шариков друг с другом
            {
                if (BallOwerlap(ball1))
                {
                    if ((dX - ball1.dX) * (X - ball1.X) + (dY - ball1.dY) * (Y - ball1.Y) < 0) //условие столкновения, учитывая dX
                    {
                        double dx = ball1.X - X, dy = ball1.Y - Y;
                        double alfa = Math.Atan2(dy, dx);
                        double angle = Math.Atan2(dY, dX) - alfa, angleBall = Math.Atan2(ball1.dY, ball1.dX) - alfa,
                        mod = Velosity, modBall = ball1.Velosity;
                        double DX = mod * Math.Cos(angle),
                        newDY = mod * Math.Sin(angle),
                        DX1 = modBall * Math.Cos(angleBall),
                        newDY1 = modBall * Math.Sin(angleBall);
                        double newDX = (DX * (Mass - ball1.Mass) + 2 * ball1.Mass * DX1) / (Mass + ball1.Mass),
                        newDX1 = (DX1 * (ball1.Mass - Mass) + 2 * ball1.Mass * DX) / (Mass + ball1.Mass);
                        angle = Math.Atan2(newDY, newDX) + alfa;
                        angleBall = Math.Atan2(newDY1, newDX1) + alfa;
                        mod = koeffColiss * Math.Sqrt(newDX * newDX + newDY * newDY);
                        modBall = koeffColiss * Math.Sqrt(newDX1 * newDX1 + newDY1 * newDY1);
                        dX = mod * Math.Cos(angle);
                        dY = mod * Math.Sin(angle);
                        ball1.dX = modBall * Math.Cos(angleBall);
                        ball1.dY = modBall * Math.Sin(angleBall);
                    }
 
                }
 
            }
        }
 
        private void button1_Click(object sender, EventArgs e) //функция-событие по нажатию кнопки
        {
            timer1.Enabled = true; //запуск таймера
 
            razmer = Convert.ToInt32(textBox2.Text); //преобразование строки из TextBox2 ("Размер атомов") в целое число
            kolvo = Convert.ToInt32(textBox1.Text); //преобразование строки из TextBox1 ("Кол-во атомов") в целое число
 
            atoms_list.Clear(); //очищение списка атомов во избежании добавления новых атомов при повторном нажатии кнопки
            int r = razmer / 2; 
            for (int i = 0; i < kolvo; i++) //цикл, пока не пройдет все атомы (kolvo)  
            {
                Atom a = new Atom(razmer, pictureBox1); //создание экземпляра конкретного атома класа Atom
                a.X = rand.Next(10, pictureBox1.Width - r); //Присваивание радномных координат созданным атомам
                a.Y = rand.Next(r, pictureBox1.Height - r); //Чтобы атом не создавался за пределами PictureBox, отнять r
                foreach (Atom a1 in atoms_list); 
                atoms_list.Add(a); //добавление созданных атомов в список
            }
 
        }
 
         public void timer1_Tick(object sender, EventArgs e) //функция, когда запускается таймер
        {
            foreach (Atom a in atoms_list) a.beginAtom();//для каждого атома из спиcка применяется функция beginAtom
            for (int i = 0; i < atoms_list.Count - 1; i++) 
                    for (int j = i + 1; j < atoms_list.Count; j++) if (atoms_list[i].BallOwerlap(atoms_list[j]))  
                    {
                    atoms_list[i].collisAtoms(atoms_list[j]); 
                    }
            pictureBox1.Invalidate(); //перерисовка PictureBox
        }
 
        private void paint(object sender, PaintEventArgs e)
        {
            foreach (Atom a in atoms_list) a.drawBall(e); // Само рисование //для каждого атома из спиcка применяется функция drawBall(e)
        }
 
    }
}

    public class LennardJonesAtom : Atom
    {
        float epsylon;
 
        public LennardJonesAtom(float radius, float epsylon) : base(radius)
        {
            this.epsylon = epsylon;
        }
 
        public override void CalcCollision(Atom other)
        {
            //calc distance between atoms
            var dist = Distance(other);
 
            if (dist > Radius * 5) //no interaction
                return;
 
            //calc force
            var k = Radius / dist;
            var f = 4 * epsylon * (float)(Math.Pow(k, 12) - Math.Pow(k, 6));//module of force (LennardJones)
 
            //to avoid infinite forces -> imitate elastic impact
            if (f > 30 || float.IsNaN(f))
            {
                base.CalcCollision(other);
                return;
            }
 
            //force vector
            var force = new PointF(X - other.X, Y - other.Y).Normalized().Mul(f);
 
            //apply force
            ApplyForce(force);
            other.ApplyForce(force.Mul(-1));
        }
    }

            //calc force
            var k = Radius / dist;

            //calc force
            var k = (Radius + other.Radius) / dist;