1 NWChem: Molecular Dynamics and QM/MM2 Molecular Dynamics Functionality Target systems: General features: biomolecules (proteins, DNA/RNA, biomembran...

Author:
Jayson Wilkinson

2 downloads 70 Views 2MB Size

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

Our partners will collect data and use cookies for ad personalization and measurement. Learn how we and our ad partner Google, collect and use data. Agree & Close