! ! Program for doing KMC simulation of a system with two ! kinds of particles: interstitials and vacancies, ! on a random lattice. ! ! Both are mobile, and the only reaction in the system ! is that when an i and a v comes to within r_recombine ! from each other, they recombine and vanish. ! ! N+N particles are generated initially randomly according to some ! distribution chosen with parameter mode with a characteristic length ! scale parameter rdistr ! ! Input: N mode rdistr Temperature rjump r_recombine tmax seed ! ! Defect jump rates are hardcoded in code ! ! For no recombination use r_recombine>=1e30 ! ! Note that this code is made specifically for educational purposes. ! Because of this some things are implemented in an as clear, rather ! than as efficient way as possible. ! module global implicit none ! module just used for generating global variables ! Can in C be replaced by ordinary global variables defined ! before main() integer :: seed type particledat double precision :: x,y,z integer :: type ! 0 vacancy, 1 interstitial double precision :: rate 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, gaussianrand end module global program intvackmc use global implicit none ! Variables read in from command line integer :: Norig,mode double precision :: rdistr,T,rjump,r_recombine,tmax ! Defect data type(particledat), allocatable :: particle(:) double precision, allocatable :: Ri(:) double precision :: Emi,Emv double precision :: wmi,wmv double precision :: ratei,ratev,Rtot,uR integer :: Nstepi,Nstepv,Nreactions,Nsteplastreaction ! Other variables integer :: i,j,Nstep,N,ilow,ihalf,ihigh double precision :: r,rsq,rsqi,rsqv,dx,dy,dz,theta,phi,time,tau,r_rec_sq character :: buf*80 integer, external :: iargc character*10 :: name(0:1) !integer, parameter :: Nstepnoreaction=10000000 ! End if no reaction in this many steps integer, parameter :: Nstepbailout=1000000000 ! Emergency bailout ! Set diffusivities. ! These values are for Si from Tang, Phys. Rev. B 55 (1997) 14279 ! w calculated from diffusivities of dv=1.18e-4 and di=1.58e-1 and r=2.35Å ! with w=6d/r^2 ! Emi=1.37; Emv=0.1; ! Units of eV wmi=1.717; wmv=0.001282; ! Units of 1/fs name(0)="V"; name(1)="I"; if (iargc() <8) then print *,'Usage: intvackmc N mode rdistr Temperature rjump r_recombine tmax(µs) seed' STOP endif call getarg(1,buf); read (unit=buf,fmt=*) Norig call getarg(2,buf); read (unit=buf,fmt=*) mode call getarg(3,buf); read (unit=buf,fmt=*) rdistr call getarg(4,buf); read (unit=buf,fmt=*) T call getarg(5,buf); read (unit=buf,fmt=*) rjump call getarg(6,buf); read (unit=buf,fmt=*) r_recombine call getarg(7,buf); read (unit=buf,fmt=*) tmax call getarg(8,buf); read (unit=buf,fmt=*) seed tmax=tmax*1e9; print *,'tmax',tmax ! Initialize random number generator call sgrnd(seed) ! Calculate rates: since T will not change, these stay constant as well ratei=wmi*exp(-Emi*e/(kB*T)); ratev=wmv*exp(-Emv*e/(kB*T)); print *,'Average jump rates [1/fs] for i, v:',ratei,ratev ! Original number of particles N=Norig*4; allocate(particle(N)); allocate(Ri(0:N)); ! Generate initial positions of particles according to mode j=1; do i=1,Norig if (mode==0 .or. mode==1 .or. mode==2) then ! Gaussian distributions with rdistr for i's thrice that of v ! Generate one v r=gaussianrand()*rdistr theta=acos(1-2.0*grnd()); phi=grnd()*2*pi particle(j)%x=r*sin(theta)*cos(phi) particle(j)%y=r*sin(theta)*sin(phi) particle(j)%z=r*cos(theta) particle(j)%type=0 j=j+1 if (mode==1 .or. mode==2) then ! and one i r=gaussianrand()*3.0*rdistr theta=acos(1-2.0*grnd()); phi=grnd()*2*pi particle(j)%x=r*sin(theta)*cos(phi) particle(j)%y=r*sin(theta)*sin(phi) particle(j)%z=r*cos(theta) particle(j)%type=1; j=j+1 if (mode==2) then ! and 2 i more r=gaussianrand()*3.0*rdistr theta=acos(1-2.0*grnd()); phi=grnd()*2*pi particle(j)%x=r*sin(theta)*cos(phi) particle(j)%y=r*sin(theta)*sin(phi) particle(j)%z=r*cos(theta) particle(j)%type=1; j=j+1 r=gaussianrand()*3.0*rdistr theta=acos(1-2.0*grnd()); phi=grnd()*2*pi particle(j)%x=r*sin(theta)*cos(phi) particle(j)%y=r*sin(theta)*sin(phi) particle(j)%z=r*cos(theta) particle(j)%type=1; j=j+1 endif endif else print *,'Unknown mode',mode STOP 'Error: unknown mode' endif enddo N=j-1; Norig=N; print *,'Doing kmc simulation of',N,'particles with mode rdistr',mode,rdistr print *,'rjump=',rjump,' T= ',T,' r_recombine',r_recombine print *,'seed',seed open(10,file='intvacmovie.xyz',status='unknown') ! --------------- Check for initial recombination ------------- r_rec_sq=r_recombine**2 Nreactions=0 if (r_recombine<1e30) then do i=1,N if (particle(i)%type < 0) cycle ! Ignore already recombined do j=1,N if (particle(j)%type < 0) cycle ! Ignore already recombined ! Only check for opposite types if (particle(i)%type == particle(j)%type) cycle dx=particle(i)%x-particle(j)%x dy=particle(i)%y-particle(j)%y dz=particle(i)%z-particle(j)%z rsq=dx*dx+dy*dy+dz*dz; if (rsq <= r_rec_sq) then ! Mark pair for removal, will be ignored in remainder ! of run !print *,'Recombine',i,j particle(i)%type=-1 particle(j)%type=-1 Nreactions=Nreactions+1 ! Since i has recombined, ignore rest of j atoms exit endif enddo enddo ! Remove recombined particles from list i=1 do if (i>N) exit !print *,'test',i,particle(i)%type,N if (particle(i)%type < 0) then ! Move everything down one step do j=i,N-1 particle(j)%x=particle(j+1)%x particle(j)%y=particle(j+1)%y particle(j)%z=particle(j+1)%z particle(j)%type=particle(j+1)%type particle(j)%rate=particle(j+1)%rate enddo N=N-1 i=i-1 ! Need to handle same index again !print *,i,j,N endif i=i+1 enddo endif print *,' Initally recombined',Nreactions,'defect pairs,',N,' particles remain' ! Main loop time=0.0; Nstep=0; Nstepi=0; Nstepv=0; Nsteplastreaction=0 do if (N==0) exit !if (Nstep > Nsteplastreaction+Nstepnoreaction) exit if (time > tmax) exit ! ------- Determine rates of all possible events --------- do i=1,N if (particle(i)%type == 0) particle(i)%rate=ratev if (particle(i)%type == 1) particle(i)%rate=ratei enddo ! Form cumulative distribution Ri function of rates ri Rtot=0.0; Ri(0)=0.0; do i=1,N Rtot=Rtot+particle(i)%rate Ri(i)=Rtot enddo !print *,(Ri(i),i=1,N),Rtot ! ------------------ Select an event ------------------- uR=grnd()*Rtot; ! Find where this element is by binary search in Ri ! Note that it is important to have Ri(0) as well ilow=0; ihigh=N; do if (ihigh <= ilow+1) exit ihalf=int((ihigh+ilow)/2+0.5) if (Ri(ihalf) > uR) then ihigh=ihalf else ilow=ihalf endif enddo i=ihigh if (i<1 .or. i>N) STOP 'i HORROR ERROR' !print *,'Chose',i if (mod(Nstep,10000)==0) then ! Output movie every now and then write(10,*) N write(10,*) 'Defect positions',Nstep,time,' fs' do j=1,N write(10,*) name(particle(j)%type),particle(j)%x,particle(j)%y,particle(j)%z,particle(j)%type enddo call flush(10) endif ! --------------- Carry out event on i -------------------- ! i.e. in this case move defect by rjump theta=acos(1-2.0*grnd()); phi=grnd()*2*pi particle(i)%x=particle(i)%x+rjump*sin(theta)*cos(phi) particle(i)%y=particle(i)%y+rjump*sin(theta)*sin(phi) particle(i)%z=particle(i)%z+rjump*cos(theta) if (particle(i)%type == 0) then Nstepv=Nstepv+1 else Nstepi=Nstepi+1 endif ! ------------------ Check for possible recombination of i ---------- if (r_recombine<1e30) then do j=1,N if (particle(j)%type < 0) cycle ! Ignore already recombined ! Only check for opposite types if (particle(i)%type == particle(j)%type) cycle dx=particle(i)%x-particle(j)%x dy=particle(i)%y-particle(j)%y dz=particle(i)%z-particle(j)%z rsq=dx*dx+dy*dy+dz*dz; if (rsq <= r_rec_sq) then ! Mark pair for removal, will be ignored in remainder ! of run !print *,'Recombine',i,j,Nreactions,N particle(i)%type=-1 particle(j)%type=-1 Nreactions=Nreactions+1 Nsteplastreaction=Nstep ! Since i has recombined, ignore rest of j atoms exit endif enddo ! Note that after this i index may no longer point to recombined particle! ! Remove recombined particles from list i=1 do if (i>N) exit !print *,'test',i,particle(i)%type,N if (particle(i)%type < 0) then ! Move everything down one step do j=i,N-1 particle(j)%x=particle(j+1)%x particle(j)%y=particle(j+1)%y particle(j)%z=particle(j+1)%z particle(j)%type=particle(j+1)%type particle(j)%rate=particle(j+1)%rate enddo N=N-1 i=i-1 ! Need to handle same index again !print *,i,j,N endif i=i+1 enddo endif ! ------------------ Increment time -------------------- tau= -log(grnd())/Rtot; time=time+tau Nstep=Nstep+1 if (mod(Nstep,10000)==0) then print *,'Nstep N tau time(µs)',Nstep,N,tau,time/1e9 endif if (Nstep == Nstepbailout) exit enddo print *,' ' write(10,*) N write(10,*) 'Defect positions',Nstep,time,' fs' do i=1,N write(10,*) name(particle(i)%type),particle(i)%x,particle(i)%y,particle(i)%z,particle(i)%type enddo close(10) rsq=0.0; rsqv=0.0; rsqi=0.0; do i=1,N rsq=rsq+particle(i)%x**2+particle(i)%y**2+particle(i)%z**2 if (particle(i)%type==0) then rsqv=rsqv+particle(i)%x**2+particle(i)%y**2+particle(i)%z**2 else rsqi=rsqi+particle(i)%x**2+particle(i)%y**2+particle(i)%z**2 endif enddo if (Nstep >= Nstepbailout) then print *,'Emergency bailout after',Nstep,' steps' endif if (N==0) then print *,'Simulation ended with complete recombination' endif !if (Nstep > Nsteplastreaction+Nstepnoreaction) then ! print *,'Simulation ended because no reactions observed over',Nstepnoreaction,' steps' !endif if (time >= tmax) then print *,'Simulation ended when maximum time reached',time/1e9,tmax endif print *,'Simulation finished at time',time/1e9,' µs after', Nstep,' steps' print *,'Number of steps for i v',Nstepi,Nstepv,' ratio',(1.0*Nstepi/max(Nstepv,1)) print *,'Recombined',Nreactions,'defect pairs,',N,' particles remain, fraction',1.0*N/Norig print *,'sqrt(rsq) sqrt(rsqv) sqrt(rsqi)',sqrt(rsq),sqrt(rsqv),sqrt(rsqi) end program intvackmc double precision function gaussianrand() implicit none ! --------------------------------------------------------- ! Random numbers with normal (Gaussian) distribution. ! Mean is 0 and standard deviation is 1 ! See W.H.Press et al., Numerical Recipes 1st ed., page 203 ! --------------------------------------------------------- double precision :: fac,v1,v2,r double precision, save :: gset integer, save :: iset=0 double precision, external :: grnd if (iset.eq.0) then 1 v1 = 2.*grnd()-1. v2 = 2.*grnd()-1. r = v1*v1+v2*v2 if (r.ge.1.) goto 1 fac = sqrt(-2.0*log(r)/r) gset = v1*fac gaussianrand = v2*fac iset = 1 else gaussianrand = gset iset = 0 endif return end ! 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