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