OUTLINE
- Intro to Simulation
- Monte Carlo
- Molecular Dynamics
- MD Implementations
- Equations
- Algorithms
- Restrictions
- Applications
OVERVIEW
- Simulations serve as a bridge between microscopic and macroscopic theory and experiment.
- Predictions are restricted by computer budget and approximations
- Possible to carry out simulations that are difficult or impossible in experimental settings
- high temperatures or pressures
- Molecular Dynamics (MD) and Monte Carlo are the two main simulation techniques
- Hybrids of the two techniques are also widely employed
- Main focus of this presentation will be on MD simulations
MONTE CARLO APPROACH
- Class of computation algorithms that rely on repeated random sampling to compute their results.
- Commonly used when simulating physical and mathematical systems
- Name is derived from the Monte Carlo Casino
- Method involves guessing random solutions and iterating from the results
- Similar to a game of battleship
- Covers a wide variety of approaches
- The Monte Carlo method is just one of many methods for analyzing propogation
- Relies on statistical mechanics rather than molecular dynamics
- Generates states according to Boltzmann probabilities.
- Stochastic in Nature
- Determines new state from previous one through Markov Processes
- Dynamics are time invariant
- Useful for situations dealing with Boltzmann distributions.
MOLECULAR DYNAMICS APPROACH
- Uses numerical methods and approximations of known physics to determine the interactions of an array of molecules.
- Provides a view of the path of the molecules/atoms over time
- Gives a route to dynamical properties of the system through equations from Newton, Hamilton, Lagrange et al.
- Transport coefficients
- Time dependent responses to peturbations
- Rheological properties
- Spectra
- Allows for study of time history dependent properties
FORCE CALCULATIONS
- Atomic forces are determined through the potential energy between the two atoms.
- The force terms are derivatives of energy expressions (Lagrangian function) in which the energy of atom is typically written as a function of its position with respect other atoms.
f = atomic force
r = radii between atoms
U = potential energy function
- Potential energy between each molecule must be calculated.
- This potential may be approximated through many different equations
- When choosing potentials the following characteristics should be considered:
- Accuracy:
- To reproduce properties of interest as closely as possible.
- Transferability:
- Can be used to study a variety of properties for which it was not fit.
- Computational speed:
- Calculations are fast with simple potentials.
- Accuracy:
- Part of the potential energy representing non-bonded interactions
- Traditionally split into 1-body, 2-body, 3-body etc
- u(r) represents an externally applied potential field – usually dropped for fully periodic simulations of bulk systems
- v(r) represents pair potentials
- Three body interactions are typically neglected
LENNARD JONES POTENTIAL
- Most commonly used form for determining pair potential
- r12 describes repulsion, r6 describes attraction
- Attractive potential derived from dispersion interactions
- To save computation time, may be truncated at r = 2.5s
WCA POTENTIAL
- Weeks-Chandler-Anderson
- Applied when attractive interactions are of less concern than the excluded volume effects which dictate molecular packing
Plot comparing Lennard Jones Potential and WCA potentials
FORCE DERIVATION
Original Equation:
Take the derivative:
And solve for the force:
INTRAMOLECULAR FORCES
- For molecules, intramolecular bonding interactions must also be considered
- Bonds involve separation between adjacent pairs of atoms in a molecular framework
Geometry of a simple chain molecule (Allen)
Where:
SOLN ALGORITHMS
- Uses Newton’s equations of motion
- Consider a system of atoms
- Solve the system of equations for the motion of each atom due to intermolecular forces.
- Perform the computationally intensive force calculations as infrequently as possible
VERLET ALGORITHM
- Method used to integrate Newton’s equations of motion
- Frequently used to calculate trajectories of particles in MD simulations
- Offers greater stability over Euler method, and time reversibility
- It is obtained from the combination of two Taylor expansions in different time directions (t+Dt), (t-Dt)
- It has a local truncation error that varies as (t4) and hence is third order.
- The acceleration is obtained from the intermolecular potential and Newton’s second law.
- Various implementations of the Verlet algorithm
- Leapfrog
- Original
- Velocity
- Examine the velocity Verlet algorithm
- Advances the coordinates and momenta over a timestep dt
- The velocity Verlet algorithm is:
- Exactly time reversible
- Volume in phase space is conserved
- Low order in time, permitting long timesteps (high derivatives not stored)
- Requires just one expensive calc per cycle
- Easy to program
PSEUDOCODE
GEAR PREDICTOR-CORRECTOR ALGORITHMS
- Gear Predictor – Corrector algorithms are composed of three steps:
- Prediction of atom positions and their time derivatives.
- Force evaluation
- Correction of atom positions and their time derivatives using the discrepancy between the predicted acceleration and that given by the evaluated force (2nd. Newton’s Law).
Predictor Step:
Force Calculation Step:
Corrector Step:
CONSTRAINTS
Treat bonds as constrained to have fixed length
Define constant bond length b, write constraint equations
In Lagrangian formulation, constraint forces will be represented thus
Where L is the undetermined multiplier and
- Boundary conditions are satisfied through forcing constraint satisfaction at the end of each time step through SHAKE and RATTLE.
- Shake calculates forces lgi to ensure that the end of step positions satisfies the first constraint equation
- Rattle calculates mgi to ensure that the end of step momenta satisfies the second constraint equation
- Introduce forces to maintain constant bond length and match constraints after each time step
PSEUDOCODE
MULTIPLE TIMESTEPS
- MD method is able to tackle systems with multiple time scales
- Adopt a timestep short enough to handle the fastest variables
- Handle fast motions with a shorter timestep
- Fast varying forces computed many times at short intervals
- Slow varying forces are used just before and just after this stage, calculated only once per timestep
- Assume two types of forces present
- “Fast” forces fi
- “Slow” forces Fi
PSEUDOCODE
LENGTH OF SIMULATION
- MD evolves a finite-sized molecular configuration forward in time, in a step-by-step fashion.
- Limits on typical time scales and length scales that can be investigated
- Simulation runs are usually short, typically 103 – 106 MD steps corresponding to a few nanoseconds of real time.
- Simulation run length t must be large compared with ta, or the correlation time.
- Correlation time may be described as the time to decay from initial value to long-time limiting value.
BOUNDARY CONDITIONS
- Free Boundary Conditions
- Atoms are allowed to move freely at the outer edges of the cell – works for a molecule or cluster in vacuum
- A free boundary condition can be also applied for fast processes when the effects of boundaries are not important due to the short time – scale of the process.
- Rigid Boundary Conditions
- Atoms at the boundary are fixed – unable to move
- Used in conjuction with other boundary conditions
- Periodic Boundary Conditions
- In cases of small sample size, significant numbers of atoms may be along the outer edges (limited interaction with other atoms)
- Solved through introducing a periodic boundary conditions – surround the volume with replicas of itself
- Each computational cell must be larger than 2 rc (minimum image criterion).
- Neighbor Listing
- Each atom stores a list of neighboring molecules which fall inside of the radius of interest
- Allows the program to quickly run through the list and perform force calculations
- Break the domain into regions
- Only search for nearby atoms (within the radii r) in the neighboring domains
- Maintain minimum image criterion (computational cells must be larger than 2 rc)
THIN FILM SILICON THERMAL CONDUCTIVITY
Determine the thermal conductivity of thin film silicon used with silicon on insulator (SOI) transistors through Molecular Dynamics
PROBLEM STATEMENT
- Trend towards miniaturization of electronic devices has led to device features in the nanometer range
- Silcon transistors are expected to reach a gate length of 28 nm by 2009.
- When the dimensions of the material are comparable to the phonon free mean path, boundary or interface scattering can limit the mean free path and affect the thermal conductivity
- Apply equilibrium molecular dynamics to the study of the thermal conductivities of silicon thin films
- Allows the prediction of both in-plane and out-of-plane thermal conductivities for comparison to existing in-plane experimental data
EQUATIONS
Stillinger-Weber potential used - accurate for determining thermal conductivity
Total Potential:
2 - Body term:
3 - Body term:
Additional Term for atoms near fixed surface or plane
- Similar to the two-body pair potential from earlier, but is dependent on distance between an atom and the plane parallel to the nearest surface of the film instead of the interatomic distance.
- A, B, p, q, l, g are constants derived by Stillinger and Weber (1985)
- e = potential well depth
- rc = cutoff radius
- rij = relative distance between atoms i, j
MPI
- Used PSC (Pittsburgh super-computing center) to reduce processing time.
- Allowed for use of multiple processors
- Split the space up into user-defined domains
- MPI stands for Message Passing Interface
- Protocol used for communication between multiple processors
- Consists of set routines in a variety of programming languages
- FORTRAN
- C
- C++
- Any other language capable of interfacing with routine libraries
- Used to divide large domain into smaller chunks to reduce processing time
- Each processor handles the force calculations for the communication cell assigned to it.
- As seen in previous figures, atoms in one communication cell may exert influence on atoms in a different communication cell.
- Therefore, the processors must pause and pass messages to each other before each force calculation.
- RESULTS
- Thermal conductivity is determined from the Green-Kubo relation between thermal conductivity and heat current autocorrelation function
- Determined from the trajectory of the particles
- kMD = thermal conductivity
- tMD = integration time
- n = dimensionality (3 for bulk, 2 for in plane directions, 1 for out of plane directions)
- V = ensemble volume
- Dt = simulation timestep
- kB = Boltzmann’s constant
- TMD = temperature of MD simulation
- M = number of time steps required for heat current autocorrelation function to decay to zero
- Ns = number of heat current autocorrelation function averages
- J = heat current
- Successfully plotted the relation between silicon thickness and thermal conductivity
- Results were found to be consistent with kinetic theory
- k = thermal conductivity
- C = specific heat
- v = phonon group velocity
- L = phonon mean free path
- Also able to determine:
- Average kinetic energy
- Internal energy
- Pressure
- Heat Capacity at constant volume
- Diffusion coefficient
- Viscosity
CONCLUSIONS
- Using equilibrium molecular dynamics, studied thermal conductivity of silicon thin films in a broad range of thicknesses and temperatures
- Found that thermal conductivity is reduced with film thickness, is proportional to the effective mean free path.
- Results are in agreement with kinetic theory and the equation for phonon radiative transfer.
- Phonons mean free path is limited by the scattering with the boundaries of the film
Pittsburgh Supercomputing Center
- Joint effort between Pitt, CMU, and Westinghouse
- Provides researchers access to several of the most powerful computing systems available for unclassified research.
- Class in computation nanomechanics being offered next semester
- Will cover essentials on MD simulations to study material behavior
- Allow for hands on experience
Texts:
Michael P. Allen. “Introduction to Molecular Dynamics Simulation” Computational Soft Matter: From Synthetic Polymers to Proteins, Lecture Notes, Norbert Attig, Kurt Binder, Helmut Grubm¨ uller, Kurt Kremer (Eds.), John von Neumann Institute for Computing, J ¨ ulich, NIC Series, Vol. 23, ISBN 3-00-012641-4, pp. 1-28, 2004. ( Read at: http://www.fz-juelich.de/nic-series/volume23/allen.pdf
)
Michael P. Allen, Glenn T. Evans, Daan Frenkel, and Bela Mulder. "Hard convex body fluids." Adv. Chem. Phys., 86:1.166, 1993.
Jayeeta Ghosh and Roland Fallera. “Comparing the density of states of binary Lennard-Jones glasses in bulk and film.” The Journal of Chemical Physics, 128
Carlos J. Gomes, Marcela Madrid, Javier Goicochea, Christina H. Amon. “In-Plane and Out-Of-Plane Termal Conductivity of Silicon Thin Films Predicted by Molecular Dynamics.” ASME Journal of Heat Transfer, Vol 128, pp 1114-1121, 2006.
G. C. Maitland, M. Rigby, E. B. Smith, and W. A. Wakeham. Intermolecular forces: their origin and determination. Clarendon Press, Oxford, 1981.
D. C. Rapaport. The Art of Molecular Dynamics Simulation. Cambridge University Press, 2004
http://books.google.com/books?id=wusLQijMy8AC
J Seyler and M N Wybourne. "The effect of surface roughness on the phonon mean free path in narrow wires." Journal of Physics: Condensed Matter, Vol 2 p 8853
A. J. Stone. The Theory of Intermolecular Forces. Clarendon Press, Oxford, 1996.
L. Verlet. “Computer Experiments on Classical Fluids - Thermodynamical Properties of Lennard-Jones molecules.” Phys. Rev., 159:98.103, 1967.
L. Verlet. “Computer Experiments on Classical Fluids- Equilibirum Correlation Functions.” Phys. Rev., 165:1.201, 1968.
Jan Zielkiewicz. “Entropy of water calculated from harmonic approximation: Estimation of the accuracy of method,” J. Chem. Phys. 128, 196101 (2008), DOI:10.1063/1.2921161
Sites:
Overview/Introduction
- http://www.fz-juelich.de/nic-series/volume23/allen.pdf
- http://books.google.com/books?id=wusLQijMy8AC
Wikipedia:
- http://en.wikipedia.org/wiki/Monte_Carlo_method
- http://en.wikipedia.org/wiki/Molecular_dynamics
- http://en.wikipedia.org/wiki/Lennard-Jones_potential
- http://en.wikipedia.org/wiki/Verlet_integration
Video Demonstration
Useful Journal Search
Local Supercomputer Center
Further information on FENE potential
More Monte Carlo Information
Stillinger-Weber Derivations
Gear Predictor-Corrector Algorithm
