分子动力学原理、算法简介与Gromacs的使用

1、分子动力学原理、算法简介与与GGromacs的使用的使用部分内容选自超算中心2012年Gromacs培训教材计算化学虚拟实验室（VLCC）计算化学虚拟实验室（VLCC）高性能计算化学应用的新型虚拟科研组织由网络中心和大连化物所于2004年发起成立联合国内外知名计算化学科学家的平台联合国内外知名计算化学科学家的平台高性能计算化学与eScience应用研究高性能计算化学应用促进高性能计算化学应用促进研发SCChem系统、ICTHPCC国际会议、用户培训应用技术支持VLCC能为您的工作提供哪些帮助？11.提供版计算化学软件提供正版计算化学软件22.计算化学高性能计算的解决方案计算化学高性能计算的解决

2、方案最近的工作蛋白质折叠的新型高精度模拟软件的开发和移植最近的工作：蛋白质折叠的新型高精度模拟软件的开发和移植AMNW在天河2节点上实现近30万核的并行执行，并行效率超过过30%30%包含106个氨基酸的蛋白质体系，3392个子任务在在142秒内实现5的蛋白质折叠计算模拟142秒内实现5ps的蛋白质折叠计算模拟33.计算程序的GPU/MIC并行移植计算程序的GPU/MIC并行移植最近的工作：最近的工作：11.耗散粒子动力学程序的GPU移植耗散粒子动力学程序的GPU移植开发了新的大规模计算算法实实现了相比于单CPU核约50倍的加速现了相比于单CPU核约50倍的加速22.计算化学软件的MIC上编译

3、实现计算化学软件的MIC上编译实现NWchemAmberAmber大家有计算化学方面的软件，硬件或合大家有计算化学方面的软件，硬件或合作问题，请联系我们：Email:微博：主要内容主要内容一：分子动力学模拟原理二二：分子动力学模拟算法分子动力学模拟算法三：Gromacs使用初步三：Gromacs使用初步一：分子动力学模拟原理1.1分子动力学模拟及其应用11.22分子动力学基本方程分子动力学基本方程1.3分子力场1.3分子力场1.4分子模拟的一般流程1.1分子动力学模拟及其应用Why MD simulations?Model phenomena that cannot be observed e

4、xperimentallyAccess 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 spheresrdu(r)0 (calculation of collision ti

5、mes)uij(r)= rdIBM-704:solid phaseliquid phaseliquid-vapour-phaseN=32: 7000 collisions / hProductionrun20000stepsN=32 6.5x105coll. 4 daysN=500: 500 collisions / hProduction run 20000 steps N=500 107coll. 2.3 years今今天的分子动力学模拟天的分子动力学模拟生命科学纳米材料凝聚态物理如：水结晶，相分离（如：蛋白质折叠）如：分子马达药物设计2013年诺贝尔化学奖年诺贝尔化学奖将实验带入信息时代

6、MartinKarplusMichaelLevittAriehWarshel复杂化学体系多尺度模型的发展算模结将实带信代量化计算与分子模拟的结合*所有图片来源于WhtdWhSliMDSiltiWhatandWhere:ScalesinMDSimulationslmesoscalecontinuummolecularMonte Carlole10-6Sexp(-E/kT)domainquantumchemistrymoleculardynamicsme Scal10-8S = F = MATim10-12S = LengthScale10-10m10-8m10

7、-6m10-4m15Length Scale1.2分子动力学基本原理Newtons equations of motionFi1.Ndtdrm2ii2i=Momentary force:dtF=U(rr)Acceleration tells us the speed. Fi=riU(r1,.,rN)Speed variation helps us to approximate positionsafter short times, this is called integrating the equations of motionCalculation for a huge number of

8、 small steps results in a trajectoryof 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 no

9、nbondedinteractions 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 fixe

10、d lists.TheForceFieldbond stretchtorsionalNon-bondedinteractionsvalence anglebend181.3.1 Bonded InteractionsUbonded=Ubondstretch+Ubondbend+UrotatealonggbondTorsion anglesAre 4-bodyNon-bondedypairpairAnglesA3bdBondsAre 3-bodyAAre 22-bboddy19Bonded Interactions: Bond StretchingHarmonicpotentialbHarmon

11、ic potentialFourth ppower ppotentialMorse potential bond stretchingCubic bond stretching potentialIn the GROMOS-96 force field, (4.35)is used forreasonsofcomputationalefficiencyIn (4.44), it should be noted that the potential is asymmetric: overstretchingleadstoinfinitelylowenergiesThefor reasons of

12、 computational efficiency.overstretching leads to infinitely low energies. The integration time step is therefore limited to 1 fs.1()2Ubdtth=KbbbKb= Force Constantb=IdealBondLengthbondstretch2b(o)pairsUKbb20bo= Ideal Bond LengthBondedInteractions:dihedral angle The angle between planes (i,j,k) and (

13、j,k,l)Figure 4.8: Principle of improper dihedral angles. Out of plane bending for rings(left), substituents of rings(middle), out of tetrahedral(right). Improper dihedralsFi49IdihdlilNotethat in the input in topology files, angles are given in degrees and force constants in kJ/mol/rad2.21Figure 4.9:

14、 Improper dihedral potential.1.3.2NonBondedInteractions:vanderWaalsThe Lennard-Jones interactionTwo frequently used formulas:(12)(6)ijijVDWCCE=(43)EVDW=4(ij)12(ij)6)(45)22VDW126ijijErrVDW()()rijrij(4.3)(4.5)Steric effectsarise from the fact that each atomwithin a moleculeoccupies a certain amount of

15、 space. If atoms are brought too close together, there is an associated cost in energydue to overlapping electron clouds(Paulior Bornrepulsion), and this may affect the molecules preferred shape (conformation)andreactivity.NonBondedInteractions:vanderWaalsBuckingham potentialFigure 4.2: The Buckingh

16、am interaction.23Non-Bonded Interactions: CoulombThe Coulomb interactionbetween two charge particles is given by:Vc(ij)fqiiqjjrijqqVrf=r1f=138935485kJmol1nme2(4.11)0138.935485f=4=kJ molnm e(4.11)1.4分子模拟的一般流程U(R)=Ubonded+Unonbonded()bondednonbonded(Esteric energy= Estretch+ Ebend+ Eimproper+ Etorsion

17、+ EvdW+ Eqq)Fi=riiU(r1,.,rN)(Potential function force)Fi1.Ndtdrm2ii2idt=(Newtons Law)25二：分子动力学模拟算法2.1非键接力的计算22.11.11短程相互作用计算短程相互作用计算2.1.2静电相互作用2.2积分、压力和温度控制22.33能量最小化能量最小化2.1.1短程相互作用计算Short-range force calculation save CPU time by using neighbor lists neighborlist:neighbor list:27cell list:most effi

18、cient is a combination of Verlet and cell list. (8 neighbor cells for2D26neighborcellsfor3D)update when a particle has for 2D, 26 neighbor cells for 3D.)madeadisplacementmade a displacement (rskin-rcut)/2大规模并行策略：Domain decompositionPartition space, instead of atoms, over nodesGood for load balancing

19、 (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算法复杂度：N2NonbondedInteractions: E

20、waldSumEwald Sum近程相互作用远程相互作用FFT is a finite and discrete Fourier transformationthspointchargesith关联相互作用Standard Ewald implementation scales like N2Most costly part: Fourier transformationtransformation, thus point charges with continuous coordiateshave to be replaced by a grid based charge density (

21、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 MethodVerletsNumericIntegrationMethodTaylor expans

22、ion:Verlets Method*Verlet., L. Computer experiments on classical fluids. i. thermodynamicalproperties 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 lea

23、p-frog algorithm uses positions rat time tand velocities vvattimett/2;itupdatespositionsandvelocitiesusingtheforcesF(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 algorithm

24、 is called Leap-Frog because rIt is equivalent to the Verletalgorithm. The algorithm is of third order in rand is time-reversible.and vare leaping like frogs over each others back.timereversible.2.2TemperaturecontrolschemesBerendsen temperature couplingThe Berendsen algorithm mimics weak couplingwit

25、h first-order kinetics to an external heat bath with given temperature T0.dTT0Tdt=gp0(23)The Berendsenthermostat 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 be

26、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 pproduce a correct ensemble is the velocityy rescalingg thermostat described below.Velocity rescaling thermostatThevelocityr

27、escalingthermostatisessentiallyaBerendsenthermostat(seeabove)withThe velocity rescaling thermostat is essentially a Berendsenthermostat (see above) with an additional stochastic term which ensures a correct kinetic energy distribution:where Kis the kinetic energy, Nfthe number of degrees of freedom

28、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 Berendsenthermostat: first order decay of temperature deviations and no oscillations.When an NVT ensemble is used, the conser

29、ved 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 Berendsenweak coupplingg alggorithm is extremelyy efficient for relaxi

30、ngg a syystem to the target temperature, butonce 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.2

31、8)ThemassparameterisasomehatakardaofdescribingcoplingstrengthespecialldetoitsdependenceonThemass 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 degrees of freedom in yo

32、ur 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 system and the rese

33、rvoir 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 however impor

34、tant 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 ittakes to relax with Nose-Hoover coupling is se

35、veral times largerthan the period of the oscillationsthat you select. These oscillations (in contrast to exponential relaxation) also means that the timeconstantnormally should be 45 times largerthan the relaxation time used with weak couplibilling,but your mileage may vary.Group temperature couplin

36、gInGROMACStemperaturecouplingcanbeperformedongroupsofatomstypicallyaproteinandIn GROMACS temperature coupling can be performed on groups of atoms, typically a protein and solvent. The reason such algorithmeswere introduced is that energy exchange between different components is not perfect, due to d

37、ifferent 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 methods (PME) these di

38、fference are much smaller but still not negligable. The parameters for temperature coupling in groups are given in the mdpfile. One special case should be mentioned: it ispossibletoT-coupleonlypartofthesystem(ornothingatallobviously)Thisisdonebyspecifyingis possible to Tcouple only part of the syste

39、m (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 Berendsenalgorithm* that scales c

40、oordinates and box vectors every step, and the extended ensemble Parrinello-Rahmanapproach. Both of these can be combined with any of the temperature coupling methods above.Berendsen pressure couplingThBdlithlthditdbttithtiThe Berendsenalgorithm rescales the coordinates and box vectors every step wi

41、th 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 comppressibilityy of the syystem. In most cases this will be a diagonal matrix, with equal elements on the diagonal, the va

42、lue of which is generally not known.Forwaterat1atmand300K:=4610-10Pa-1=4610-5Bar-1Berendsenet 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 PParrinelloRahmanpressurecouplingParri

43、nello-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 ensemblilldfidfhklihble is not well-defined for the weak coupling scheme.For this reason, G

44、ROMACS also supports constant-pressure simulations using the Parrinello-Rahmanapproach*, which is similar to the Nose-Hoover temperature coupling. WWiththPillRhbttthbttdbthtibbith the Parrinello-Rahmanbarostat, the box vectors as represented by the matrix bobey the matrix equation of motion.The volu

45、me of the box is denoted V, and Wis a matrix parameter that determines the strenggth of the coupplingg. The matrices Pand Prefare the current and reference pressures, refp,respectively.The (inverse) mass parameter matrix W-1determines the strength of the coupling, and how theboxcanbedeformedthe box

46、can be deformed.Ylhidhiihlibiliidh40You only have to provide the approximate isothermal compressibilitiesand 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 particles are also changed, just as for the Nose-

47、Hoover coupling. In most cases you would combine the Parrinello-Rahmanbarostatwith the Nose-Hooverthermostat,buttokeepitsimpleweonlyshowtheParrinello-RahmanNoseHoover thermostat, but to keep it simple we only show the ParrinelloRahmanmodification here:JustasfortheNose-Hooverthermostatyoushouldrealiz

48、ethattheParrinello-RahmantimeJust as for the Nose-Hoover thermostat, you should realize that the Parrinello-Rahmantime constant is notequivalent to the relaxation time used in the Berendsenpressure coupling algorithm.IIn mostt cases you willdt45tiltittithPillRhill need to use a 45 times larger time

49、constantwith Parrinello-Rahmancoupling. If your pressure is very far from equilibrium, the Parrinello-Rahmancoupling may result in very large box oscillationsthat could even crash your run. In that case you wouldhavetoincreasethetimeconstantor(better)usetheweakcouplingschemetoreachwould have to incr

50、ease the time constant, or (better) use the weak coupling scheme to reach the target pressure, and then switch to Parrinello-Rahmancoupling once the system is in equilibrium.412.3 能量最小化EnergyMolecular dynamicsuses thermal energyto explore the energysurffaceState AState BEnergyminimizationEnergy mini

51、mizationstops at local minimaSteepest descent (SD,最陡下降算法):The simplest iteration scheme consists of following the “steepestdescent” direction:xk+1=xkf(xk)(sets the minimumalong the line definedbythegradient)Usually, SD methods leads to improvement quickly, but then exhibitslow progress toward a solu

52、tion.Thlddfiitiliiitiittiby the gradient)They are commonly recommended for initial minimization iterations,when the starting function and gradient-norm values are very large.43Illustration of the steepest (gradient) descent method.DisadvantagesAs per the examples, the disadvantages of this method ar

53、e as follows:it is relatively slow close to the minimum;the linear search may cause problems;integrator = steep emtol= 100.0 ; in kJ mol-1 nm-1emstep= 0.002 ; in nm44it might zigzag down valleys.Conjugate gradients (CG):In each step of conjugate gradient methods, a search vector pkisdefined by a rec

54、ursive formula:()The corresponding new position is found by line minimization alongpk+1=f(xk)+k+1pkpk:xk+1=xk+kpkthe CG methods differ in their definition of :FletcherReeves:FR=f(xk+1)f(xk+1)-Fletcher-Reeves:-Polak-Ribierek+1=f(xk)f(xk)()()(1)(1)()1PRkkkkfffxfxfxPolakRibiere+=+-Hestenes-Stiefelk+1f(

55、xk)f(xk)HS=f(xk+1)f(xk+1)f(xk)45HestenesStiefel()()11kkkk+=pfx+fxConjugate gradients (CG)：Acomparisonoftheconvergenceofgradientdescent(ingreen)andconjugategradient(inred)Ingeneral,steepestdescentswillbringyouclosetothenearestlocalbringyouclosetothenearestlocalminimumveryquickly,whileconjugategradien

56、tsbringsyouveryconjugategradientsbringsyouveryclosetothelocalminimum,butperformsworsefarawayfromtheminimum.46Newtons methods:Newtons methodis a popular iterative method for finding the 0 of a one-dimensional function: f(x) = 0()()kkkgxx1=x()kxk+1xkgxx33x22x11x00L-BFGS is particularly well suited for

57、 optimization problems with a large number of dimensions.In practice this seems to converge faster than Conjugate Gradients, but due to the correction steps necessary it is not (yet) parallelized.三：Gromacs使用初步3.1Gromacs的安装33.22输入文件准备输入文件准备3.3作业运行、提交3.3作业运行、提交3.4练习3.1Gromacs的安装Install Gromacs4.5.5Che

58、ck compiler and mpienvironmentwhich ifortwhich iccwhichmpiccwhich mpicc.49INSTALLATIONInstall Gromacs4.5.5 Serial versionsource /opt/intel/Compiler/11.1/073/bin/ifortvars.sh intel64source /opt/intel/Compiler/11.1/073/bin/iccvars.sh intel64export MKL_PATH=/opt/intel/Compiler/11.1/073/mklexport CC=icc

59、./configure -prefix=/opt/gromacs-4.5.5makemake installmakke lik(til)links (optional)Note: default is to use “-enable-shared”50INSTALLATIONInstall Gromacs 4.5.5 Paralle version# before configgure, issue make distcleansource /opt/intel/Compiler/11.1/073/bin/ifortvars.sh intel64source /opt/intel/Compil

60、er/11.1/073/bin/iccvars.sh intel64export MKL_PATH=/opt/intel/Compiler/11.1/073/mklexport PATH=/opt/mpich2-1.4.1p1-intel-11.1-shared/bin:$PATHexport LD_LIBRARY_PATH=/opt/mpich2-1.4.1p1-intel-11.1-shared/lib:$LD_LIBRARY_PATHexport CC=mpiccexport LIBDIR=/opt/mpich2-1.4.1p1-intel-11.1-shared/libexport L