#Lennard Jones gas, 2d, Monte Carlo. Units: Length=sigma, # Energy=Epsilon, Mass=Atomic mass. Periodic boundary conditions from __future__ import division, print_function from visual import * from visual.graph import * import numpy as np def forces(ux,uy,N,L): # gives accelerations, all particles have equal masses=1 rx=ux - ux.reshape(N,1) #matrix of x differences trimx=greater(rx,L/2) #locates differences >L/2 rx=rx-L*trimx # subtract L from those elements trimx=less(rx,-L/2) #locates differences < -L/2 rx=rx+L*trimx #add L to those elements ry=uy - uy.reshape(N,1) trimy=greater(ry,L/2) ry=ry-L*trimy trimy=less(ry,-L/2) ry=ry+L*trimy r2=rx*rx+ry*ry+eye(N,N) #add and subtract the identity matrix (avoid 1/0 for diagonal) rm2=1./r2-eye(N,N) trim=greater(rm2,0.160) #cuts off the forces at 2.5*sigma, i.e. 1/r**2 < 0.16 set to 0. rm2=rm2*trim rm8=rm2**4 rm14=rm2**7 fff=24.*(2.*rm14-rm8) ffx=rx*fff ffy=ry*fff fx = sum(ffx,axis=0) fy = sum(ffy,axis=0) return fx,fy def pe(ux,uy,N,L): #potential energy rx=ux - ux.reshape(N,1) #matrix of x differences trimx=greater(rx,L/2) #locate differences >L/2 rx=rx-L*trimx #subtract L from those elements trimx=less(rx,-L/2) #locate differences < -L/2 rx=rx+L*trimx #add L to those elements ry=uy - uy.reshape(N,1) trimy=greater(ry,L/2) ry=ry-L*trimy trimy=less(ry,-L/2) ry=ry+L*trimy r2=rx*rx+ry*ry+eye(N,N) rm2=1./r2-eye(N,N) trim=greater(rm2,0.160) #cuts off the potential at 2.5*sigma rm2=rm2*trim rm6=rm2**3 rm12=rm2**6 r6_12=4.*(rm12-rm6) pot=0.5*sum(r6_12) return pot #Intialize parameters N=64 #Number of particles L=15. #Edge of box sp=1.6 #Initial spacing tries=2000*N delta=.6 T=0.7 #setup print ('T ',T, 'tries ', tries,'delta ',delta) Lsp=int(L/sp) #initial spacing Nsqrt=ceil(sqrt(N)) Atoms=[] ux=zeros(N) uy=zeros(N) #Graphics setup r=0.5 kskip=10*N # display after 10 tries per atom scene = display(title="Lennard-Jones MC", x=30, y=0, foreground=color.black,background=color.yellow) #Movie square = curve(pos=[(-L/2-r,-L/2-r),(-L/2-r,L/2+r), (L/2+r,L/2+r),(L/2+r,-L/2-r),(-L/2-r,-L/2-r)]) #container g2=gdisplay(x=500,width=500,title="PE", foreground=color.black,background=color.white) #energies potential=gcurve(color=color.blue) for i in arange(N): #Initial conditions. x,y=divmod(i,Nsqrt) ux[i]=-L/4+sp*(x + y/2.) if ux[i] > L/2-r: ux[i]=-L+2*r+ux[i] uy[i]=-L/4+sp*y*0.866 c=color.blue Atoms=Atoms+[sphere(pos=(ux[i],uy[i],0), radius=r, color=c)] PE0=pe(ux,uy,N,L) print("PE(0)=",PE0) PE=PE0 #Start dynamics k=0 # counts acceptances virial=0. PEave=0. count=0 #counts averaging steps for t in arange(tries): rate(10000) j=random.randint(N) #try to move this atom dx=delta*(random.random()-0.5) xo=ux[j] #old x xn=xo+dx#new x if xn>L/2: xn=xn-L if xn<-L/2: xn=xn+L ux[j]=xn dy=delta*(random.random()-0.5) yo=uy[j] yn=yo+dy if yn>L/2: yn=yn-L if yn<-L/2: yn=yn+L uy[j]=yn PEn=pe(ux,uy,N,L) DE=PEn-PE w=exp(-DE/T) u=random.random() if u < w: #keep step PE=PEn k=k+1 else: #discard step ux[j]=xo uy[j]=yo if mod(t,kskip) ==0: #Plot every kskip steps, and average potential.plot(pos=(t,PE)) for i in arange(N): Atoms[i].pos=(ux[i],uy[i],0.) if t > tries/3: count=count+1 ax,ay=forces(ux,uy,N,L) virial=virial+0.5*(np.dot(ax,ux) +np.dot(ay,uy)) PEave=PEave+PE virial=virial/count PEave=PEave/count print('acceptances',k,'ratio',k/tries) print('Temp=KE/N= ',T,' /N= ',PEave/N ) print('Virial= ',virial) Pi=N*T/L**2 print('Ideal P= ', Pi) print('Pressure= ', Pi+virial/L**2)