- 147 Views
- Uploaded on

Download Presentation
## PowerPoint Slideshow about ' Atomistic Simulations' - zudora

**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.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 - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript

Atomistic Simulations

Byeong-Joo Lee

Dept. of MSE

Pohang University

of Science and Technology

(POSTECH)

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

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

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

General Aspects of Atomistic Simulation– force & potential energy

j

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

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

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

General Aspects of Atomistic Simulation– radial cutoff

Cutoffradius

General Aspects of Atomistic Simulation– potential truncation

General Aspects of Atomistic Simulation– Neighbor list

skin

Cutoffradius

Neighbor List is updated whenever maximum atomic displacement > ½ skin

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

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

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

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

General Aspects of Atomistic Simulation

• Kinetic Energy

• Temperature

• Pressure

Virial equation

• Stress Tensor

General Aspects of Atomistic Simulation – 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

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

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

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

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

Analysis – to confirm equilibrium

Analysis – to confirm maintenance of structure

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

Analysis – Identification of structure

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

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

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

Molecular Dynamics – Analyzing: localatomistic structure

Molecular Dynamics – Analyzing: localatomistic structure

Analysis – Identification of 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– 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

Download Presentation

Connecting to Server..