NWChem: Molecular Dynamics and QM/MM
Molecular Dynamics Functionality
Target systems:
biomolecules (proteins, DNA/RNA, biomembranes)
General features:
energy evaluation (SP) energy minimization (EM) molecular dynamics simulation (MD) free energy evaluation (MSTP & MCTI) quantum molecular dynamics (QMD) hybrid molecular dynamics (QM/MM)
Force Fields Classical empirical force fields Bond-stretching, angle-bending, torsional, out-of-plane-bending, electrostatic, and van der Waals non-bonded interactions
Electronic polarization First order Self-consistent induced fields
Smooth particle-mesh Ewald electrostatics Effective pair potential MD only
Self-guided MD
Simulation
Integration Newton’s Equations of Motion (leapfrog, Browne-Clark)
Constant Temperature and Pressure (Berendsen weak coupling)
Periodic Boundary Conditions (minimum image convention)
Geometry Optimization (steepest descent, conjugate gradient)
Twin-Range Verlet Neighbor Lists (cell index method)
Constraints (SHAKE)
Free Energy Evaluation Single Step Thermodynamic Perturbation (SSTP)
∆G = G1 - G0 = - RT ln < exp ( - ∆ H/RT)>0 with ∆ H=H1-H0
Multiple Step Thermodynamic Perturbation (MSTP)
DG = Si ( Gi+1 - Gi ) = -Si RT ln < exp ( - ∆ Hi /RT)>i with ∆ Hi=Hi+1 – Hi
Multiconfiguration Thermodynamic Integration (MCTI) DG = Si ( Gi+1 - Gi ) = Si < ∂ H(λ) / ∂λ >i Dλ i
Single and dual topology Hamiltonians Double-wide sampling Separation-shifted scaling (SSS) Potentials of mean force over multiple processors Statistical error correlation analysis
Database Directories ffield_lvl
ffield
force field e.g. amber, charmm lvl 1, 2, 3, 4, 5, 6, 7, 8, or 9
Directories supplied with NWChem are named ffield_ℓ where ℓ is one of: s x q u t c
standard parameters as published for the force field extensions as published in open literature contributed parameters by NWChem team user contributed parameters temporary current
Defined in input file or ~/.nwchemrc ffield amber amber_1 /software/nwchem/share/data/amber/amber_s/ amber_2 /software/nwchem/share/data/amber/amber_x/ amber_3 /home/newton/data/amber/amber_q/ spce /software/nwchem/share/data/solvents/spce.rst
File Names
system.ext system ext
user defined molecular system name file type, e.g. pdb PDB file top topology file seq sequence file
system_calc.ext system calc ext
user defined molecular system name user defined identification for the calculation identification, e.g. em, md002, tiA file type, e.g. out output file rst restart file qrs energy minimized restart file prp property file trj trajectory file gib free energy file
PREPARE Functionality and requirements frg seq
sgm top
pdb
rst
Functionality:
Requirements
•Topology and restart file generation • Coordinates from pdb formatted file or geometry input • Fragment and segment file generation from coordinates • Solvation • Potential of mean force functions • Topology modification for free energy and QM/MM calculations • File format conversion, e.g. from rst to pdb
• PDB format, i.e. IUPAC atom names, residue names, etc. • Automated atom typing based on force field typing rules • Force field parameters from par file(s) • Partial atomic charges from frg or sgm files
PREPARE Input Example PDB format, i.e. IUPAC atom names, residue names, etc. Automated atom typing based on force field typing rules Force field parameters from par file(s) Partial atomic charges from frg or sgm files memory noverify heap 1 mb stack 48 mb global 24 mb start crown basis "ao basis" print H library sto-3g Define basis sets C library sto-3g calculated O library sto-3g Na library sto-3g end prepare system crown_em modify atom 7:Na set 3 type K end task prepare
in case partial charges need to be
Read coordinates from crown.pdb Specify atom type change in topology Generate topology file crown.top Restart file crown_em.rst
MD Input Example md system crown_md data 1000 isotherm isobar record rest 1000 coord 100 prop 10 end task md dynamics task shell “cp crown_md.rst crown_ti.rst” md system crown_ti equil 1000 data 2000 over 1000 step 0.002 isotherm isobar new forward 21 of 21 print step 100 stat 1000 record rest 1000 free 1 end task md thermodynamics
Molecular dynamics input: NpT ensemble, 1000 steps
Copy restart file.
Free energy simulation
MCTI: 21*(1000+2000) NpT ensemble
ANALYZE input memory heap 8 mb stack 64 mb global 24 mb start job analysis system system_calc reference system_calc.rst file crown_md?.trj 0 10 select _C? _O? essential project 1 crown_vec1 project 2 crown_vec2 end
# ? Is wild card, replaced by series 0-10 # select subset of atoms # Perform essential dynamics analysis Projection onto specified vector
task analyze
Files required:
system.top ( molecular topology ) system_calc.rst ( coordinates )
Relative solvation ∆G Input Example memory noverify heap 10 mb stack 32 mb global 32 mb start crown basis "ao basis" print H library sto-3g C library sto-3g O library sto-3g Na library sto-3g end prepare system crown_em modify atom 7:Na set 3 type K end task prepare md system crown_em sd 1000 noshake end task md optimize prepare read crown_em.qrs write crown_em.pdb solvate write crown_md.rst end task prepare
Define basis sets
Read coordinates from crown.pdb Specify atom type change in topology Generate topology file crown.top Restart file crown_em.rst
Perform energy minimization
Solvate the solute Generate crown_md.rst
Relative solvation ∆G Input Example (continued) md system crown_md data 1000 isotherm 298.15 trelax 0.1 0.1 isobar record rest 1000 scoor 100 prop 10 end task md dynamics task shell “cp crown_md.rst crown_ti.rst” md system crown_ti equil 1000 data 2000 over 1000 step 0.002 isotherm 298.15 trelax 0.1 0.1 isobar new forward 21 of 21 print step 100 stat 1000 record rest 100 free 1 end task md thermodynamics
Molecular dynamics input: NpT ensemble, 1000 steps
Copy restart file.
Free energy simulation input:
MCTI: 21*(1000+2000) NpT ensemble
QM/MM approach
Combines two different descriptions – quantummechanical and classical
Cytochrome IFC3, Schewanella
The level theory changes based on a particular region Reactive regions – quantum mechanical description (QM) Regions where no chemical changes occur (or important) are treated a the classical molecular mechanics level (MM) 14
QM Region
109 atoms, DFT
MM region, 80,000 atoms
Structure of QM/MM energy functional
Eqm / mm (r, R;ψ ) = Eqm (r, R;ψ ) + Emm (r, R )
All QM-dependencies are in the first term
int ext E= r , R ; ψ E r , R ; ψ + E [ ] [ ] qm qm qm [r , R; ρ ]
Internal QM energy (theory dependent)
∑∫ I
All other classical terms are in the second term
15
Coulomb interactions with MM atoms
Bonded, angle, dihedral Coulomb interactions Vdw interactions
Z I ρ ( r ') RI − r '
dr '
QM/MM Interface in NWChem MM QM
QM/MM
Modular implementation sits on top of molecular mechanics (MM) and quantum (QM) modules Generic interface driven by function calls Manages data flow domain decomposition of coordinates in MM module replicated geometry data in QM module all data transfers happen in core Dispatches high level operations (e.g. optimization)
16
Summary of QM/MM Capabilities
Ground and excited state properties Structural optimization Reaction Pathway Calculations Dynamical Simulations Statistical Sampling (free energies)
17
Example of QM/MM Single Point Calculation start capk charge 2
QM region definition
prepare source capk.pdb modify segment 358 quantum … ignore update lists write capk.rst end task prepare
Active Site of cAPK kinase of QM/MM DFT calculation 75Scaling atoms, 817 basis 300 functions 250
Time (s)
200
md system capk cutoff 1.5 end
100 50
basis * library "Ahlrichs pVDZ" end qmmm link_atoms hydrogen end
150
0
Hydrogen Link Atoms
8
16
Commence QM/MM DFT Calculation task qmmm dft energy Valiev, M., et al.JPC B, 2007. 111(47): p. 13455-13464. 18
32
64
Processors
96
128
Example: QM/MM Excited State Calculations of cytosine base in DNA
• QM/MM
coupled-cluster CR-EOMCCSD(T) calculations of two lowest excited states • Protein environment has a significant influence on the excitation leading to a 0.4 eV stabilization of the ππ* excited state compared to gas phase • M. Valiev, K. Kowalski, JCP, 125(21), (2006) 19
QM/MM Optimization Large system sizes (103-106) makes direct optimizations impractical even with QM/MM approximation Key observations Most of degrees of freedom are in MM region The structure of far away MM regions has a little effect on the structure of QM region Small displacements of MM atoms affect little the electronic structure of QM region.
Decouple optimization of QM and MM regions
20
Optimization Algorithm 1. QM Region Optimization a. b. c.
MM region is fixed Typically 20-30 steps Requires solution of Schrödinger Equation
2. MM Region Optimization a. b. c.
QM region is fixed 1000-3000 steps QM region is represented either as effective point charges or static electron density distribution
3. Repeat the cycle until convergence
21
Example: Optimization of Zn-porphyrin in solution
QM region Zn-porphyrin (37 atoms) DFT/B3LYP MM region – 869 SPC/E waters 4.5 hours on 48 processors versus direct optimization would take ~ 2 days
22
…. qmmm region qm solvent maxiter 10 3000 method lbfgs sd ncycles 20 density espfit end task qmmm dft optimize
Input File
Reaction Pathway II
Nudged Elastic Band Method Pathway approximated discrete set of intermediate structures Beads represent different snapshots of reactive QM region along the pathway Forces on beads are calculated at the relaxed solvent configuration
23
Free Energy Differences 1 − β E ( A) − E ( B ) ) ∆WAB = − ln e (
β
A
B A
Accurate quantum mechanical description is a major challenge especially for high level methods ( 104 – 105 energy evaluations) The solution Introduce intermediate less expensive representation(s) Redistribute sampling using thermodynamic cycles
24
QM/MM Representations Different QM/MM representations Computational Efficiency
CC/MM
DFT/MM
MM/MM
Complexity and Accuracy
Example of MM/MM representation - QM atoms are replaced by effective point charges Qi reproducing correct field Eqm = ∑ i ,I
25
Z I Qi R I − ri
Z I ρ (r ') Z I Qi ∑I ∫ R − r' dr' = ∑ i , I R I − ri I
Free Energy Ladder 1 − β E ( A) − E ( B ) ) ∆WAB = − ln e (
Representation
β
A
B A
MM/MM
EMM/MM(A)
EMM/MM(B)
DFT/MM
EDFT/MM(A)
EDFT/MM(B)
A
B
Configuration
DFT →MM DFT →MM MM − ∆W BB + ∆W ∆W AB = (∆W AA ) AB
Valiev et al JCP 127, 051102 (2007) 26
Free Energy Ladder 1 − β E ( A) − E ( B ) ) ∆WAB = − ln e (
Representation
β
A
B A
MM/MM
EMM/MM(A)
EMM/MM(B)
DFT/MM
EDFT/MM(A)
EDFT/MM(B)
CC/MM
ECC/MM(A)
ECC/MM(B)
A
B
Configuration
MM DFT →MM DFT →MM CC →DFT CC →DFT + ∆W − ∆W + ∆W − ∆W BB ∆W AB = (∆W AA ) AB ) ( AA BB
27
Calculation of MM/MM free energy EMM/MM(A)
EMM/MM(B)
Can use any of the methods developed for classical free energy calculations Transformation between A and B configurations rλ = (1− λ )rA + λrB Qλ = (1− λ )QA + λQB
Free Energy Perturbation ∆W
ESP AB
= −∑ i
28
1
β
ln e
− β∆EλESP i →λi +1
λi
Calculation of DFT->MM free energy EMM/MM(A) EDFT/MM(A)
∆W
DFT → MM AA
=−
1
β
ln e (
DFT →MM − β ∆E AA
) MM / MM
DFT → MM ∆E AA = EDFT / MM (rA , R;ψ A ) − EMM / MM (rA , R; QA )
“Vertical” change of transformation (fixed QM region) MM representation is closely tailored to DFT by point charge fitting Approximate by the energy difference Can utilize free energy perturbation approach by resampling MM/MM trajectory
29
Example : Free energy barrier for the phosphorylation reaction in kinase protein
~15 kcal/mol reaction barrier
• cAPK protein kinase catalyze the transfer of the g-phosphoryl • Determination of the reaction pathway using NEB QM/MM ap • Calculation of free energy using effective charge approximatio • Valiev et al J. Phys Chem B. 111(47):13455-64. (2007). 30
Importance of free energy (Cheng et al JACS, 127,1553,2005)
PKA
31
Questions?
32
