0
$\begingroup$

I am simulating a Lennard-Jones fluid using MC simulation. The code always uses a reduced unit. I want to find the potential energy of the system. Periodic boundary condition implemented. I have simulated with 2 particles to find out if the simulation gives correct result. For 2 particles the equilibrium separation is $2^{1/6}$ , so energy is -1. (Take $\sigma =1$) . The code is below.

import numpy as np
import random
import matplotlib.pyplot as plt


# reduced units:
#T(reduced) = kT/epsilon | r(reduced) = r/sigma | U(reduced) = U/epsilon
# General Parameters
DIM=2
npart=2
L=10
volume=L**DIM
density=npart/volume
print("volume = ", volume, " density = ", density,"Number of atoms =",npart)
print ("L is " ,L)
T = 2; Nsteps = 50000; maxdr =0.0001;printfeq=100;DIM=2
#System parameters
beta=1/T



def E(dr2):
    """Returns LJ 6-12 interaction energy for a particular distance between
       2 particles"""

    return  4*((dr2)**(-6) - (dr2)**(-3)) # r is given in unit of sigma. dr2 is distance^2

def P(x):
    """ gives boltzman factor for position at x"""
    return np.exp(-(beta*E(x))) 

def PBC(L,pos):       

    """PBC check for dim dimension system with equal length L in all dimension. 
    INPUT : position array,length,dimension
    OUTPUT: New position
    """
    for k in range(DIM):
            if (pos[k]>0.5):
                pos[k]=pos[k]-1
            if (pos[k]<-0.5):
                pos[k]=pos[k]+1


    return (pos)


def distance(current_position):
    """Takes the current position array of the configuration and finds out 
    distance between each pairs. Neglectd if distance > rcutoff
    INPUT: Array of current position of each particles
    OUTPUT: Array containing distances between each pair of LJ particles.
    """
    Distances=[]
    for i in range (npart):
        for j in range (i+1,npart):
            dr=(current_position[i]-current_position[j])*L
            dr2=np.dot(dr,dr)

            if (dr2!=0):
                Distances.append(dr2)

    return Distances



Energy=[0 for _ in range (Nsteps)]
Distances=[0 for _ in range (Nsteps)]
current_position=np.zeros([npart,DIM])



#------------------ Initialise the Setup with particles distributed uniformly ------------

ip=-1
x=0
y=0
lim=int(np.sqrt(npart))+1
for i in range(0,lim):
    for j in range(0,lim):
        if(ip<npart):
            x=i*(1/lim)
            y=j*(1/lim)
            current_position[ip]=np.array([x,y])
            ip=ip+1
        else:
            break
MassCentre = np.sum(current_position,axis=0)/npart
current_position=current_position-MassCentre


Distances[0]=distance(current_position)



for i in Distances[0]:

    Energy[0]+=E(i)
print(Energy[0])


# -------------------------MC Simulation ----------------------------
rejected=0
for step in range(1,Nsteps):
    if (step%printfeq==0):
        print ("Completed ",step,"steps")
    trial_position=np.zeros([npart,DIM])

    trial_energy=0

    for i in range (npart):

        displacex=(random.uniform(0,1)-0.5)*maxdr 
        displacey=(random.uniform(0,1)-0.5)*maxdr 

        pos=current_position[i]+np.array([displacex,displacey])
        trial_position[i]=PBC(L,pos)
    Distances[step]=distance(trial_position)
    #print(trial_position)

    for i in ((Distances[step])):
        trial_energy+= E(i)
    if (trial_energy<Energy[step-1]):
        current_position=trial_position
        Energy[step]=trial_energy 
    else:
        delta=trial_energy-Energy[step-1]
        if (random.random()<P(delta)):

            current_position=trial_position
            Energy[step]=trial_energy 

        else:
            rejected+=1

            Energy[step]=Energy[step-1]



print(Energy)           
print ("Rejected moves ",rejected,"out of",Nsteps)
steps=np.arange(0,Nsteps)            
plt.figure(1)    
plt.plot(steps, Energy,'o',label='simulation result')
#plt.plot(steps, Distances,'b-',label='simulation result')
plt.xlabel(' steps')
plt.ylabel('Energy of configuration')
plt.show()


I have used 50000 steps with no separate steps for equilibrium. The output I get is like enter image description here

I want to know if this behavior is normal. It doesn't seem so. The sudden jump of potential energy to -1. Or am I doing something wrong? In MC I don't have info on kinetic energy, so I can't check energy conservation.Thank you in advance.

$\endgroup$
1
$\begingroup$

There are few issues in your code. The main one is that, in the P(x) function, return np.exp(-(beta*E(x))) should be return np.exp(-(beta*x)).

Then you should increase the maxdr value. I reckon a value of 0.01 should be enough. These two changes are enough to get a plot that looks realistic.

A couple of additional points:

  • In general, attempting to move many particles at once dramatically reduces the MC acceptance probability. Usually only single-particle moves are carried out.
  • You should apply PBC to distances, not only to the coordinates themselves, or particles that are close to opposite edges won't interact.
$\endgroup$
  • $\begingroup$ I did not get why E(x) should be x . It is the boltzman factor and it should contain the energy not displacement. $\endgroup$ – Endeavour Mar 28 at 17:28
  • 1
    $\begingroup$ Sure, but you call it like this: P(delta), where delta is an energy difference. $\endgroup$ – lr1985 Mar 28 at 17:33
  • $\begingroup$ Yeah,you are right. I had the P(delta) function for that purpose. $\endgroup$ – Endeavour Mar 29 at 2:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.