1 / 47

Atomistic Simulations

Atomistic Simulations. Byeong-Joo Lee Dept. of MSE Pohang University of Science and Technology (POSTECH) calphad@postech.ac.kr. Atomistic Simulation. I. General Aspects of Atomistic Modeling II. Fundamentals of Molecular Dynamics III. Fundamentals of Monte Carlo Simulation

zudora
Download Presentation

Atomistic Simulations

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

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

  2. Atomistic Simulation I. General Aspects of Atomistic Modeling II. Fundamentals of Molecular Dynamics III. Fundamentals of Monte Carlo Simulation IV. Semi-Empirical Atomic Potentials (EAM, MEAM) V. Analysis, Computation of Physical Quantities VI. Exercise using Simulation Code

  3. General Aspects of Atomistic Simulation • Objectives • · Equilibrium configuration of a system of interacting particles, • Analysis of the structure and its relation to physical properties. • · Dynamic development • Applications • · Crystal Structure, Point Defects, Surfaces, Interfaces, • Grain Boundaries, Dislocations, Liquid and glasses, … • Challenges • · Time and Size Limitations • · Interatomic Potential

  4. General Aspects of Atomistic Simulation • Equation of Motion • Potential Energy(e.g. pairwise interaction) • Potential Truncation · Radial Cutoff · Long-range Correction · Neighbor List • Periodic Boundary Condition · Minimum image criterion · Free surfaces

  5. General Aspects of Atomistic Simulation– force & potential energy j

  6. General Aspects of Atomistic Simulation –equation of motion @ time = to, assume that the force will not change during the next ∆t period @ time = to+ ∆t, based on the newly obtained atom positions, forces on all individual atoms are newly calculated. and, assume again that the forces will not change during the next ∆t period

  7. General Aspects of Atomistic Simulation– size of Δt How large should be the ∆t ? More than several tens of time steps should be spent during a vibration ∆t≈10-15 sec. for most solid elements

  8. General Aspects of Atomistic Simulation • Equation of Motion • Potential Energy(e.g. pairwise interaction) • Potential Truncation · Radial Cutoff · Long-range Correction · Neighbor List • Periodic Boundary Condition · Minimum image criterion · Free surfaces

  9. General Aspects of Atomistic Simulation– radial cutoff Cutoffradius

  10. General Aspects of Atomistic Simulation– potential truncation

  11. General Aspects of Atomistic Simulation– Neighbor list skin Cutoffradius Neighbor List is updated whenever maximum atomic displacement > ½ skin

  12. Example – Potential & Force Table do i=1,Ncomp eps(i,i) = epsilon(i) sig(i,i) = sigma(i) sig12(i,i) = sig(i,i)**12 sig6(i,i) = sig(i,i)**6 enddo c cutoff potential Phicutoff = 4.d0*eps*(sig12/(Rcutoff**12) - sig6/(Rcutoff**6)) c c Make a table for pair potential and its derivative between each pair RsqMin = Rmin**2 DeltaRsq = ( Rcutoff**2 - RsqMin ) / ( LineTable - 1 ) OverDeltaRsq = 1.d0 / DeltaRsq c do k=1,LineTable Rsq = RsqMin + (k-1) * DeltaRsq rm2 = 1.d0/Rsq rm6 = rm2*rm2*rm2 rm12 = rm6*rm6 c c 4*eps* [s^12/r^12 - s^6/r^6] - phi(Rc) PhiTab(k,:,:) = 4.d0*eps * (sig12*rm12 - sig6*rm6) - phicutoff c c The following is dphi = -(1/r)(dV/dr) = 24*eps*[2s^12/r^12-s^6/r^6]/r^2 DPhiTab(k,:,:) = 24.d0*eps*rm2 * ( 2.d0*sig12*rm12 - sig6*rm6 ) enddo

  13. Example – Neighbor List RangeSq = Range*Range L = 1 c do i = 1,Natom Mark1(i) = L istart = i+1 do j = istart,Natom Sij = pos(:,i) - pos(:,j) where ( abs(Sij) > 0.5d0 .and. Ipbc.eq.1 ) Sij = Sij - sign(1.d0,Sij) end where c go to real space units Rij = BoxSize*Sij Rsqij = dot_product(Rij,Rij) if ( Rsqij < RangeSq ) then List(L) = j L = L + 1 endif enddo Mark2(i) = L - 1 enddo

  14. General Aspects of Atomistic Simulation • Equation of Motion • Potential Energy(e.g. pairwise interaction) • Potential Truncation · Radial Cutoff · Long-range Correction · Neighbor List • Periodic Boundary Condition · Minimum image criterion · Free surfaces

  15. General Aspects of Atomistic Simulation– Periodic Boundary Condition do i=1,Natom read(2,*) PosAtomReal(i) pos(:,i) = PosAtomReal/BoxSize enddo ------------------------------------- Sij = pos(:,i) - pos(:,j) where ( abs(Sij) > 0.5d0 .and. Ipbc.eq.1 ) Sij = Sij - sign(1.d0,Sij) end where When PBC is applied, Rcutoff should be < than the half of sample dimension

  16. General Aspects of Atomistic Simulation • Kinetic Energy • Temperature • Pressure Virial equation • Stress Tensor

  17. General Aspects of Atomistic Simulation – Stress Tensor Equilibrium system @ 0K Equilibrium system @ a finite temperature Thermally induced stress

  18. Molecular Dynamics – Time Integration of equations of motion Integration Algorithm Computing time and memory Numerical error Stability Energy/momentum conservation Reversibility The Verlet algorithm Leapfrog algorithm

  19. Molecular Dynamics – Time Integration of equations of motion Velocity Verlet algorithm Integration Algorithm Computing time and memory Numerical error Stability Energy/momentum conservation Reversibility Gear’s Predictor-corrector algorithm

  20. Example – Time integration of equations of motion c dr = r(t+dt) - r(t) and updating of Displacement deltaR = deltat*vel + 0.5d0*deltatsq*acc pos = pos + deltaR c BoxSize rescale for constant P if(ConstantP) then BoxSize = BoxSize + deltat*VelBoxSize +.5d0*deltatsq*AccBoxSize Volume = product(BoxSize) Density = dble(Natom) / Volume deltaV = Volume - Vol_Old Recipro(1) = BoxSize(2) * BoxSize(3) Recipro(2) = BoxSize(3) * BoxSize(1) Recipro(3) = BoxSize(1) * BoxSize(2) VelBoxSize = VelBoxSize + 0.5d0*deltat*AccBoxSize endif c velocity rescale for constant T -> v(t+dt/2) if(ConstantT .and. (temperature > 0) ) then chi = sqrt( Trequested / temperature ) vel = chi*vel + 0.5d0*deltat*acc else vel = vel + 0.5d0*deltat*acc endif c a(t+dt),ene_pot,Virial call Compute_EAM_Forces(Ipbc,f_stress) c add Volume change term to Acc when ConstantP if(ConstantP) then do i=1,Natom acc(:,i) = acc(:,i) - 2.0d0 * vel(:,i) * VelBoxSize/BoxSize enddo endif c v(t+dt) vel = vel + 0.5d0*deltat*acc

  21. Molecular Dynamics – Running, measuring and analyzing • 1. Build-up : initial configuration and velocities • 2. Speed-up : Radial cut-off and Neighbor list • 3. Controlling the system • · MD at constant Temperature • · MD at constant Pressure • · MD at constant Stress • 4. Measuring and analyzing • · Average values of physical quantities • · Physical Interpretation from statistical mechanics

  22. Molecular Dynamics – MD at constant temperature Scaling of velocities : To : Target value of temperature T : instantaneous temperature v(t+Δt/2) = sqrt (To/T)v(t) + ½ a(t)Δt

  23. Molecular Dynamics – Lagrangian equation of motion q: generalized coordinates : time derivative of q K: Kinetic energy V: Potential energy

  24. Molecular Dynamics – MD at constant pressure

  25. Molecular Dynamics – MD at constant stress

  26. Molecular Dynamics – MD at constant stress: Parrinello and Rahman

  27. Analysis – to confirm equilibrium

  28. Analysis – to confirm maintenance of structure

  29. Example – Computation of order parameter, Lambda if(Lh_fun.and.mod(mdstep,Nsampl).eq.0 .and. mdstep.le.Nequi) then Alambda = dcos(FourPi*pos(1,1)*xcell) & + dcos(FourPi*pos(2,1)*ycell) & + dcos(FourPi*pos(3,1)*zcell) Alam1 = 0.d0 do i=2,Natom Alambda = Alambda + dcos(FourPi*pos(1,i)*xcell) & + dcos(FourPi*pos(2,i)*ycell) & + dcos(FourPi*pos(3,i)*zcell) Alam1 = Alam1 + dcos(FourPi*(pos(1,i)-pos(1,1))*xcell) & + dcos(FourPi*(pos(2,i)-pos(2,1))*ycell) & + dcos(FourPi*(pos(3,i)-pos(3,1))*zcell) enddo Alambda = Alambda / dble(3*Natom) Alam1 = Alam1 / dble(3*Natom-3) endif

  30. Analysis – Identification of structure

  31. Molecular Dynamics – Analyzing: radial distribution function • Partial RDFs of Cu50Zr50 • during a rapid cooling • Cu–Cu pair • Cu–Zr pair • (c) Zr–Zr pair

  32. Molecular Dynamics – Analyzing: radial distribution function

  33. Molecular Dynamics – Analyzing: angle distribution function

  34. Example – Sampling data for g( r) and bond-angle distribution do i = 1,Natom do j = istart,Natom Sij = pos(:,i) - pos(:,j) call Rij_Real(Rij,Sij) Rsqij = dot_product(Rij,Rij) NShell = idint(sqrt(Rsqij)/delR + 0.5d0) Ngr(Nshell) = Ngr(Nshell) + 1 c if ( Rsqij < RmSq ) then do k = j+1,Natom Sik = pos(:,i) - pos(:,k) call Rij_Real(Rik,Sik) Rsqik = dot_product(Rik,Rik) if ( Rsqik < RmSq ) then Num_Ang = Num_Ang + 1 costheta = dot_product(Rik,Rij) /dsqrt(Rsqij) /dsqrt(Rsqik) degr = 180.d0 * dacos(costheta) / pi NShell = idint(degr + 0.5d0) Nag(Nshell) = Nag(Nshell) + 1 endif enddo endif enddo enddo

  35. Example – Printing data for g( r) and bond-angle distribution c print g(r) function values every Ng_fun's time step with normalization c Ntotal = Num_Update * Natom totalN = dble(Num_Update * Natom) / 2.d0 c do i=1,160 gr = 0.d0 if(Ngr(i) .gt. 0) then radius = delR * dble(i) vshell = FourPiDelR * radius * radius gr = dble(Ngr(i)) / (totalN * Density * vshell) write(9,'(1x,f8.3,i10,f13.6)') radius,Ngr(i),gr endif enddo c c print bond-angle distribution c do i=1,180 if(Nag(i) .gt. 0) then gr = dble(Nag(i)) / dble(Num_Ang) write(9,'(1x,i8,f16.4)') i,gr endif enddo

  36. Molecular Dynamics – Analyzing: localatomistic structure

  37. Molecular Dynamics – Analyzing: localatomistic structure

  38. Analysis – Identification of structure

  39. Molecular Dynamics – Computing Physical Properties • 1. Bulk Modulus & Elastic Constants • 2. Surface Energy, Grain Boundary Energy & Interfacial Energy • 3. Point Defects Properties • · Vacancy Formation Energy • · Vacancy Migration Energy • · Interstitial Formation Energy • · Binding Energy between Defects • 4. Structural Properties • 5. Thermal Properties • · Thermal Expansion Coefficients • · Specific Heat • · Melting Point • · Enthalpy of melting, …

  40. Computing Physical Properties– Elasticity Tensor Most general linear relationship between stress tensor and strain tensor components 81 constants 36 constants 21 constants for the most general case of anisotropy Further reduction of the number of independent constants can be made by rotating the axis system, but the final number depends on the crystal symmetry . Short-hand matrix notation: tensor 11 22 33 23 31 12 matrix 1 2 3 4 5 6 examples C1111 = C11, C1122 = C12, C2332 = C44 for Isotropic Crystals for Cubic Crystals for Hexagonal Crystals C11(= C22), C12, C13(= C23), C33, C44(= C55), C66, (= (C11- C12)/2)

  41. Computing Physical Properties– Elastic Moduli Change of the potential energy and total volume upon application of the strain

  42. Computing Physical Properties– Elastic Moduli In practical calculations the evaluation of the elastic moduli is done by applying suitable small strains to the block of atoms, relaxing, and evaluating the total energy as a function of the applied strain. In this calculation the internal relaxations are automatically included. The elastic moduli are then obtained by approximating the numerically calculated energy vs. strain dependence by a second or higher order polynomial and taking the appropriate second derivatives with respect to the strain.

  43. Computing Physical Properties– Bulk Moduli epsilon = 1.0d-03 fac = 1.6023d+00 / epsilon / epsilon Volume = product(BoxSize) AV = Volume / dble(NAtom) c c B : increase BoxSize() by a factor of +- epsilon c BoxSize = Old_BoxSize * (1.d0 + epsilon) call Compute_EAM_Energy(Ipbc,ene_sum,0) ene_ave_p1 = ene_sum / dble(Natom) c BoxSize = Old_BoxSize * (1.d0 - epsilon) call Compute_EAM_Energy(Ipbc,ene_sum,0) ene_ave_m1 = ene_sum / dble(Natom) c B0 = fac * (ene_ave_p1 + ene_ave_m1 - 2.d0*ene_ave_0) / AV / 9.d0 write(*,*) ' B = ', B0, ' (10^12 dyne/cm^2 or 100Gpa)'

  44. Computing Physical Properties– Elastic Moduli for Cubic for HCP

  45. Computing Physical Properties– Ni3Al/Ni Interfacial Energy σ = {E( ) – [E( )+E( )]} / 2A

  46. Computing Physical Properties– Grain Boundary Energy

  47. Computing Physical Properties– Melting Point (Interface method) • Prepare two samples, one is crystalline and the other is liquid • Estimate the melting point • Perform an NPT MD run for solid at the estimated MP, and determine sample dimensions (apply 3D PBC) • Perform an NPT MD run for liquid, keeping the same sample dimensions in y and z directions with the solid • Put together the two samples in the x direction • Perform NVT MD runs at various sample size in the x direction, then perform NVE MD runs • ※ Partial melting or solidification and change of temperature will occur depending on the estimated mp. in comparison with real mp. • The equilibrium temperature between solid and liquid phases at zero external pressure is the melting point of the material • ※If full melting or solidification occurs the step 2-6 should be repeated with newly estimated melting point

More Related