Transcript CIT-Molecular-Mechanics-Lect-1

Fundamentals of Molecular
Dynamics for Nano-technology
Applications
Mario Blanco
Materials and Process Simulation Center
California Institute of Technology
IPAM, October 16, 2002
Outline
•
•
•
•
•
•
•
•
Hierarchical Multi-scale Modeling
Brief Review of QM and MD
Force Field Development
Long Range Atomic Potentials
Large Scale MD simulations
Fluid Control Nano-valve
NEMS
Challenges
Hierarchical Simulations
Chemist
Multi-scale
Years
Yards
Cancer
Research
Genetic
Engineering
Seconds
Inches
2   
Fossil Energy
Fuel Cells
Nanotechnology
C1
Chemistry
Organelle
Modeling
Receptor
Modeling
Pharmaceuticals
Polymers
Electronic
& Optical
Ceramics
Materials
Specialty
Chemicals
Metal
& Catalysts
Alloys
Materials
Catalysis
Microseconds
Microns

  
Biochemistry
Molecular
Self-Assembly
Picoseconds
Nanometers
Material
Science
Chemistry
Equilibrium
& Rate
Constants
Design
Meso-scale
Modeling
Molecules
F=ma
Molecular Dynamics
Force Fields
Femptoseconds
Angstroms
H = E 
QUANTUM
MECHANICS
Atoms
Electrons
© W.A. Goddard III, M. Blanco, 1998
Chemical Engineer
ELECTRONS => ATOMS => SEGMENTS => GRIDS
ENGINEERING
DESIGN
TIME
UNIT
PROCESS
DESIGN
years
hours
FINITE
ELEMENT
ANALYSIS
minutes
s
PROCESS
SIMULATION
seconds
MESOSCALE
DYNAMICS
= E 
microsec
Hierarchical
Modeling
Viscosity
Thermal Conductivity
Friction Coefficients
Activity Coefficients
Phase Diagram
nanosec
MOLECULAR
DYNAMICS
F= M A
Segment Averages
Group Additivities
Solubilities
QSPR
Equation of State
picosec
QUANTUM
MECHANICS
H=E
ASOG
UNIQUAC
Yield 
Young's 
Force Field
Charges
femtosec
1A
1 nm
10 nm
micron
mm
yards DISTANCE
Mechanical Engineer
time
ELECTRONS ATOMS
GRAINS
GRIDS
Continuum simulations of
real devices and materials
hours
Continuum
(FEM)
minutes
seconds
microsec
nanosec
MD
picosec
femtosec
Micromechanical modeling
MESO
Deformation and Failure
(dislocations, cracks, etc.)
QM
distance
Å
nm
micron
mm
cm
Transport properties
(diffusion, thermal transport, etc.)
meters
Accurate calculations for bulk phases and
molecules (EOS, dissociation curves, etc.)
New generation reactive force fields
based purely on first principles
Optical properties
For metals, oxides, organics.
Describes: mechanical properties, chemistry,
charge transfer, atomic polarization, etc.
MSC/Caltech
Quantum Mechanics (Schrodinger)
• Periodic Systems (3D & 2D)
• High Accuracy for Large Systems
• Fast Results for Large Systems
• Solvation (Poisson-Boltzmann)
Force Fields
• Polarizable, Charge Transfer
• Variable Bond Orders
• Phase Transitions
• Mixed Metal, Ceramic, Polymer
MesoScale Dynamics
• Coarse Grained FF
• Gas Diffusion
• Hybrid MD and Meso Dynamics
• Tribology
• Epitaxial Growth
Utilization: Web-based user-friendly
tools for nonexperts
Molecular Dynamics
• Large Systems (CMM, Parallel)
• Non-Equilibrium Dynamics
– Viscosity
– Thermal Conductivity
• Solvation (Schrodinger)
– surface tension, contact angles
• Hybrid QM/MD
• Hierarchical NEIMO
• Plasticity
– Formation Twins, Dislocations
– Crack Initiation
• Interfacial Energies
• Generalized Hildebrand Solubilities
Process Simulation
• Vapor-Liquid Equilibria
• Reaction Networks
Nanotechnology
Nano-electronic Devices
Nano electro-mechanical systems (NEMS)
The Hamiltonian
2
2 2
2
H  
 
 
k 2M k
i 2mi
k

l k
Z k Zl e2

Rkl
i

j i
Z k e2
e2

rij
rkj
kj
Born-Oppenheimer
Z
1
    2   k  
i 2
kj rkj
i

j i
  E
1
rij
E(R) = Energy Eigenvalues
Function of nuclear coordinates
n
   i
i 1
n electron trial
solution
M
 i   Ci 
1 e orbitals=LCAO

M sufficiently large leads to Hartree-Fock limit
Correlation Energy = Difference between experimental and HF energy
Slater Determinants
Pauli Exclusion Principle
Two electrons: anti-symmetrized linear combination
 
1
2
2 1s (1)1s (2) (1)  (2)   (1) (2)
1    
 2    
etc…
n electrons: Slater Determinant
electrons
orbitals
1 11 2 ...1 n 
 2 1 2 2 ... 2 n 
 3 1 3 2 ... 3 n 
  .
.
.
 n 1 n 2 ... n n 
1
n!
Potential Energy Surface (PES)
Bond
R
E ( R) 
Kb
R  Ro 2
2
Ro = 0.9 – 2.2 A
K = 700 Kcal/mol-A2
Potential Energy Surface (PES)
Non-Bond
R
Do=0.01 to 10 Kcal/mol
Ro=1.0 to 4.0 A
0.04
Van der Waals (induced dipoles)
0.03
E(R)= A e-CR - B/R6
E(R)=Do{[Ro/R)12]}- 2[Ro/R)6]}
0.02
Buckinham
Lennard-Jones 12-6
Exponential-6
0.01
0
EvdW ,ij
-0.01
 6 
 
R      Rij
 exp  1  ij    

 DvdW 
R


6
   6 

 RvdW
vdW







6



Electrostatics (point charges)
-0.02
-0.03
2
3
4
5
6
Ecoul ( Rij )  
j

i j
qi q j
Rij
q=-1 to 1
Force Fields
Molecular Mechanics
q = Ri nuclear coordinates (PES)
V(q) from Quantum Mechanics
Hessian Matrix: second derivatives
 2V
Vij 
Ri R j
Quantum Mechanics
K
ij
, Ro , o , qi 
Minimization
Yields equilibrium
structures
Molecular Dynamics
L  T  V (q)
 dR 
T  1 / 2 mi vi  1 / 2 mi  i 
 dt 
j
j
2
2
In cartesian coordinates
H   pi qi  L
j
H
  p i
qi
Generalized forces
H
 qi
pi
Generalized momentum
pi 
L
qi
q(t  t )  2q(t )  q(t  t ) 
p (t ) 2
t
m
Verlet algorithm
 (t ) 4
150 deg
C6F14 Dihedral QM energy (B3LYP, 631G*)
3.5
3
Energy (Kcal/mol)
2.5
164 deg
2
1.5
1
0.5
0
145
150
155
160
165
170
175
180
185
-0.5
Dihedral angle
180 deg
Energy Components
Molecular origin of helicity in Teflon
3
torsion
2.5
Energy (Kcal/mol)
vdw
2
electrostatic
Force Field
1.5
QM B3LYP
1
O
0.5
0
150
-0.5
-1
155
160
165
170
 (degrees)
175
180
Predictions
The simulated monoclinic (M phase) structure of C20F42,
looking through the crystal along the C direction on the left
and along the B-direction on the right.
(A = 9.65 Å; B = 5.70 Å; C = 28.3 Å,  = 97.2o;  =  = 90 o)
High Pressure forms of C20F42
The figure shows predicted stable helical conformations for C20F42. From left to right t+, t-, g+, g-, h+, and
h- enantiomeric pair conformations. The atoms are colored to facilitate the viewing of their helical nature.
The tighter the dihedral angle (from 164 to 60) the shorter the molecule gets. Fluorine atoms of each color
would be located on the same side if the molecule were prepared in the all-trans conformation.
MD simulation of Uniaxial Tension in Crystalline Cu nano-wire
(Tahir, Strachan, Goddard, 2000
•N=1370 atoms
•strain rate = 0.5% / 10 ps
= 5x109 1/s
•T = 300 K
•Failure at ~ 100%
•up to 25 %
6
5
4
Stres[GPa]
3
2
1
0
0 51
01
52
02
53
03
54
04
5
S
t
r
a
i
n
[
%
]
Failure: Impact Spallation
(Strachan & Goddard)
N=10,000 Ta atoms
Velocity= 6 km/s
MD allows the study of
spallation of metals vs. T
And impact speeds
Deformation: dislocation mobility from MD
simulation in Ta
Screw Dislocations Dipole annihilation (PBC)
MD simulations using First
Principles Force Fields:
•Core energy and structure of
dislocations
•Dislocation mobility (Peierls
stress, kink energies, etc.)
•Provides an accurate atomistic
description of the fundamental
unit mechanisms that control
plasticity
MD @ T = 0.001 K
5670 atoms
Transport properties:
gas diffusion in polymers
• High-frequency
short-wave-length
modes
of the polymer
are not critical for
diffusion
• The penetrant
molecule does not
affect
the polymer
dynamics
We developed the
Multiple Time Step
MD method to study
gas diffusion in
polymer.
• The model captures
the correct:
Thermal fluctuations
Time correlations of
the polymer
Features: New FF terms
Bonds
Angles
Dihedrals
Van der Waals
Coulomb
E  [k / 2 sin( 0 )][cos( )  cos( 0 )]2
1
E   k [1  d cos(n )]
2
E  [1/ 2k / sin(0 )][cos( )  cos(0 )]2
Cos harmonic angle
p
Dreiding Dihedral
Umbrella Inversion
Exponential-6
Exocyclic bonds
,n
n 1
6
E  D [(
)e 
6 
R
o
( 1 R / Ro )

R
(
)( ) ]
 6 R
o
6
Sandia Lammps
• parallelism via a spatial-decomposition algorithm
• long-range Coulombic interactions
– Ewald or PPPM (particle-mesh Ewald): Smeared of charges on 3d grid
followed by FFT transform of charge density, solution of Poisson
Equation, Differentiation (forces),remapping to atoms
• Force Fields:
– harmonic molecular potentials (bond, angle, torsion, improper)
– class II (cross-term) molecular potentials
• NVE, NVT, NPT dynamics
• constraints on atoms or groups of atoms
• rRESPA long-timescale integrator
• energy minimizer (Hessian-free truncated Newton method)
MSC-Lammps Added Features
• Extended Force Field Energy Functions
• Smooth Cubic Spline switch functions for non-bond
• Preprocessing of input files for automatic FF and EEXP
generation
• Trajectory File Generation
• PostProcessing of Thermodynamic Functions by VAC
Wilson’s Method
• Atomistic Energy Partitioning and Visualization
MSC-Lammps Programming Paradigm
Preprocessin
g
MD
Run
MD
P1Run
MD
P1Run
MD
P1Run
P1
Post
processing
Periodic Boundary Conditions
1x1x1
2x2x2
Ewald Trick
Ewald Equation
Charge Density
Ewald Forces
Basic PPPM Parallel Algorithm
•Interpolate smeared charges to grid
•Sum charge from grid ghost points
•Send charge values from spatial-to FFT-decomposition
•Solve Poisson’s equation on grid via FFTs.
•Send field values from FFT-to spatial-decomposition.
•Acquire fields for grid ghost points.
•Interpolate field values on grid to atoms.
Particle Mesh
Mesh Forces
Benchmarks (7,000 atoms)
Lammps 7000 Atom case
Parallelism (7,000 Atoms)
10
120
1
1
10
100
1000
% Efficiency
CPU Time (Secs)
100
80
60
40
0.1
20
0
1
0.01
1
7K
4.42
100.0
2
7K
2.20
100.5
4
7K
1.13
97.8
8
7K
0.602
91.9
100
Number of Processors
Number of Processors N
Procs
N
CPU
|| Eff
10
16
7K
0.324
85.3
32
64
128
256
512
7K
7K
7K
7K
7K
0.17 0.106 0.066 0.0465 0.0453
81.3 65.2 52.3
37.2
19.1
1000
Design of a Nanomechanical Fluid
Control Valve Based on the Deflection
of a Functionalized Silicon Cantilever
Coupling of Molecular Mechanics and
Mechanical Engineering Methods
Santiago Solares, Mario Blanco, and William A. Goddard III
Materials and Process Simulation Center, Caltech
July, 2002
micro ()
storage
Nano-Valve
Flow Control
nano-assembly
(nm)
SOME EARLY DESIGN LEARNINGS
“MOLDING” Vs. “CRIMPING”
USE OF SPLINES AND CUTOFFS
TWISTING OF Si(111) CANTILEVER
Si(100) CANTILEVER STRAIN
ENERGY ANALYSIS
Si100 CANTILEVER - Strain Energy Vs. Curvature
(corrected to have the plane of zero stress on the bottom face)
600.0
500.0
Energy, kcal/mol
y = 6E+07x2 + 15761x + 15.003
R2 = 0.9947
400.0
300.0
200.0
100.0
0.0
0
0.0005
0.001
0.0015
0.002
0.0025
0.003
Curvature, 1/Ang.
STRAIN ENERGY FOR A DEFLECTED FLAT SLAB = YwLoH3/24Rc2
AVERAGE CALCULATED YOUNG’S MODULUS: 76.7 Gpa (EXPERIMENTAL, BULK 47 GPa)
17,17 CARBON NANOTUBE
STRAIN ENERGY ANALYSIS
17,17 NANOTUBE - Strain Energy Vs. Curvature
(plane of zero stress bisects the tube)
1800.0
1600.0
1400.0
y = 345828x + 93.835
R2 = 0.9908
Energy, kcal/mol
1200.0
BUCKLING POINT
1000.0
SMOOTH
BENDING
BUCKLED
800.0
600.0
400.0
y = 1E+08x2 + 14067x - 1.6189
R2 = 1
200.0
0.0
0
0.0005
0.001
0.0015
0.002
0.0025
0.003
0.0035
0.004
0.0045
0.005
Curvature, 1/Ang
STRAIN ENERGY FOR A DEFLECTED HOLLOW ROD = pYLo(Ro4 – Ri4)/8Rc2
AVERAGE CALCULATED YOUNG’S MODULUS: 1719 Gpa, Experimental from vibrational
frequencies 1250 GPa, graphite 630 Gpa, CNT tensile molecular simulations 640 – 673 GPa
CARBON NANOTUBE “CRIMPING”
STRAIN ENERGY ANALYSIS
NEW DESIGN
76,500 ATOMS
Measuring Muscle Strength
with an AFM Tip
Atomic Force Microscopy (AFM) tip can be
attached to a surface-bound molecule
Cantilever deflection proportional to “pull” of molecular muscle
Dimerization vs. Complexation —
Not Necessarily the Same
Differentiate by NMR and mass spectrometry
Free Energy Analysis
• Density of Vibration State (Berens and Wilson,1981)
3N
1
1
S (v)  4p  m j
lim
2p   2
j 1



dt exp( i 2pvt)v j (t )
2


0
S (v)dv  3N
• Partition Function: Harmonic approximation

ln Q   dvS(v) ln q HO (v)
0
q HO (v) 
exp(  hv / 2)
1  exp(  hv)
classical limit
q HO (v) 
1
hv
Density Density
of States
Distribution
of States
Distribution
• Thermodynamics
5
1


0
Q
E
dvS(v)W (v)

S Q  k B  dvS(v)WSQ (v)
0
 hv

hv
WEQ (v)  


2
exp(

hv
)

1


4
Low freq vibs


hv
WSQ (v)  
 ln 1  exp(  hv) CH
3  exp( hv)  1

AQ  V0   1  dvS(v)WAQ (v)
 1  exp(  hv) 
WAQ (v)  ln

2
 exp(  hv / 2) 
V0  E MD  E C
E MD

0
1o
2c
3c
stretching
S(v) (cm)
E  V0  
Q
1
Total energy from MD


0
0
E C   1  dvS(v)WEC (v)   1  dvS(v)
0
0
1000
2000
v (cm-1)
3000
4000
Relative Stability Between Chain and Loop
AA
EE
-TS
500
100
61
23
0
7
-100
-133
+
+
+
+
+
+
-700
+
+
+
+
+
+
+
+
rotaxane type
1l
1o
+
+
+
+
+
+
-254
+
+
+
+
-248
+
-144
-300
-500
+
261
157
144
+
A (kcal/mol)
300
309
+
+
+
+
+
2c
2l
+
+
+
+
+
+
+
+
+
+
+
A(l-o) = 23 kcal/mol
+FSM = 29 kcal/mol
A(l-c) = -53 kcal/mol
+FSM = -50 kcal/mol
+
+
+
+
+
Free Energy of Concatenated Rotaxanes
Free Energy of Assembly
100
50
A
DA
E
DE
0 0
0 -TS
-DTS
34
7
61
21
67
40
0
-27
-48
+
+
+
rotaxane type
+
+
-250
+
+
+
+
+
+
+
+
-115
+
+
+
+
+
+
+
+
+
-200
+
+
+
+
+
+
+
-150
+
-100
+
-50
+
Energy (kcal/mol)
150
Opto-mechanical Muscles
•Azobenzene monomer response to light
•365 nm trans
cis
•420 nm cis
trans
extended
•Tension of a single strand ~ 205 pN
•Length change per monomer ~ 0.25 nm
•Not to exceed 50 nm in length for good
efficiency in Quantum yield
•Determination of mechanical properties
•Integration into NEMS
hn 365 nm)
contracted
hn 420 nm)
One designed example: A B C D groups key design issue
B
N
D
O
N
N
O
S
S
C
D
O
D
+
O
S
S
N
A
N
A
N
B
N
C
N
Challenges
• Multi-scale algorithms
– Particle-Continuum Algorithms
• Energy/Mass/Flux conservation
– More General embedding algorithms
• Atoms to “beads”
• Atoms to “strings”
• Atoms to complex topological objects
–
–
–
–
Space filling
Grain boundary properties
Embedding mechanical (elastic)
Viscoelastic, Thermal, and Chemical Properties
Challenges…
• Long-Range Interactions
– Beyond particle mesh methods ?
90
80
70
% CPU time (2 processors)
60
50
40
30
20
Task
Other
I/O
Fcomm
Comm
Exch
Nay-2
Nay-1
Dihed
Angle
Bond
0
Nbond=Short+Long
10
18,482 atoms
Nafion Fuel Cell
80 C
100 ps MD
148 hours 2 proc
Non-bond=127 hours
Acknowledgements
• Goddard Group
– Tahir Cagin, Alejandro Strachan
– Shiang-Tai Lin
– Santiago Solares
• Stoddard Group
• Andres Jaramillo