#Lennard Jones gas, 2d, Monte Carlo. Units: Length=sigma, # Energy=Epsilon, Mass=Atomic mass. Periodic boundary conditions # Gibbs ensemble; method of Frenkel & Smit from __future__ import division, print_function from visual import * from visual.graph import * from random import random def pe(uxk,uyk,Nk,Lk): #potential energy rx=uxk - uxk.reshape(Nk,1) #matrix of x differences trimx=greater(rx,Lk/2) #locate differences >L/2 rx=rx-Lk*trimx #subtract L from those elements trimx=less(rx,-Lk/2) #locate differences < -L/2 rx=rx+Lk*trimx #add L to those elements ry=uyk - uyk.reshape(Nk,1) trimy=greater(ry,Lk/2) ry=ry-Lk*trimy trimy=less(ry,-Lk/2) ry=ry+Lk*trimy r2=rx*rx+ry*ry+eye(Nk,Nk) rm2=1./r2-eye(Nk,Nk) trim=greater(rm2,0.160) #cuts off the potential at 2.5*sigma rm2=rm2*trim rm6=rm2**3 - .0040906 #make pe continuous at cutoff rm12=rm2**6 - .0000167 r6_12=4.*(rm12-rm6) pot=0.5*sum(r6_12) return pot def move(j,uxk,uyk,Nk,Lk,PE,acc): #move particle j, reset ux,uy,PE,acc dx=delta*(random()-0.5) xo=uxk[j] #old x xn=xo+dx #new x if xn>Lk/2: xn=xn-Lk if xn<-Lk/2: xn=xn+Lk uxk[j]=xn dy=delta*(random()-0.5) yo=uyk[j] yn=yo+dy if yn>Lk/2: yn=yn-Lk if yn<-Lk/2: yn=yn+Lk uyk[j]=yn PEn=pe(uxk,uyk,Nk,Lk) DE=PEn-PE w=exp(-DE/T) u=random() if u < w: #keep step PE=PEn acc=acc+1 else: #discard step uxk[j]=xo uyk[j]=yo return uxk,uyk,PE,acc def expand(ux0,uy0,N0,L0,PE0,ux1,uy1,N1,L1,PE1,acc): #change volumes, reset V0,V1, ux, uy, L,PE,acc V0=L0**2 V1=V-V0 lnVn=log(V0/V1)+(random()-0.5)*dlnV #random walk in ln(V) V0n=V*exp(lnVn)/(1.+ exp(lnVn)) #new vol 0 V1n=V-V0n #new vol 1 L0n=V0n**(0.5) #new box 0 length L1n=V1n**(0.5) #new box 1 length f0=L0n/L0 f1=L1n/L1 ux0n=f0*ux0; uy0n=f0*uy0 #rescale positions 0 ux1n=f1*ux1; uy1n=f1*uy1 #rescale positions 1 PE0n=pe(ux0n,uy0n,N0,L0n) PE1n=pe(ux1n,uy1n,N1,L1n) DE0=PE0n-PE0-(N0+1)*T*(log(V0n/V0)) DE1=PE1n-PE1-(N1+1)*T*(log(V1n/V1)) w=exp(-DE0/T)*exp(-DE1/T) u=random() if u < w: #keep step PE0=PE0n PE1=PE1n ux0=ux0n ; uy0=uy0n ux1=ux1n ; uy1=uy1n L0=L0n L1=L1n acc=acc+1 return ux0,uy0,L0,PE0,ux1,uy1,L1,PE1,acc def exchange(ux0,uy0,N0,L0,PE0,ux1,uy1,N1,L1,PE1,acc): #move between boxes, reset ux,uy,N0,N1,PE,acc V0=L0**2 V1=V-V0 if random() < 0.5: #add to 0 prefactor=(N1*V0)/((N0+1)*V1) N0n=N0+1 N1n=N1-1 x=(random() -0.5)*L0; y= (random() -0.5)*L0 ux0n=append(ux0,x); uy0n=append(uy0,y) #add at random in0 j=int(N1*random()) ux1n=delete(ux1,j); uy1n=delete(uy1,j) # delete j box 1 else: #add to 1 prefactor=(N0*V1)/((N1+1)*V0) N1n=N1+1 N0n=N0-1 x=(random() -0.5)*L1; y= (random() -0.5)*L1 ux1n=append(ux1,x); uy1n=append(uy1,y) #add at random in 1 j=int(N0*random()) ux0n=delete(ux0,j); uy0n=delete(uy0,j) # delete j from 0 PE0n=pe(ux0n,uy0n,N0n,L0) PE1n=pe(ux1n,uy1n,N1n,L1) DE0=PE0n-PE0 DE1=PE1n-PE1 w=prefactor*exp(-DE0/T)*exp(-DE1/T) u=random() if u < w: #keep step PE0=PE0n PE1=PE1n ux0=ux0n ; uy0=uy0n ux1=ux1n ; uy1=uy1n N0=N0n N1=N1n acc=acc+1 return ux0,uy0,N0,PE0,ux1,uy1,N1,PE1,acc #Parameters N=80 #Total number of particles N0=int(N/2) #initial numbers N1=N-N0 Nswap=N L0=10. # Initial Edge of boxes L1=25. V=L0**2+L1**2 #Total volume T=0.55 #Temperature sp=1.6 #Initial spacing tries=9000*N delta=.6 #max particle move dlnV=0.5 #max change in ln(V), i.e. V' = V exp(rnd*dlnV) r=0.5 #display radius of particle diam=2*r density0=N0/L0/L0 #initial densities density1=N1/L1/L1 print ('T=',T) print(' initial density0=',density0,' density1=',density1) print(' tries/N=', tries/N,' delta=',delta) #Graphics setup s0 = display(title="Lennard-Jones MC, box0", x=30, y=0, foreground=color.black,background=color.white) #Movie0 mybox0=box(pos=(0,0,-1),height=L0+diam,width=.01,length=L0+diam, color=color.yellow) s1= display(title="Lennard-Jones MC, box1", x=500, y=0, foreground=color.black,background=color.white) #Movie1 mybox1=box(pos=(0,0,-1),height=L1+diam,width=.01,length=L1+diam, color=color.yellow) g2=gdisplay(x=500,y=500,width=500,title="N's, 0=blue,1=green", foreground=color.black,background=color.white) #numbers number0=gcurve(color=color.blue) number1=gcurve(color=color.green) g3 = gdisplay(y=500,width=500, title="densities, 0=blue,1=green", foreground=color.black,background=color.white) #density dens0=gcurve(color=color.blue) dens1=gcurve(colot=color.green) kskip=10*N # display after 10 tries per atom #Initial conditions N0s=ceil(sqrt(N0)) N1s=ceil(sqrt(N1)) Atoms0=[] Atoms1=[] ux0=zeros(N0); uy0=zeros(N0) ux1=zeros(N1); uy1=zeros(N1) for i in arange(N0): x,y=divmod(i,N0s) ux0[i]=-L0/4+sp*(x + y/2.) if ux0[i] > L0/2: ux0[i]=-L0+ux0[i] uy0[i]=-L0/4+sp*y*0.866 for i in arange(N1): x,y=divmod(i,N1s) ux1[i]=-L1/4+sp*(x + y/2.) if ux1[i] > L1/2: ux1[i]=-L1+ux1[i] uy1[i]=-L1/4+sp*y*0.866 PE0=pe(ux0,uy0,N0,L0) PE1=pe(ux1,uy1,N1,L1) print("initial PE0=",PE0, " PE1=",PE1) PE=PE0 #Start dynamics acc=0 # counts acceptances PE0ave=0. PE1ave=0. dens0ave=0. dens1ave=0. count=0 #counts averaging steps for t in arange(tries): rate(10000) j=int((N0+N1+Nswap)*random()) if j N0old: for i in arange(N0old,N0): Atoms0=Atoms0+[sphere(pos=(ux0[i],uy0[i],0), radius=r, color=color.blue)] elif N0 N1old: for i in arange(N1old,N1): Atoms1=Atoms1+[sphere(pos=(ux1[i],uy1[i],0), radius=r, color=color.blue)] elif N1 tries/2: count=count+1 PE0ave=PE0ave+PE0 PE1ave=PE1ave+PE1 dens0ave=dens0ave+density0 dens1ave=dens1ave+density1 PE0ave=PE0ave/count dens0ave=dens0ave/count dens1ave=dens1ave/count print('acceptances',acc,'ratio',acc/tries) print('Temp= ',T, 'final =', dens0ave,' =',dens1ave ) print('final N0=',N0,' N1=',N1) print('final PE0=', PE0,' PE1=',PE1)