0

I am trying to simulate noble gas interactions based on the Lennar Jones Potential. And I ran into the problem that the interactions of the particles is apparently calculated in discrete steps so that the acceleration tends to either jump to unreasonably high values or doesn't change significantly when two particles approach each other. The detection algorithm that I'm using looks as follows:

def accerlation(self):
   index = 0
   for particle, particle2 in itertools.combinations(self.particles, 2):
       index += 1
       
       x_diff = (particle.r[0] - particle2.r[0])*S
       y_diff = (particle.r[1] - particle2.r[1])*S
       distance = np.sqrt(x_diff**2 + y_diff**2)
       
       if distance < criticaldistance1 and distance> criticaldistance2:
           print("Interaction") 
           particle.a[0] -= 1/(Particle().m) *(LennardJones_Force(x_diff))
           particle.a[1] -= 1/Particle().m *(LennardJones_Force(y_diff))
           print(particle.a)

(The print commands are for debugging) The Lennard Jones Force is simply a function that returns (24*epsilon*sigma**6*(Distance**6 - 2*sigma**6))/Distance**13, where epsilon and sigma are constant values.

If you can't spot a mistake there, the problem might also be in this part:

def Changeovertime(self): 
       self.Collision_Box() #elaxtic collsions with box
       for particle in self.particles:
           particle.r += self.dt * particle.v + self.dt**2*particle.a 
           particle.a_prior = particle.a 
           self.accerlation() #
           particle.a = np.zeros(2) #
           particle.v = (particle.v + self.dt/2 * (particle.a + particle.a_prior)) *self.scaling()

When I print the acceleration of a pair of particles during the animation it looks as follows:

[-3.21599405e-01 -1.05489024e-18]\
Interaction\
[-3.35299415e-14  2.61407475e-19]\
Interaction\
[-2.52825200e+31 -1.05489024e-07]\
Interaction\
[-6.70598831e-14  5.22814950e-19]\
Interaction\
[ 1.57188229e-01 -5.51566856e-19]\
Interaction\
[-6.70598831e-03  5.22814950e-08]\
Interaction\
[ 3.14376458e-01 -1.10313371e-18]\
Interaction\
[-3.35299416e-14  2.72195193e-19]\
Interaction\
[-6.70598831e-14  5.44390386e-19]\
(Note the jump from entry 2 to entry 3)
The Fabio
  • 5,369
  • 1
  • 25
  • 55
Titus
  • 11
  • 2

2 Answers2

1

Without commenting on the rest of the code, the problem you see is due to the lines:

particle.a[0] -= 1/(Particle().m) *(LennardJones_Force(x_diff))
particle.a[1] -= 1/Particle().m *(LennardJones_Force(y_diff))

The force, and the acceleration, are vectors: you can not break them into x and y components like that. You have to calculate the force and project it into x and y (you can work with acceleration if you want).

force = LennardJones_Force(distance)
particle.a[0] -= 1/(particle.m) * force * x_diff/distance
particle.a[1] -= 1/(particle..m) * force * y_diff/distance

Please check the signs (they should be opposite for particle and particle2).

The huge jump you see is because two particles had a very similar x or y coordinate (even if they were far apart), and the miscalculated force was huge.

azelcer
  • 1,383
  • 1
  • 3
  • 7
0

It seems that physics is not obeyed.
Verify that Newton's third law is implemented.\

Note that the moment these particles interact, an action-reaction pair is formed, so to ensure conservation of energy the force acting on particle j must be the same force acting on particle i (in direction and absolute value): Fx[i] = -Fx[j] (imagine i red and j blue!).

Here's a note: if the masses of these particles are equal and their values ​​are 1, you can assume that acceleration equals force Fx[i] = ax[i]/1. Then, all that remains is to decompose these forces (Fr[i] and Fr[j]) into their components (Fx[i], Fy[i] and Fz[i] for Fr[i] and Fx[j], Fy[j] and Fz[j] for Fr[j]), always remembering that the force that particle j does on particle i is equal to the force that particle i does on particle j, being different only by their signs (different senses).

In computer practice, you do it this way:

def acceleration():
    for i in range(0,N_particles):
        for j in range(i+1,N_particles):
            rij = np.sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2 + (z[i] - z[j])**2)
            if (rij < rmaxlj):
                flj = (48*epsilon*(sigma**12/rij**13) + 24*self.epsilon*(sigma**6/rij**7))/rij**2
                Fx[i] += flj*(x[i] - x[j])
                Fy[i] += flj*(y[i] - y[j])
                Fz[i] += flj*(z[i] - z[j])
                Fx[j] -= flj*(x[i] - x[j])
                Fy[j] -= flj*(y[i] - y[j])
                Fz[j] -= flj*(z[i] - z[j])
    return Fx, Fy, Fz

Note that all forces that are added to each component of particle j are "removed" from each component of particle i. This is actually Newton's third law in action!
PS.: these rmaxlj works as if you were using criticaldistance and also you can use i as 0 and j as 1

Test this and tell me later!

JailsonKbT
  • 3
  • 3