1、分子动力学原理、算法简介与与G Gromacs的使用的使用部分内容选自超算中心2012年Gromacs培训教材计算化学虚拟实验室（VLCC）计算化学虚拟实验室（VLCC）高性能计算化学应用的新型虚拟科研组织由网络中心和大连化物所于2004年发起成立联合国内外知名计算化学科学家的平台联合国内外知名计算化学科学家的平台高性能计算化学与eScience应用研究高性能计算化学应用促进高性能计算化学应用促进 研发SCChem系统、ICTHPCC国际会议、用户培训 应用技术支持VLCC能为您的工作提供哪些帮助？1 1.提供版计算化学软件提供正版计算化学软件2 2.计算化学高性能计算的解决方案计算化学高性能
2、计算的解决方案最近的工作蛋白质折叠的新型高精度模拟软件的开发和移植最近的工作：蛋白质折叠的新型高精度模拟软件的开发和移植AMNW在天河2节点上实现近30的并行执行，并行效率超过过30%30%包含106个氨基酸的蛋白质体系，3392个子任务 在在142秒内实现5的蛋白质折叠计算模拟142秒内实现5ps的蛋白质折叠计算模拟3 3.计算程序的GPU/MIC并行移植计算程序的GPU/MIC并行移植最近的工作：最近的工作：1 1.耗散粒子动力学程序的GPU移植耗散粒子动力学程序的GPU移植开发了新的大规模计算算法 实实现了相比于单CPU核约50倍的加速现了相比于单CPU核约50倍的加速2 2.计算化学软
3、件的MIC上编译实现计算化学软件的MIC上编译实现NWchem AmberAmber大家有计算化学方面的软件，硬件或合大家有计算化学方面的软件，硬件或合作问题，请联系我们：Email:微博：主要内容主要内容一：分子动力学模拟原理二二：分子动力学模拟算法分子动力学模拟算法三：Gromacs使用初步三：Gromacs使用初步一：分子动力学模拟原理1.1分子动力学模拟及其应用1 1.2 2分子动力学基本方程分子动力学基本方程1.3 分子力场1.3分子力场1.4分子模拟的一般流程1.1分子动力学模拟及其应用Why MD simulations?Model phenomena that cannot b
4、e observed experimentallyAccess to thermodynamics quantities(free energies,binding energies,)Simulators and experimentalists can have a synergistic relationship,leading to new insights into materials propertiesFirstmoleculardynamicssimulation(1957/59)Hard disks and spheresr du(r)0 (calculation of co
5、llision times)uij(r)=rdIBM-704:solid phaseliquid phaseliquid-vapour-phaseN=32:7000 collisions/hProduction run 20000 stepsN=32 6.5x105coll.4 daysN=500:500 collisions/hProduction run 20000 steps N=500 107coll.2.3 years今今天的分子动力学模拟天的分子动力学模拟生命科学纳米材料凝聚态物理如：水结晶，相分离（如：蛋白质折叠）如：分子马达药物设计2013年诺贝尔化学奖年诺贝尔化学奖将实验带入
6、信息时代MartinKarplusMichaelLevittAriehWarshel复杂化学体系多尺度模型的发展算模结将实带信代量化计算与分子模拟的结合*所有图片来源于www.nobelprize.orgWh td WhSli MD Sil tiWhatandWhere:ScalesinMD SimulationslmesoscalecontinuummolecularMonte Carlole10-6Sexp(-E/kT)domainquantumchemistrymoleculardynamicsme Scal10-8S=F=MATim10-12S=Length Scale10-10m10
7、-8m10-6m10-4m15Length Scale1.2分子动力学基本原理Newtons equations of motionFi1.Ndtd rm2ii2i=Momentary force:dtF=U(rr)Acceleration tells us the speed.Fi=riU(r1,.,rN)Speed variation helps us to approximate positions after short times,this is called integrating the equations of motionCalculation for a huge numb
8、er of small steps results in a trajectory of particlesU(rr)U(r1,rN)(very important)1.3 分子力场1.Bonded:covalent bond-stretching,angle-bending,and 1.3分子力场dihedrals.These are computed on the basis of fixed lists.2.Non-bonded:Lennard-Jones or Buckingham,and Coulomb or modified Coulomb.The nonbondedinterac
9、tions are computed on the basis of a neighbor list(a list of non-bonded atoms within a certain radius),in which exclusions are already removed.3.Restraints:position restraints,angle restraints,distance restraints,orientation restraints and dihedral restraints,all 17based on fixed lists.TheForceField
10、bond stretchtorsionalNon-bondedinteractionsvalence anglebend181.3.1 Bonded InteractionsUbondedUbondstretchUbondbendUrotatealong gbondTorsion anglesAre 4-bodyNon-bondedypairpairAnglesA3 b dBondsAre 3-bodyA Are 2 2-b bod dy19Bonded Interactions:Bond StretchingHarmonic potentialbHarmonic potentialFourt
11、h p power p potentialMorse potential bond stretchingCubic bond stretching potentialIn the GROMOS-96 force field,(4.35)is used for reasons of computational efficiencyIn(4.44),it should be noted that the potential is asymmetric:overstretching leads to infinitely low energies Thefor reasons of computat
12、ional efficiency.overstretching leads to infinitely low energies.The integration time step is therefore limited to 1 fs.1()2Ubdtt h=KbbbKb=Force Constantb=Ideal Bond Lengthbondstretch2b(o)pairsUKbb20bo=Ideal Bond LengthBondedInteractions:dihedral angle The angle between planes(i,j,k)and(j,k,l)Figure
13、 4.8:Principle of improper dihedral angles.Out of plane bending for rings(left),substituents of rings(middle),out of tetrahedral(right).Improper dihedralsFi4 9 Idih d li lNote that in the input in topology files,angles are given in degrees and force constants in kJ/mol/rad2.21Figure 4.9:Improper dih
14、edral potential.1.3.2NonBondedInteractions:vanderWaalsThe Lennard-Jones interactionTwo frequently used formulas:(12)(6)ijijVDWCCE=(4 3)EVDW=4(ij)12(ij)6)(4 5)22VDW126ijijErrVDW()()rijrij(4.3)(4.5)Steric effects arise from the fact that each atom within a molecule occupies a certain amount of space.I
15、f atoms are brought too close together,there is an associated cost in energy due to overlapping electron clouds(Pauli or Born repulsion),and this may affect the molecules preferred shape(conformation)and reactivity.NonBondedInteractions:vanderWaalsBuckingham potentialFigure 4.2:The Buckingham intera
16、ction.23Non-Bonded Interactions:CoulombThe Coulomb interaction between two charge particles is given by:Vc(ij)fqi iqj jr ijq qVrf=r1f=138 935485kJ mol1nm e2(4.11)0138.935485f=4=kJ molnm e(4.11)1.4分子模拟的一般流程U(R)UbondedUnon bonded()bondednon bonded(Esteric energy=Estretch+Ebend+Eimproper+Etorsion+EvdW+
17、Eqq)Firi iU(r1,.,rN)(Potential function force)Fi1.Ndtd rm2ii2idt=(Newtons Law)25二：分子动力学模拟算法2.1非键接力的计算2 2.1 1.1 1短程相互作用计算短程相互作用计算2.1.2静电相互作用2.2积分、压力和温度控制2 2.3 3能量最小化能量最小化2.1.1短程相互作用计算Short-range force calculation save CPU time by using neighbor lists neighbor list:neighbor list:27 cell list:most effi
18、cient is a combination of Verlet and cell list.(8 neighbor cells for 2D 26 neighbor cells for 3D)update when a particle has for 2D,26 neighbor cells for 3D.)made a displacementmade a displacement (rskin-rcut)/2大规模并行策略：Domain decompositionPartition space,instead of atoms,over nodesGood for load balan
19、cing(uniform system,what about non-uniform system?)Bad for communication bandwidthEach node imports coordinate and exports forces from neighbors within a sphere with radius=cutoff.2.1.2 静电相互作用周期性边界条件下静电相互作用的计算2.1.2静电相互作用11NNqiqjUcoul=40211|+|=xyzcoulnnnijijnxnynz=+rnnxyz算法复杂度：N2Nonbonded Interaction
20、s:Ewald SumEwald Sum近程相互作用远程相互作用FFT is a finite and discrete Fourier transformation th s point chargesith关联相互作用Standard Ewald implementation scales like N2Most costly part:Fourier transformationtransformation,thus point charges with continuous coordiates have to be replaced by a grid based charge de
21、nsity(step 1).Process can be speeded up by a method called“Fast Fourier Transformation”(FFT)Particle-Mesh-Ewald算法B-样条插值法扩散电荷无法显示图像。计算机可能没有足够的内存以打开该图像，也可能是该图像已损坏。请重新启动计算机，然后重新打开该文件。如果仍然显示红色“x”，则可能需要删除该图像，然后重新将其插入。样条插值法扩散电荷PME的并行策略：2.2 Numeric Integration MethodVerlets NumericIntegrationMethodTaylor e
22、xpansion:Verlets Method*Verlet.,L.Computer experiments on classical fluids.i.thermodynamical properties of lennard-jones molecules.Phys.Rev.159:98103,1934LeapfrogNumericIntegrationMethodThe MD program utilizes the so-called leap-frogalgorithm*for the integration of the equations of motion.The leap-f
23、rog algorithm uses positions rat time tand velocities v vat timett/2;it updates positions and velocities using the forcesF(t)determinedat time t-t/2;it updates positions and velocities using the forces F(t)determined by the positions at time t:Figure 3.6:The Leap-Frog integration method.The algorith
24、m is called Leap-Frog because rIt is equivalent to the Verlet algorithm.The algorithm is of third order in rand is time-reversible.and vare leaping like frogs over each others back.time reversible.2.2TemperaturecontrolschemesBerendsen temperature couplingThe Berendsen algorithm mimics weak coupling
25、with first-order kinetics to an external heat bath with given temperature T0.dTT0Tdt=gp0(23)The Berendsen thermostat suppresses the fluctuations of the kinetic energy.This means that,strictly speaking,one does not generate a proper canonical ensemble.For very small systems the sampling will indeed b
26、e incorrect.But for larger systems most properties will not be affected significantly,except for the distribution of the kinetic energy itself.A similar thermostat which does p produce a correct ensemble is the velocity y rescaling g thermostat described below.Velocity rescaling thermostatThe veloci
27、ty rescaling thermostat is essentially a Berendsen thermostat(see above)withThe velocity rescaling thermostat is essentially a Berendsen thermostat(see above)with an additional stochastic term which ensures a correct kinetic energy distribution:where Kis the kinetic energy,Nfthe number of degrees of
28、 freedom and dWa Wiener process.There are no additional parameters,except for a random seed.This thermostat produces a correct canonical ensemble and still has the advantage of the Berendsen thermostat:first order decay of temperature deviations and no oscillations.When an NVT ensemble is used,the c
29、onserved energy quantity is written to the energy and log file.Berendsen et al.Molecular dynamics with coupling to an external bath.J.Chem.Phys.81:3684(1984)2.2TemperaturecontrolschemesNose-Hoover(chain)temperature couplingThe Berendsen weak coup pling g alg gorithm is extremely y efficient for rela
30、xing g a sy ystem to the target temperature,but once your system has reached equilibrium it might be more important to probe a correct canonical ensemble.This is unfortunately not the case for the weak coupling scheme,although the difference is usually negligible.22iiiidddt=mdtrFr(3.27)01d(TT)dtQ=(3
31、.28)The mass parameter is a some hat a k arda of describing co pling strength especiall d e to its dependence onThe mass parameter is a somewhat awkward way of describing coupling strength,especially due to its dependence on reference temperature(and some implementations even include the number of d
32、egrees of freedom in your system when defining Q).To maintain the coupling strength,one would have to change Qin proportion to the change in reference temperature.For this reason,we prefer to let the GROMACS user work instead with the period Tof the oscillations of(3.29)kinetic energy between the sy
33、stem and the reservoir instead.It is directly related to Qand T0viaThis provides a much more intuitive way of selecting the Nose-Hoover coupling strength(similar to the weak coupling relaxation),and in additionTis independent of system size and reference temperature.2.2TemperaturecontrolschemesIt is
34、 however important to keep the difference between the weak coupling scheme and the Nose-Hoover algorithm in mind:Using weak coupling you get a strongly damped exponential relaxation,while the Nose-Hoover approach produces an oscillatory relaxation.The actual time it takes to relax with Nose-Hoover c
35、oupling is several times larger than the period of the oscillations that you select.These oscillations(in contrast to exponential relaxation)also means that the timeconstant normally should be 45 times larger than the relaxation time used with weak couplibilling,but your mileage may vary.Group tempe
36、rature couplingIn GROMACS temperature coupling can be performed on groups of atoms typically a protein andIn GROMACS temperature coupling can be performed on groups of atoms,typically a protein and solvent.The reason such algorithmes were introduced is that energy exchange between different componen
37、ts is not perfect,due to different effects including cutoffs etc.If now the whole system is coupled to one heat bath,water(which experiences the largest cutoff noise)will tend to heat up and the protein will cool down.Typically 100 K differences can be obtained.With the use of proper electrostatic m
38、ethods(PME)these difference are much smaller but still not negligable.The parameters for temperature coupling in groups are given in the mdp file.One special case should be mentioned:it is possible to T-couple only part of the system(or nothing at all obviously)This is done by specifyingis possible
39、to T couple only part of the system(or nothing at all obviously).This is done by specifying zero for the time constant Tfor the group of interest.2.2 Simulating at constant PIn the same spirit as the temperature coupling,the system can also be coupled to a pressure bath.GROMACS supports both the Ber
40、endsen algorithm*that scales coordinates and box vectors every step,and the extended ensemble Parrinello-Rahman approach.Both of these can be combined with any of the temperature coupling methods above.Berendsen pressure couplingTh Bdlithlthditd bttitht iThe Berendsen algorithm rescales the coordina
41、tes and box vectors every step with a matrix,which has the effect of a first-order kinetic relaxation of the pressure towards a given reference pressure P0:The scaling matrix is given byHere is the isothermal comp pressibility y of the sy ystem.In most cases this will be a diagonal matrix,with equal
42、 elements on the diagonal,the value of which is generally not known.For water at 1 atm and 300 K:=4 610-10Pa-1=4 610-5Bar-1Berendsen et al.Molecular dynamics with coupling to an external bath.J.Chem.Phys.81:3684(1984)For water at 1 atm and 300 K:=4.61010Pa1=4.6105Bar12.2 Simulating at constant PParr
43、inello Rahman pressure couplingParrinello-Rahman pressure couplingIn cases where the fluctuations in pressure or volume are important per se(e.g.to calculate thermodynamic properties)it might at least theoretically be a problem that the exact ensembl ill d fid fhklihble is not well-defined for the w
44、eak coupling scheme.For this reason,GROMACS also supports constant-pressure simulations using the Parrinello-Rahman approach*,which is similar to the Nose-Hoover temperature coupling.WWith th PillR hbt t th btt d b tht ibbith the Parrinello-Rahman barostat,the box vectors as represented by the matri
45、x bobey the matrix equation of motion.The volume of the box is denoted V,and Wis a matrix parameter that determines the streng gth of the coup pling g.The matrices Pand Prefare the current and reference pressures,refp,respectively.The(inverse)mass parameter matrix W-1determines the strength of the c
46、oupling,and how the box can be deformedthe box can be deformed.Yl hidhiihlibili id h40You only have to provide the approximate isothermal compressibilities and the pressure time constant pin the input file(Lis the largest box matrix element)2.2 Simulating at constant PThe equations of motion for the
47、 particles are also changed,just as for the Nose-Hoover coupling.In most cases you would combine the Parrinello-Rahman barostat with the Nose-Hoover thermostat,but to keep it simple we only show the Parrinello-RahmanNose Hoover thermostat,but to keep it simple we only show the Parrinello Rahmanmodif
48、ication here:Just as for the Nose-Hoover thermostat you should realize that the Parrinello-Rahman timeJust as for the Nose-Hoover thermostat,you should realize that the Parrinello-Rahman time constant is notequivalent to the relaxation time used in the Berendsen pressure coupling algorithm.I In most
49、 t cases you willd t4 5 tiltittith PillR hill need to use a 45 times larger time constant with Parrinello-Rahmancoupling.If your pressure is very far from equilibrium,the Parrinello-Rahman coupling may result in very large box oscillations that could even crash your run.In that case you would have t
50、o increase the time constant or(better)use the weak coupling scheme to reachwould have to increase the time constant,or(better)use the weak coupling scheme to reach the target pressure,and then switch to Parrinello-Rahman coupling once the system is in equilibrium.412.3 EnergyMolecular dynamicsuses
