- 118 Views
- Uploaded on
- Presentation posted in: General

Atomistic Simulations

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.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Byeong-Joo Lee

Dept. of MSE

Pohang University

of Science and Technology

(POSTECH)

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

- 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

• 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

j

@ 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

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

• 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

Cutoffradius

skin

Cutoffradius

Neighbor List is updated whenever maximum atomic displacement > ½ skin

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

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

• 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

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

• Kinetic Energy

• Temperature

• Pressure

Virial equation

• Stress Tensor

Equilibrium system @ 0K

Equilibrium system @ a finite temperature

Thermally induced stress

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

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

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

- 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

Scaling of velocities :

To : Target value of temperature

T : instantaneous temperature

v(t+Δt/2) = sqrt (To/T)v(t) + ½ a(t)Δt

Molecular Dynamics – Lagrangian equation of motion

q: generalized coordinates

: time derivative of q

K: Kinetic energy

V: Potential energy

Molecular Dynamics – MD at constant pressure

Molecular Dynamics – MD at constant stress

Molecular Dynamics – MD at constant stress: Parrinello and Rahman

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

Molecular Dynamics – Analyzing: radial distribution function

- Partial RDFs of Cu50Zr50
- during a rapid cooling
- Cu–Cu pair
- Cu–Zr pair
- (c) Zr–Zr pair

Molecular Dynamics – Analyzing: radial distribution function

Molecular Dynamics – Analyzing: angle distribution function

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

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

Molecular Dynamics – Analyzing: localatomistic structure

Molecular Dynamics – Analyzing: localatomistic structure

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, …

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)

Computing Physical Properties– Elastic Moduli

Change of the potential energy and total volume upon application of the strain

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.

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)'

Computing Physical Properties– Elastic Moduli

for Cubic for HCP

Computing Physical Properties– Ni3Al/Ni Interfacial Energy

σ =

{E( ) – [E( )+E( )]}

/ 2A

Computing Physical Properties– Grain Boundary Energy

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