Simulation of a gas in 2-D using a Verlet algorithm

[gibatom]

Homework Statement

Hello everybody,

I’m currently working on an university project, aiming to understand better the kinetic of a gas of atoms (in 2D).

To do so, i’m using a Verlet algorithm, which basically consists of improving the classical method of integrating the equations of motions of all the atoms making up the gas (with Newton’s 2nd law). With this algorithm, the energy should be conservered (that’s why it was developped in the first place) - which is effectively what i see with my Python code when the gas is not confined.

Neverthelsesss i would like to modelize a gas confined in a box, in order to study values such as the pressure, the temperature at equilibrum or the repartition of atom speeds. To do so, i use periodic conditions (instead of rebounds with the box, which i thought would be more complex) : i.e when an atom leaves the box, it appears on the other side, like PacMan.
Now, the problem i have is that as soon as an atom leaves the box, comes in the other side and interacts frankly with other atoms, the energy in the box grows all of a sudden.

Has any of you already encountered this problem ? Or has someone an idea to solve this problem of the non-conservation of the energy when i use periodic conditions ?

Thank you very much in advance for your help,
(If necessary i can show you my code)

Relevant Equations

-Equations of motions (2nd law of thermodynamics)

Verlet Algorithm with periodic conditions

 

[BvU]

Hello @gibatom ,

$\qquad$!​

So what have you done so far to pinpoint the problem and to exclude programming errors ?
The description

the energy in the box grows all of a sudden.

is rather vague and our telepathic capabilities are severely limited

Also check out the PF guidelines: they 'ask' you to post your best attempt

$\$

 

[DrClaude]

The problem is that you are not actually implementing periodic boundary conditions. The latter requires your atom to interact not only with atoms on its side of the wall, but also with the image of the atoms "on the other side of wall," brought about by the PBC.

If you add images, you need to introduce a cut-off distance beyond which you do not consider atom-atom interactions to be significant anymore, otherwise you end up with an infinite number of pair interactions.

Adding the walls to the simulation is not that difficult and will give you more realistic results. (If at a given time step the atom would pass a wall, simply make it bounce with the proper angle.)

 

[gibatom]

Hello @gibatom ,

$\qquad$!​

So what have you done so far to pinpoint the problem and to exclude programming errors ?
The description

is rather vague and our telepathic capabilities are severely limited

Also check out the PF guidelines: they 'ask' you to post your best attempt

$\$

Thank you for your answer BvU. Yes i am new here, i am not still used to the functioning of the forum: i solemnly promise to read the PF guidelines
Sorry for the unclearness, i will give more details in my answers to DrClaude if you want to read it.
I already identified a problem with my simulation : it is the non conservation of energy and i see that it occurs when a particle leaves the box from one side, enters it from the other side and interacts with other atoms with a certain speed. I don't know how to adress this issue. Do you think i shoud post the code here ?

 

[gibatom]

The problem is that you are not actually implementing periodic boundary conditions. The latter requires your atom to interact not only with atoms on its side of the wall, but also with the image of the atoms "on the other side of wall," brought about by the PBC.

If you add images, you need to introduce a cut-off distance beyond which you do not consider atom-atom interactions to be significant anymore, otherwise you end up with an infinite number of pair interactions.

Adding the walls to the simulation is not that difficult and will give you more realistic results. (If at a given time step the atom would pass a wall, simply make it bounce with the proper angle.)

Thank you very much DrClaude for your answer.

Sorry, reading your answer, i realize i was not accurate enough describing my code. I effectively use the image of the atoms as you proposed, and calculate the interactions between those atoms and the ones on the box, only if the distance is less than 3 times the equilibrum distance of two atoms (using the minimization of the Lennard-Jones potential). I know that this number 3 is a specific number but i plan to make it evolve and sees how the gas evolve.
Therefore, at the end of the day, i make as if there were 8 artifical boxes surrounding my box, my actual box being in the middle and i use a cut off distance. (8 boxes are enough since the box is much bigger than the cut-off distance).

And i effectively plan to simulate the box using the rebounds on the walls but i need to also implement the periodic boundary conditions because i would like to do a comparaison. According to statistical thermodynamics the macroscopic values should be the same, considering that the number of atoms is great, with both methods.

By the way, do you think i should post my code, so you can have a look ?

 

[DrClaude]

By the way, do you think i should post my code, so you can have a look ?

Yes, please do, because if you are using images correctly, I don't see how the energy could jump when an atom crosses a wall.

 

[gibatom]

Yes, please do, because if you are using images correctly, I don't see how the energy could jump when an atom crosses a wall.

Here is my code. I translated almost everything into english but i let some variables in french if they are understandable (like "masse" for mass for example). I tried to make helpful commentaries as well to facilitate your reading of the code but if something is still unclear in my code, please let me know.

Thank you very much for the time you take to help me

 

[bob012345]

Is it possible to define a box without walls such that whenever an atom crosses the boundary of the imaginary box it just re-enters the opposite side like the box was the limits of a finite but unbounded universe? You could also make the box a sphere (or circle in 2D).

 

[berkeman]

Here is my code.

Where?

 

[gibatom]

Where?

I tried to attach it on the discussion but the forum won't publish it altough it seems to work when i try to upload it. Do you know how i can fix it ?

 

[gibatom]

Is it possible to define a box without walls such that whenever an atom crosses the boundary of the imaginary box it just re-enters the opposite side like the box was the limits of a finite but unbounded universe? You could also make the box a sphere (or circle in 2D).

Yes it should be possible if we consider that the atoms in the box interacts with fictive atoms around the box, who are images of the atoms in the box.

I also plan to modelize the rebounds on the box as you suggest and yes, indeed, working with a circle box is a good idea as it should help to avoid angle problems.

 

[berkeman]

I tried to attach it on the discussion but the forum won't publish it altough it seems to work when i try to upload it. Do you know how i can fix it ?

What is the file type? Can you attach it as a text file? You can also paste it into the forum window and enclose it with code tags like this:

[ code ]
<<your code>>
[ /code ]

(but leave out the spaces in the code tags)

 

[gibatom]

It's just a normal text file.
I am trying what you said :

from math import *
from numpy import *
from random import *
import matplotlib.pyplot as plt
N=4 #number of atoms in the gas

#using natural units
masse=1 #mass of an atom
dt=0.01 #time discretization
re = 2**(1/6) # equilibrum distance of the gas
length = 10 # length of the box
rc= 3 * re # cut_off distancedef aleatory_conditions (N): # return a random initial disposition and speeds
                             # for a gas of atoms
    l = []
    for i in range (N):
            
            a = uniform (0,length)
            b= uniform (0,length)
            c = uniform (-10,10) 
            l.append ([a,b,c,d])     
    return l

#Here is an initial configuration with which i saw the problem of the non-conservation
#of the energy : to launch the simulation with 50 images you can write simulation (a,50) 
#The numbers appearing before each simulation is the energy (defined below)

a= [[0.12380704056059177,
  0.6400585085502852,
  -6.7628327284249945,
  6.749532700539124],
 [3.812100505520707, 4.278142550643636, 9.01799159869551, 1.830313934178493],
 [4.079411610252755, 3.497598765109152, -9.0583500309088, 9.620341248119676],
 [7.6770984176024335,
  5.507714560653829,
  9.558704480961538,
  -7.114344763321643],
 [9.056707641094103, 8.099909642177849, 4.1853720606214, -0.2568527731476937],
 [7.0244623741647985,
  9.511569037976678,
  -0.2562828129841783,
  -5.262070353630726],
 [8.54483045652798, 4.711184438020645, 3.8934052821016767, -9.270864515219447],
 [1.015324476971885, 6.000336731042237, 8.909315796456141, -8.3327988101696]]

'''
a= aleatory_conditions (N)
'''#### FORCE CALCULATION ####

def force (r1,r2): #calculate the force one atom in r1 ([x1,x2])
                   #creates on a seconde atom in r2 ([x2,y2])
    x1,y1 = r1
    x2,y2 = r2 
    r = sqrt ((x1-x2)**2 + (y1-y2)**2)
    
    return - 4 * (-12*(1/r)**13 + 6*(1/r)**7)
    
def force_total (i,l): #from the list containing all the positions and speed of all atoms
                       # (x1,y1,v1,v2), it returns the total force felt by the atom i
                       # taking into account rc
    Fx=0                   
    Fy = 0
    xi = l[i][0]
    yi = l[i][1]
    for j in range (len(l)):
        if j != i:
            liste = [-1,0,1] 
            for a in liste :
             for b in liste : #the cases a=-1, a=1, (a=0,b=-1) and (a=1,b=1) are 
                              #counting the interactions with the atoms 
                              #"outside the box" 
                    xj = l[j][0] + a * length
                    yj = l[j][1] + b * length
                    r= sqrt ((xj-xi)**2 + (yj-yi)**2)
                    if xj != xi :
                     teta = atan (abs((yj-yi)/(xj-xi))) 
                    else :
                        teta = pi / 2
                    Ftot = force ([xi,yi],[xj,yj])
                    
                    if r < rc : 
                                
                        if a==-1 :
                           
                            Fx += Ftot * cos (teta)
                            
                            if b==-1 :
                                Fy += Ftot * sin (teta)
                                
                                
                                
                            elif b==0 :
                                if yj > yi :
                                     Fy -= Ftot * sin (teta)
                                     
                                else :
                                     Fy += Ftot * sin (teta)
                                     
                            else : # b=1
                                 Fy -= Ftot * sin (teta)
                                 
                        
                        elif a == 0 :
                            
                            if b == -1:
                                Fy += Ftot * sin(teta)
                                
                                if xj > xi :
                                    Fx -= Ftot * cos (teta)
                                    
                                else :
                                    Fx += Ftot * cos (teta)
                                   
                            if b == 0 :
                                
                                 if (xj > xi) :
                                  Fx -= Ftot*cos(teta)
                                  
                                 else  : 
                                  Fx += Ftot*cos(teta)
                                  
                                  
                                 if (yj <yi):
                                   Fy += Ftot*sin(teta)
                                 
                                 else : 
                                  Fy -= Ftot * sin (teta)
                                  
                            
                            if b==1 :
                                 Fy -= Ftot * sin (teta)
                                 
                                 if xj > xi :
                                    Fx -= Ftot * cos (teta)
                                    
                                 else :
                                    Fx += Ftot * cos (teta)
                                   
                        
                        else : #a=1
                         
                           Fx -= Ftot * cos (teta)
                           
                           if b==-1 :
                              Fy += Ftot * sin (teta)  
                              
                              
                           elif b==0 :
                               
                                if yj > yi :
                                    Fy -= Ftot * sin (teta) 
                                    
                                else :
                                    Fy += Ftot * sin (teta) 
                                  
                                  
                           else : #b=1
                             Fy -= Ftot * sin (teta)                      
    return (Fx,Fy)

#### ATOMS EVOLUTION ####

def incrementation (l): #gives the new list of atoms positions and speeds at t+dt
    l_p1 = []
    l_p2=[]
    N=len(l) #nombre de particules
    
    for i in range (N):
       
        x_t, y_t, v_x_t, v_y_t = l[i]
        fx, fy = force_total (i,l) 
        
        #we calculat new positions and speeds:
            
        x_t_plus_dt = mod (x_t + v_x_t * dt + (1/2)*(fx/masse)*dt**2, length)
        y_t_plus_dt = mod (y_t + v_y_t * dt + (1/2)*(fy/masse)*dt**2,length)
        #we work with mod (...,length) in order to confine the atoms
        v_x_t_plus_dt = (1/2) * (fx/masse) * dt + v_x_t
        v_y_t_plus_dt = (1/2) * (fy/masse) * dt + v_y_t
        l_p1.append ([ x_t_plus_dt,y_t_plus_dt, v_x_t_plus_dt, v_y_t_plus_dt])      
    for i in range (N):
        
        fx,fy = force_total (i,l_p1)
        #Verlet method :
            
        x_t, y_t, v_x_t, v_y_t = l_p1[i]
        v_x_t_plus_dt = (1/2)*(fx/masse) * dt + v_x_t
        v_y_t_plus_dt = (1/2)*(fy/masse) * dt + v_y_t
        l_p2.append ([ x_t,y_t, v_x_t_plus_dt, v_y_t_plus_dt])
    
    return l_p2    
    
      

#### SIMULATION PLOTTING ####def plotting (l): #plot the position of all atoms from a list of positions and speeds (given
                   #by aleatory conditions for example
                
    x = []
    y = []
    for i in range (len (l)):
        a = l[i][0]
        b = l[i][1]
        x.append(a)
        y.append(b)
    plt.scatter (x,y)
    plt.xlim(0,length)
    plt.ylim(0,length)
    plt.show()

   
def simulation (l,N): # plot how the gas evolve for a time N*dt starting with the 
                      #initial configuration l
    m=l
    for i in range (N):
        print (energy (m))
        plotting (m)
        m = incrementation (m)
    b= str(m)
    file = open("data.txt", "w") 
    file.write(b) #to stock the last configuration of atoms
    file.close()

        

  

    
#### ENERGY CALCULATION ####

def potential (r1,r2): #calculate the potential energy of two atoms in r1 et r2

    x1,y1 = r1
    x2,y2= r2 
    r = sqrt ((x1-x2)**2 + (y1-y2)**2)
    return 4 * ((1/r)**12-(1/r)**6)

def energy (l):
    Ec_totale = 0 #kinetic energy
    Ep_totale= 0 #potential energy
    for i in range (len(l)):
        x_t, y_t, v_x_t, v_y_t = l[i]
        Ec= (1/2)*masse*(v_x_t**2 + v_y_t**2)
        Ec_totale +=Ec
    for i in range (len(l)):
        for j in range (len(l)):
            if j !=i:
                x1,y1 = l[i][0], l[i][1]
                x2,y2 = l[j][0], l[j][1]
                Ep= potential ([x1,y1],[x2,y2])
                Ep_totale += Ep
    return  (Ec_totale + Ep_totale/2)

 

[gibatom]

What is the file type? Can you attach it as a text file? You can also paste it into the forum window and enclose it with code tags like this:

[ code ]
<<your code>>
[ /code ]

(but leave out the spaces in the code tags)

it worked, thank you !

 

[gibatom]

Here is my code. I translated almost everything into english but i let some variables in french if they are understandable (like "masse" for mass for example). I tried to make helpful commentaries as well to facilitate your reading of the code but if something is still unclear in my code, please let me know.

Thank you very much for the time you take to help me

Here is the missing code :

from math import *
from numpy import *
from random import *
import matplotlib.pyplot as plt
N=4 #number of atoms in the gas

#using natural units
masse=1 #mass of an atom
dt=0.01 #time discretization
re = 2**(1/6) # equilibrum distance of the gas
length = 10 # length of the box
rc= 3 * re # cut_off distancedef aleatory_conditions (N): # return a random initial disposition and speeds
                             # for a gas of atoms
    l = []
    for i in range (N):
            
            a = uniform (0,length)
            b= uniform (0,length)
            c = uniform (-10,10) 
            l.append ([a,b,c,d])     
    return l

#Here is an initial configuration with which i saw the problem of the non-conservation
#of the energy : to launch the simulation with 50 images you can write simulation (a,50) 
#The numbers appearing before each simulation is the energy (defined below)

a= [[0.12380704056059177,
  0.6400585085502852,
  -6.7628327284249945,
  6.749532700539124],
 [3.812100505520707, 4.278142550643636, 9.01799159869551, 1.830313934178493],
 [4.079411610252755, 3.497598765109152, -9.0583500309088, 9.620341248119676],
 [7.6770984176024335,
  5.507714560653829,
  9.558704480961538,
  -7.114344763321643],
 [9.056707641094103, 8.099909642177849, 4.1853720606214, -0.2568527731476937],
 [7.0244623741647985,
  9.511569037976678,
  -0.2562828129841783,
  -5.262070353630726],
 [8.54483045652798, 4.711184438020645, 3.8934052821016767, -9.270864515219447],
 [1.015324476971885, 6.000336731042237, 8.909315796456141, -8.3327988101696]]

'''
a= aleatory_conditions (N)
'''#### FORCE CALCULATION ####

def force (r1,r2): #calculate the force one atom in r1 ([x1,x2])
                   #creates on a seconde atom in r2 ([x2,y2])
    x1,y1 = r1
    x2,y2 = r2 
    r = sqrt ((x1-x2)**2 + (y1-y2)**2)
    
    return - 4 * (-12*(1/r)**13 + 6*(1/r)**7)
    
def force_total (i,l): #from the list containing all the positions and speed of all atoms
                       # (x1,y1,v1,v2), it returns the total force felt by the atom i
                       # taking into account rc
    Fx=0                   
    Fy = 0
    xi = l[i][0]
    yi = l[i][1]
    for j in range (len(l)):
        if j != i:
            liste = [-1,0,1] 
            for a in liste :
             for b in liste : #the cases a=-1, a=1, (a=0,b=-1) and (a=1,b=1) are 
                              #counting the interactions with the atoms 
                              #"outside the box" 
                    xj = l[j][0] + a * length
                    yj = l[j][1] + b * length
                    r= sqrt ((xj-xi)**2 + (yj-yi)**2)
                    if xj != xi :
                     teta = atan (abs((yj-yi)/(xj-xi))) 
                    else :
                        teta = pi / 2
                    Ftot = force ([xi,yi],[xj,yj])
                    
                    if r < rc : 
                                
                        if a==-1 :
                           
                            Fx += Ftot * cos (teta)
                            
                            if b==-1 :
                                Fy += Ftot * sin (teta)
                                
                                
                                
                            elif b==0 :
                                if yj > yi :
                                     Fy -= Ftot * sin (teta)
                                     
                                else :
                                     Fy += Ftot * sin (teta)
                                     
                            else : # b=1
                                 Fy -= Ftot * sin (teta)
                                 
                        
                        elif a == 0 :
                            
                            if b == -1:
                                Fy += Ftot * sin(teta)
                                
                                if xj > xi :
                                    Fx -= Ftot * cos (teta)
                                    
                                else :
                                    Fx += Ftot * cos (teta)
                                   
                            if b == 0 :
                                
                                 if (xj > xi) :
                                  Fx -= Ftot*cos(teta)
                                  
                                 else  : 
                                  Fx += Ftot*cos(teta)
                                  
                                  
                                 if (yj <yi):
                                   Fy += Ftot*sin(teta)
                                 
                                 else : 
                                  Fy -= Ftot * sin (teta)
                                  
                            
                            if b==1 :
                                 Fy -= Ftot * sin (teta)
                                 
                                 if xj > xi :
                                    Fx -= Ftot * cos (teta)
                                    
                                 else :
                                    Fx += Ftot * cos (teta)
                                   
                        
                        else : #a=1
                         
                           Fx -= Ftot * cos (teta)
                           
                           if b==-1 :
                              Fy += Ftot * sin (teta)  
                              
                              
                           elif b==0 :
                               
                                if yj > yi :
                                    Fy -= Ftot * sin (teta) 
                                    
                                else :
                                    Fy += Ftot * sin (teta) 
                                  
                                  
                           else : #b=1
                             Fy -= Ftot * sin (teta)                      
    return (Fx,Fy)

#### ATOMS EVOLUTION ####

def incrementation (l): #gives the new list of atoms positions and speeds at t+dt
    l_p1 = []
    l_p2=[]
    N=len(l) #nombre de particules
    
    for i in range (N):
       
        x_t, y_t, v_x_t, v_y_t = l[i]
        fx, fy = force_total (i,l) 
        
        #we calculat new positions and speeds:
            
        x_t_plus_dt = mod (x_t + v_x_t * dt + (1/2)*(fx/masse)*dt**2, length)
        y_t_plus_dt = mod (y_t + v_y_t * dt + (1/2)*(fy/masse)*dt**2,length)
        #we work with mod (...,length) in order to confine the atoms
        v_x_t_plus_dt = (1/2) * (fx/masse) * dt + v_x_t
        v_y_t_plus_dt = (1/2) * (fy/masse) * dt + v_y_t
        l_p1.append ([ x_t_plus_dt,y_t_plus_dt, v_x_t_plus_dt, v_y_t_plus_dt])      
    for i in range (N):
        
        fx,fy = force_total (i,l_p1)
        #Verlet method :
            
        x_t, y_t, v_x_t, v_y_t = l_p1[i]
        v_x_t_plus_dt = (1/2)*(fx/masse) * dt + v_x_t
        v_y_t_plus_dt = (1/2)*(fy/masse) * dt + v_y_t
        l_p2.append ([ x_t,y_t, v_x_t_plus_dt, v_y_t_plus_dt])
    
    return l_p2    
    
      

#### SIMULATION PLOTTING ####def plotting (l): #plot the position of all atoms from a list of positions and speeds (given
                   #by aleatory conditions for example
                
    x = []
    y = []
    for i in range (len (l)):
        a = l[i][0]
        b = l[i][1]
        x.append(a)
        y.append(b)
    plt.scatter (x,y)
    plt.xlim(0,length)
    plt.ylim(0,length)
    plt.show()

   
def simulation (l,N): # plot how the gas evolve for a time N*dt starting with the 
                      #initial configuration l
    m=l
    for i in range (N):
        print (energy (m))
        plotting (m)
        m = incrementation (m)
    b= str(m)
    file = open("data.txt", "w") 
    file.write(b) #to stock the last configuration of atoms
    file.close()

        

  

    
#### ENERGY CALCULATION ####

def potential (r1,r2): #calculate the potential energy of two atoms in r1 et r2

    x1,y1 = r1
    x2,y2= r2 
    r = sqrt ((x1-x2)**2 + (y1-y2)**2)
    return 4 * ((1/r)**12-(1/r)**6)

def energy (l):
    Ec_totale = 0 #kinetic energy
    Ep_totale= 0 #potential energy
    for i in range (len(l)):
        x_t, y_t, v_x_t, v_y_t = l[i]
        Ec= (1/2)*masse*(v_x_t**2 + v_y_t**2)
        Ec_totale +=Ec
    for i in range (len(l)):
        for j in range (len(l)):
            if j !=i:
                x1,y1 = l[i][0], l[i][1]
                x2,y2 = l[j][0], l[j][1]
                Ep= potential ([x1,y1],[x2,y2])
                Ep_totale += Ep
    return  (Ec_totale + Ep_totale/2)

 

[gibatom]

Hello,

Did anyone have the time to look at my code ?