Atomistic simulations
This presentation is the property of its rightful owner.
Sponsored Links
1 / 47

Atomistic Simulations PowerPoint PPT Presentation


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

Atomistic Simulations. Byeong-Joo Lee Dept. of MSE Pohang University of Science and Technology (POSTECH) [email protected] Atomistic Simulation. I. General Aspects of Atomistic Modeling II. Fundamentals of Molecular Dynamics III. Fundamentals of Monte Carlo Simulation

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

Atomistic Simulations

Byeong-Joo Lee

Dept. of MSE

Pohang University

of Science and Technology

(POSTECH)

[email protected]


Atomistic simulation

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

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 simulation1

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

General Aspects of Atomistic Simulation– force & potential energy

j


General aspects of atomistic simulation equation of motion

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

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 simulation2

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

General Aspects of Atomistic Simulation– radial cutoff

Cutoffradius


General aspects of atomistic simulation potential truncation

General Aspects of Atomistic Simulation– potential truncation


General aspects of atomistic simulation neighbor list

General Aspects of Atomistic Simulation– Neighbor list

skin

Cutoffradius

Neighbor List is updated whenever maximum atomic displacement > ½ skin


Example potential force table

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

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 simulation3

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

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 simulation4

General Aspects of Atomistic Simulation

• Kinetic Energy

• Temperature

• Pressure

Virial equation

• Stress Tensor


General aspects of atomistic simulation stress tensor

General Aspects of Atomistic Simulation – Stress Tensor

Equilibrium system @ 0K

Equilibrium system @ a finite temperature

Thermally induced stress


Atomistic simulations

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


Atomistic simulations

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

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

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

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


Atomistic simulations

Molecular Dynamics – Lagrangian equation of motion

q: generalized coordinates

: time derivative of q

K: Kinetic energy

V: Potential energy


Atomistic simulations

Molecular Dynamics – MD at constant pressure


Atomistic simulations

Molecular Dynamics – MD at constant stress


Atomistic simulations

Molecular Dynamics – MD at constant stress: Parrinello and Rahman


Analysis to confirm equilibrium

Analysis – to confirm equilibrium


Analysis to confirm maintenance of structure

Analysis – to confirm maintenance of structure


Example computation of order parameter lambda

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

Analysis – Identification of structure


Atomistic simulations

Molecular Dynamics – Analyzing: radial distribution function

  • Partial RDFs of Cu50Zr50

  • during a rapid cooling

    • Cu–Cu pair

    • Cu–Zr pair

    • (c) Zr–Zr pair


Atomistic simulations

Molecular Dynamics – Analyzing: radial distribution function


Atomistic simulations

Molecular Dynamics – Analyzing: angle distribution function


Example sampling data for g r and bond angle distribution

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

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


Atomistic simulations

Molecular Dynamics – Analyzing: localatomistic structure


Atomistic simulations

Molecular Dynamics – Analyzing: localatomistic structure


Analysis identification of structure1

Analysis – Identification of structure


Atomistic simulations

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


Atomistic simulations

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)


Atomistic simulations

Computing Physical Properties– Elastic Moduli

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


Atomistic simulations

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.


Atomistic simulations

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


Atomistic simulations

Computing Physical Properties– Elastic Moduli

for Cubic for HCP


Atomistic simulations

Computing Physical Properties– Ni3Al/Ni Interfacial Energy

σ =

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

/ 2A


Atomistic simulations

Computing Physical Properties– Grain Boundary Energy


Atomistic simulations

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


  • Login