! ! Metropolis MC code written by Kai Nordlund Feb 20-21, 2002 ! ! To compile: f90 -O -o metropolismc metropolismc.f90 -lU77 -lfio ! module global implicit none integer :: seed type particledat double precision :: x,y,z end type particledat double precision, parameter :: pi=3.14159265358979 double precision, parameter :: e=1.6021892d-19 double precision, parameter :: kB=1.380662d-23 double precision, external :: grnd end module global program metropolismc use global implicit none ! Variables read in from command line integer :: N,mode double precision :: T,size,djump ! Atom variables type(particledat), allocatable :: atom(:) ! Other variables integer :: i,j,k,l,iat,istep,isteplastsuccess,nfail,nsuccess double precision :: Eorig,Etot,Eiold,Einew,dx,dy,dz,a character*80 :: buf integer, external :: iargc integer, parameter :: Nfailmax=100000 integer, parameter :: istepmax=1000000 type(particledat) :: basis(4) logical :: accept ! Variables for physical results double precision :: Esum,Eave integer, parameter :: istepcollect=istepmax/2 if (iargc() <6) then print *,'Usage: intvackmc N size mode T djump seed' STOP endif call getarg(1,buf); read (unit=buf,fmt=*) N call getarg(2,buf); read (unit=buf,fmt=*) size call getarg(3,buf); read (unit=buf,fmt=*) mode call getarg(4,buf); read (unit=buf,fmt=*) T call getarg(5,buf); read (unit=buf,fmt=*) djump call getarg(6,buf); read (unit=buf,fmt=*) seed ! Initialize random number generator call sgrnd(seed) allocate(atom(N)) if (mode==0) then ! Random atom positions in range -size/2 to size/2 do i=1,N atom(i)%x = (grnd()-0.5)*size atom(i)%y = (grnd()-0.5)*size atom(i)%z = (grnd()-0.5)*size enddo print *,'Generated',N,' random atom positions' else if (mode==1) then N=256 a=3.62538 basis(1)%x=0.0; basis(1)%y=0.0; basis(1)%z=0.0; basis(2)%x=0.5; basis(2)%y=0.5; basis(2)%z=0.0; basis(3)%x=0.5; basis(3)%y=0.0; basis(3)%z=0.5; basis(4)%x=0.0; basis(4)%y=0.5; basis(4)%z=0.5; iat=0; do i=0,3 do j=0,3 do k=0,3 do l=1,4 iat=iat+1; atom(iat)%x=(i+basis(l)%x+0.25)*a-size/2; atom(iat)%y=(j+basis(l)%y+0.25)*a-size/2; atom(iat)%z=(k+basis(l)%z+0.25)*a-size/2; enddo enddo enddo enddo print *,'generated 4x4x4 FCC lattice with 256 atoms' else STOP 'Unknown mode' endif call energy(atom,N,size,Etot,0) Eorig=Etot print *,'Initial total energy',Etot open(10,file='mcmovie.xyz',status='unknown') istep=0; nsuccess=0; nfail=0; isteplastsuccess=0; Esum=0.0; Eave=0.0; do ! Trial jump ! Select atom i=int(grnd()*N)+1; if (i>N) i=N; ! Get old energy for atom i call energy(atom,N,size,Eiold,i) ! Move it randomly dx=(1.0-2.0*grnd())*djump; dy=(1.0-2.0*grnd())*djump; dz=(1.0-2.0*grnd())*djump; atom(i)%x=atom(i)%x+dx atom(i)%y=atom(i)%y+dy atom(i)%z=atom(i)%z+dz if (atom(i)%x >= size/2) atom(i)%x=atom(i)%x-size if (atom(i)%x < -size/2) atom(i)%x=atom(i)%x+size if (atom(i)%y >= size/2) atom(i)%y=atom(i)%y-size if (atom(i)%y < -size/2) atom(i)%y=atom(i)%y+size if (atom(i)%z >= size/2) atom(i)%z=atom(i)%z-size if (atom(i)%z < -size/2) atom(i)%z=atom(i)%z+size ! Get new energy for atom i call energy(atom,N,size,Einew,i) accept=.false. if (Einew <= Eiold) then accept=.true. else if (grnd() < exp(-(Einew-Eiold)*e/(kB*T))) accept=.true. endif if (accept) then isteplastsuccess=istep nsuccess=nsuccess+1 ! Follow evolution of total energy Etot=Etot-Eiold+Einew else atom(i)%x=atom(i)%x-dx atom(i)%y=atom(i)%y-dy atom(i)%z=atom(i)%z-dz if (atom(i)%x >= size/2) atom(i)%x=atom(i)%x-size if (atom(i)%x < -size/2) atom(i)%x=atom(i)%x+size if (atom(i)%y >= size/2) atom(i)%y=atom(i)%y-size if (atom(i)%y < -size/2) atom(i)%y=atom(i)%y+size if (atom(i)%z >= size/2) atom(i)%z=atom(i)%z-size if (atom(i)%z < -size/2) atom(i)%z=atom(i)%z+size nfail=nfail+1 endif istep=istep+1 ! Collect physical results if (istep>istepcollect) then Esum=Esum+Etot Eave=Esum/(N*(istep-istepcollect)) endif ! Print out results every now and then if (mod(istep,10000)==1) then ! Get total energy from scratch to prevent numerical ! drift call energy(atom,N,size,Etot,0) print *,"E",istep,1.0*nsuccess/istep,Etot/N,Eave write(10,*) N write(10,*) 'Atom positions',istep do i=1,N write(10,*) "Cu",atom(i)%x,atom(i)%y,atom(i)%z,1 enddo call flush(10) endif if (istep>= isteplastsuccess+Nfailmax) exit if (istep>= istepmax) exit enddo close(10) end program metropolismc subroutine energy(atom,N,size,Etot,iat) use global implicit none ! iat=0: calculate Etot ! else calclate Ei type(particledat) :: atom(*) integer :: N,iat double precision :: size,Etot integer :: i,j,j0 double precision :: dx,dy,dz,r,rsq,rcut,rcutsq,V,help1,help2 double precision :: morsealpha,morseDe,morseRe ! Cu params morseDe = 0.3429 morsealpha = 1.3588 morseRe = 2.866 rcut=7.0; rcutsq=rcut*rcut; Etot=0.0 i=0 ! Simplistic N^2 energy calculation. Could be much optimized do i=i+1 if (i>N) exit j0=i+1 if (iat>0) then i=iat j0=1 endif do j=j0,N if (i==j) cycle dx=atom(i)%x-atom(j)%x if (dx < -size/2) dx=dx+size; if (dx >= size/2) dx=dx-size; dy=atom(i)%y-atom(j)%y if (dy < -size/2) dy=dy+size; if (dy >= size/2) dy=dy-size; dz=atom(i)%z-atom(j)%z if (dz < -size/2) dz=dz+size; if (dz >= size/2) dz=dz-size; rsq=dx*dx+dy*dy+dz*dz if (rsq <= rcutsq) then ! i-j interaction within cutoff, evaluate it r=sqrt(rsq) ! Utilize the fact that one of the Morse terms is the ! square of the other ! help1=exp(-morsealpha*(r-morseRe)); help2=help1*help1 ! Calculate potential energy V=morseDe*(help2-2.0*help1) Etot=Etot+V endif enddo if (iat>0) exit enddo end subroutine energy ! Mersenne twister random number generator !************************************************************************ subroutine sgrnd(seed) ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) ! dimension mt(0:N-1) ! the array for the state particle common /block/mti,mt save /block/ ! ! setting initial seeds to mt[N] using ! the generator Line 25 of Table 1 in ! [KNUTH 1981, The Art of Computer Programming ! Vol. 2 (2nd Ed.), pp102] ! mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue ! return end !************************************************************************ double precision function grnd() ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) ! constant vector a parameter(UMASK = -2147483648) ! most significant w-r bits parameter(LMASK = 2147483647) ! least significant r bits ! Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) ! dimension mt(0:N-1) ! the array for the state vector common /block/mti,mt save /block/ data mti/N1/ ! mti==N+1 means mt[N] is not initialized ! dimension mag01(0:1) data mag01/0, MATA/ save mag01 ! mag01(x) = x * MATA for x=0,1 ! TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) ! if(mti.ge.N) then ! generate N words at one time if(mti.eq.N+1) then ! if sgrnd() has not been called, call sgrnd(4357) ! a default initial seed is used endif ! do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif ! y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) ! if(y.lt.0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif ! return end