280 likes | 546 Views
Friday of Week 3. Monday: force fields and Verlet integration Building a bpti simulation Wednesday: Review of sims setup on dropbox Water,salt-water,bpti Running namd and analysis Friday: Tom's Review/Overview Student reports on progress. 5E Discovery Approach.
E N D
Friday of Week 3 • Monday: • force fields and Verlet integration • Building a bpti simulation • Wednesday: • Review of sims setup on dropbox • Water,salt-water,bpti • Running namd and analysis • Friday: • Tom's Review/Overview • Student reports on progress
5E Discovery Approach http://www.nasa.gov/audience/foreducators/nasaeclips/5eteachingmodels/index.html Engage: Demo VMD/NAMD usage w/out much background Explore: you "play" with the VMD/NAMD scripts I've provided MAKE use of the VMD and NAMD user's guides and Tutorials Explain: I present theory, science, algorithms underlying what NAMD/VMD are doing You explain with your presentations what you have done/learned Elaborate/Extend: Make an hypothesis and test it. Find a feature of VMD/NAMD... show and tell. Evaluate: Present results to class AND discuss.
Organization of Sims Directory bin bpti salt-water water
Updated Salt-Water 01-alt-make-saltwater.vmd display projection Orthographic package require solvate solvate -o OUT/tmp -minmax {{ -20 -20 -20 } { 20 20 20 } } package require autoionize autoionize -psf OUT/tmp.psf -pdb OUT/tmp.pdb -o OUT/saltwater -sc .150 set all [atomselect [molinfo top ] " all " ] $all moveby [ vecscale -1 [measure center $all ] ] set minmax [ measure minmax $all ] set box [vecsub [lindex $minmax 1] [lindex $minmax 0]] molinfo [molinfo top] set {a b c alpha beta gamma} "$box 90.0 90.0 90.0" animate write pdb {OUT/ion.pdb} beg 0 end 0 skip 1 [molinfo top] mol delete 0 set cellfp [ open "OUT/saltwater.xsc" w ] puts -nonewline $cellfp [format "0 " ] puts -nonewline $cellfp [format "%12.2f %12.2f %12.2f " [lindex $box 0 ] 0 0 ] puts -nonewline $cellfp [format "%12.2f %12.2f %12.2f " 0 [lindex $box 1 ] 0 ] puts -nonewline $cellfp [format "%12.2f %12.2f %12.2f " 0 0 [lindex $box 0 ] ] puts -nonewline $cellfp [format "%12.2f %12.2f %12.2f " 0 0 0 ] close $cellfp quit
Molecular Dynamics Papers See the Doc-and-Man pdf files for Amber-Case2005.pdf Charmm-Brooks2009.pdf Namd-Phillips2005.pdf These provide technical details and references.
Molecular Dynamics Overview Potential Energy: Force: Dynamics: Numerical Integration: Velocity Verlet
Molecular Dynamics in Practice 1) Build the system 2) minimize the system Goal here is eliminate forces that break the integrator 3) equilibrate the system Goal here is to obtain correct N,T,P values 4) collect dynamics/sampling Goal here is to do your science
# input set SYS water coordinates OUT/$SYS.pdb extendedSystem OUT/$SYS.xsc structure OUT/$SYS.psf parameters IN/par_all22_prot.inp paratypecharmm on # output set output OUT/min-eq outputname $output dcdfile ${output}.dcd veldcdfile ${output}.dvd xstFile ${output}.xst veldcdfreq 500 dcdfreq 500 xstFreq 500 binaryoutput yes binaryrestart yes outputEnergies 100 restartfreq 0 fixedAtoms off # Basic dynamics exclude scaled1-4 1-4scaling 1 COMmotion no dielectric 1.0 # Simulation space partitioning switching on switchdist 9 cutoff 10 pairlistdist 12 # Multiple timestepping firsttimestep 0 timestep 1 stepspercycle 20 nonbondedFreq 2 fullElectFrequency 4 # Temperature control set temperature 298 temperature $temperature; # initial temperature # PBC wrapAll on dcdUnitCell yes PME yes PMEGridSpacing 1 # Scripting minimize 1000 reinitvels $temperature run 10000 Let's Walk Thru a NAMD Config File
Inputs set SYS water coordinates OUT/$SYS.pdb extendedSystem OUT/$SYS.xsc structure OUT/$SYS.psf parameters IN/par_all22_prot.inp paratypecharmm on
Outputs set output OUT/min-eq outputname $output dcdfile ${output}.dcd veldcdfile ${output}.dvd xstFile ${output}.xst veldcdfreq 500 dcdfreq 500 xstFreq 500 binaryoutput yes binaryrestart yes outputEnergies 100 restartfreq 0
Basic Dynamics exclude scaled1-4 1-4scaling 1 COMmotion no dielectric 1.0 switching on switchdist 9 cutoff 10 pairlistdist 12 set temperature 298 temperature $temperature; # initial temperature
Cut-off and Switching Long range Lennard Jones Correction and be used to fill in the long range part
Long-Range Approximations PME: particle mess Ewald FullDirect: Fast Multipole Approximation:
PBC and PME # PBC extendedSystem OUT/$SYS.xsc wrapAll on dcdUnitCell yes PME yes PMEGridSpacing 1
Multiple Time Stepping # Multiple timestepping firsttimestep 0 timestep 1 stepspercycle 20 nonbondedFreq 2 fullElectFrequency 4
Water and Hydrogens Special algorithms are used for water and for hydrogens since “H” is the lightest (read fastest moving atom) WATER: make it rigid → bonds and angles fixed Settle or Shake algorithms. Idea move waters then correct positions. Other Bonds: move heavy atom then “fix” hydrogen
NVE No special switches this is the default
NVT langevin on; # do langevin dynamics langevinDamping 1; # damping coefficient (gamma) of 1/ps langevinTemp $temperature; # bath temperature langevinHydrogen no; # don't couple langevin bath to hydrogens seed 12345 WIKIPEDIA EVEN WARNS:If the main objective is to control temperature, care should be exercised to use a small damping constant g. As g grows, it spans the inertial all the way to the diffusive (Brownian) regime. The Langevin dynamics limit of non-inertia is commonly described as Brownian dynamics.
NTP = NVToptions + P control(scale the box size) # Pressure control langevinPiston on langevinPistonTarget 1.01325; # in bar -> 1.01325 bar = 1 atm langevinPistonPeriod 200 langevinPistonDecay 100 langevinPistonTemp $temperature useFlexibleCell no useGroupPressure no fixedAtomsForces off
Velocity Verlet Algorithm Derivative Taylor Verlet Just add: Rearrange: Use Physics:
Velocity Verlet in Practice Half-kick Drift Force eval Half-kick You can also add extra forces, restrict the maximum move, modify velocities (e.g. rescale) etc... At the various steps as needed
RMSD analysis N is the number of atoms being compared. Nt is the number of time steps being used in the calculation. r(tj) is the position of atom at time t and <r> is the average value of the position of atom . The RMSD should converge if the atoms are not “flowing” but rather fluctuating around some position.
Boltzman H-function from Velocities in NVT (canonical) Fraction of atoms with velocity between v and v + dv For any component For any function A that depends on velocities (e.g. K.E.) Boltzmann's H-function as measure equilibration. Hx is negative and should decrease in time to by f(vx)dvx above “Molecular Dynamics Simulations, Elem.Methods”, J.M.Haile, Chap 2.7
Diffusion Green-Kubo expression for integrated velocity autocorrelation function. (see McQuarrie Stat. Mechanics) Or the usual expression
Backcalculating a Force Constant • If the system is equilibrated then equipartition holds • Each oscillator has 1/2kT of energy so for all bonds which have 1/2Ko(R-Ro)^2 then nbonds/2kT = 1/2Ko<R-R0>^2 or <R-Ro>^2*kT/nbonds = Ko NOTE in a typical system not all bonds have same value Ko so you must be careful to select only one type of bond in you analysis.
Some Exercises to Consider 0)time evolution of RMSD of bpti backbone differs depending on ensemble used: NVE,NTP, NVT... do they converge to the same value after sufficient simulation time? what differs in the structure between bpti in the systems? 1) check that radial distribution function in simulation of water gives proper values for gO-O and gO-H and compare to values in literature. 2) check that radial distribution function for IONS in the salt-water sims is correct and compare to values in literature. Note you can use the "Ioninize" plug in to check different ions and different concentrations of ion. 1a and 2a ) instead of radial distribution look at values of Diffusion observed in the systems. Does presence of salt alter the diffusion of waters in the simulations? 1b and 2b) how do simulation parameters, namely langevin damping term, alter the observed values of diffusion? 3) Performance tests: benchmarks of timing vs. system size for various size boxes of water... how does performance scale with number of cpus? how does performance sale with system size? what simulation options most strongly affect simulation speed? 4) benchmarks of one of the systems showing how to optimization a simulation with GPU enabled version of NAMD 5) vary the Integration time step and see how conserved quantities e.g total energy behaves in the NVE ensemble, are affected by the timestep. Can you prove that the integrator is 4th order by varying the integration step? What's max value of dt that can use before the system fails. 6) how long does a system of just water take to converge to the correct density if you start off with in correct boundary conditions? what is the best minimization and equilibration protocol for "fixing" it? 7) any of the above as a function of simulation parameter? note that VMD/NAMD come with a couple of different flavors of the charmm force field so you can even see if differences in force field affect your answer. (first check that for your particular system there are differences in force field parameters)