Atomistic Simulations

1. Atomistic Simulations Byeong-Joo Lee Dept. of MSE Pohang University of Science and Technology (POSTECH) calphad@postech.ac.kr

2. What is Monte Carlo ? y implicit integer (i-n) implicit double precision (a-h,o-z) call random_seed i_circle = 0 do i = 1, 100000000 call random_number(x) call random_number(y) r2 = x*x + y*y if(r2 .lt. 1.d0) i_circle = i_circle + 1 if(mod(i,1000000) .eq. 0) then pi = 4.0d0 * dble(i_circle) / dble(i) write(*,*) i, pi endif enddo c stop end 1 x

3. Monte Carlo – Fundamentals (from Lecture Note of Prof. V. Vitek, 2002)

4. Monte Carlo – Distribution of states in real systems f(x) a b x

5. Monte Carlo – General Monte Carlo Algorithm

6. Monte Carlo – Canonical Ensemble (N,V,T)

7. Monte Carlo – Metropolis Method

8. Monte Carlo – Examples of Metropolis Method

9. Monte Carlo – Magnetization vs. Temperature by Ising model of spins

10. Monte Carlo Simulation - General Scheme (Potts Model: Srolovitz)

11. Abnormal Grain Growth– Computer Simulation

12. Monte Carlo – Isothermal-Isobaric Canonical Ensemble (N,P,T)

13. Monte Carlo –Canonical ensemble without exchanging atoms do j = 1, Natom call random_number(random) ; k = (Natom-Nsleep)*random + 1 + Nsleep oldpos(:,1) = pos(:,k); OldBoxSize = BoxSize call random_number(random) if(Natom*random.ge.3) then call MC_List(1,k,0,NListMC) ** Select neighbor atoms of the k atom call Compute_EAM_Energy(Ipbc,ene_old,NListMC) call random_number(random); deltaR(1,k) = delpos*(random-0.5d0) call random_number(random); deltaR(2,k) = delpos*(random-0.5d0) call random_number(random); deltaR(3,k) = delpos*(random-0.5d0) pos(:,k) = pos(:,k)+deltaR(:,k)/BoxSize DisplaceList(:,k) = DisplaceList(:,k) + deltaR(:,k) call Compute_EAM_Energy(Ipbc,ene_new,NListMC) else call Compute_EAM_Energy(Ipbc,ene_old,0) c Box size variation: First select one axis then select amount call random_number(random); ix = idint(3.d0*random) + 1 call random_number(random) BoxSize(ix) = BoxSize(ix) * (1.0d0 + dellen*(random-0.5d0)) call Compute_EAM_Energy(Ipbc,ene_new,0) endif

14. Monte Carlo –Canonical ensemble without exchanging atoms iaccept = 1 if(ene_new .gt. ene_old) then P01 = dexp((ene_old-ene_new)/boltz/TRequested) call random_number(random) if(random.gt.P01) then ** Return to the original configuration pos(:,k) = oldpos(:,1); BoxSize = OldBoxSize iaccept = 0 DisplaceList(:,k) = DisplaceList(:,k) + deltaR(:,k) endif endif if(iaccept.eq.1) then Dk = dsqrt(dot_product(DisplaceList(:,k),DisplaceList(:,k))) if( Dk > Skin ) ListUpdateRequested=1 if(ListUpdateRequested.eq.1) then call Update_List(Range,delR,Lg_fun,0,Nequi,Nsave,Ipbc) ListUpdateRequested = 0 endif delta_ene = ene_new - ene_old ene_sum = ene_sum + delta_ene endif enddo

15. Monte Carlo –Canonical ensemble exchanging atoms do j = 1, Natom call random_number(random); k = (Natom-Nsleep)*random + 1 + Nsleep call random_number(random); l = (Natom-Nsleep)*random + 1 + Nsleep itemp1 = IAtomType(k); itemp2 = IAtomType(l) oldpos(:,1) = pos(:,k); oldpos(:,2) = pos(:,l) OldBoxSize = BoxSize call random_number(random) if(Natom*random.ge.3) then call MC_List(2,k,l,NListMC) ** Select neighbor atoms of the k and l atoms call Compute_EAM_Energy(Ipbc,ene_old,NListMC) IAtomType(k) = itemp2; IAtomType(l) = itemp1 call random_number(random); deltaR(1,k) = delpos*(random-0.5d0) call random_number(random); deltaR(2,k) = delpos*(random-0.5d0) call random_number(random); deltaR(3,k) = delpos*(random-0.5d0) pos(:,k) = pos(:,k)+deltaR(:,k)/BoxSize; DisplaceList(:,k) = DisplaceList(:,k) + deltaR(:,k) call random_number(random); deltaR(1,l) = delpos*(random-0.5d0) …… …… pos(:,l) = pos(:,l)+deltaR(:,l)/BoxSize; DisplaceList(:,l) = DisplaceList(:,l) + deltaR(:,l) call Compute_EAM_Energy(Ipbc,ene_new,NListMC) else

16. Phase Separation and Surface Segregation in Cu-Ni Thin Films 373 K 298 K 473 K 800 K

17. Order/Disorder in Co-Pt System - S.I. Park et al., Scripta Mater., 2001. Property Pt3Co PtCo PtCo3 Cohesive Energy 5.500 5.215 4.873 (eV/atom) 5.555±0.017 5.228±0.005 Lattice Constant a=3.833 3.754, c/a=.98 3.625 (Å) a=3.831 3.745, c/a=.98 3.668 Transition 1070-1080 970-980 760-770 Temperature (K) 1000 1100 840

18. Effect of Interface Energy Anisotropy on Morphological Evolution Pure W W + 0.4wt% Ni Vacuum Annealing Pure W W-14at%Ni

19. System in contact with a heat reservoir – Canonical ensemble Number of NVT systems with Ei: Average Energy of a system:

20. Open System with a heat reservoir – Grand Canonical ensemble Number of μVT systems with Ei & Nr: Average Energy of a system: Average # of ptl. of a system:

21. Monte Carlo –Grand Canonical Ensemble (μ,V,T) • The number of particles is variable • fluctuations of the concentration are allowed • Generally, select randomly one of the procedures • Change of the configuration • exp(-ΔEp/kT) • Creation of a particle • exp[-(ΔEp-μ)/kT] • Destruction of particles • exp[-(ΔEp+μ)/kT]

22. Monte Carlo –Modified Grand Canonical Ensemble (μ,V,T)

23. Monte Carlo –Grand Canonical ensemble do j = 1, Natom call random_number(random); k = Natom*random+1 itemp0 = IatomType(k); oldpos(:,1) = pos(:,k); OldBoxSize = BoxSize call random_number(random) if(Natom*random.ge.3) then call MC_List(1,k,0,NListMC); call Compute_EAM_Energy(Ipbc,ene_old,NListMC) if(random.gt.0.5) then; itemp1 = 2 else; itemp1 = 1; endif IatomType(k) = itemp1 call random_number(random); deltaR(1,k) = delpos*(random-0.5d0) call random_number(random); deltaR(2,k) = delpos*(random-0.5d0) call random_number(random); deltaR(3,k) = delpos*(random-0.5d0) pos(:,k) = pos(:,k)+deltaR(:,k)/BoxSize call Compute_EAM_Energy(Ipbc,ene_new,NListMC) else call Compute_EAM_Energy(Ipbc,ene_old,0) call random_number(random); ix = idint(3.d0*random) + 1 call random_number(random) BoxSize(ix) = BoxSize(ix) * (1.0d0 + dellen*(random-0.5d0)) call Compute_EAM_Energy(Ipbc,ene_new,0) endif

24. Monte Carlo –Grand Canonical ensemble prefac = (Amass(itemp0)/Amass(itemp1))**(-1.5) if(itemp1.ne.itemp0) then if(itemp1.gt.itemp0) then; delmu1 = delmu else; delmu1 = (-1) * delmu; endif else delmu1 = 0 endif P01 = prefac*dexp((delmu1-ene_new+ene_old)/boltz/TRequested) call random_number(random) if(random.gt.P01) then IatomType(k) = itemp0 pos(:,k) = oldpos(:,1) BoxSize = OldBoxSize else Dk = dsqrt(dot_product(DisplaceList(:,k),DisplaceList(:,k))) if( Dk > Skin ) call Update_List(Range,delR,Lg_fun,0,Nequi,Nsave,Ipbc) delta_ene = ene_new - ene_old ene_sum = ene_sum + delta_ene endif enddo

25. Monte Carlo – Applications Segregation to an interface Order-disorder transition in binary alloy Morphology Evolution of an Alloy Particle MD+MC+MS Hybrid Simulation

26. Monte Carlo – Deposition or Etching

27. Monte Carlo – Determination of Phase Diagram A Modified Embedded Atom Method Interatomic Potential for the Cu-Ni SystemByeong-Joo Lee and Jae-Hyeok Shim, CALPHAD 28, 125-132 (2004).

28. Monte Carlo – Surface or Grain Boundary Segregation Grand Canonical Monte Carlo

29. MD: Tensile Impact on Fe, Fe-H, C, N

30. Limitation of MD Simulation • MD simulation time is • too short (a few ns) • to allow diffusion of atoms • D~10-12 m2/sec • t~10-8 sec • (Dt)1/2~10-10 m

31. MC+MD Hybrid Simulation • Effects of Interstitial Solute Atoms cannot be investigated using MD Simulations • because of not being able to allow probable diffusion of interstitial solute atoms • As a means to overcome the above limitation • MC+MD Hybrid Simulation can be considered

32. MC+MS Hybrid Simulation MC MS Possibility of atomic exchange becomes much natural

33. MC+MS: Diffusion of N in Fe-N

34. MC+MS+MD: Tensile Loading Fe-C Fe-N