分子动力学模拟与药物研发
著：雅格布·D·德兰特，J·安德鲁·麦坎蒙，译：张浩晨
摘要：本文讨论了例如蛋白质等大分子结构的受体及其相应的小分子配体的多原子计算机模拟在药物研发中的作用。包括对于隐秘或存在变构现象的结合位点的识别，对于传统虚拟筛选方法的改进增强，以及对于小分子结合所需能量的直接预测。目前所用模拟方法的局限性，即较高的计算成本和所要求分子模拟力场的近似性，也在文中提及。随着计算机计算能力的提升和分子动力学算法的改进，计算机辅助药物设计的前景充满希望，其中分子动力学模拟技术将充当越来越重要的作用。
关键字：分子动力学模拟，计算机辅助药物设计，隐秘结合位点，变构结合位点，虚拟筛选，
自由能预测
介绍
1965年诺贝尔物理学奖获得者理查德·费曼，曾经有句名言说：“如果我们要做出一种假设来引导我们去探索生命的奥秘，那就是一切都是由原子构成的，而生命体所做的一切行为都可以理解为原子间有规律的或温和或剧烈的跃动。”过去的50年里，很多生物物理学一直致力于更好地了解原子的这种摇动和摆动的性质。微观世界的量子力学运动规律令那些只熟悉宏观动力学的人们感到惊讶。微观运动没有特定的规律，而是存在概率性；化学键的形成不是机械连接，而是通过改变以波和粒子状态存在的电子云。正如费曼所说的那样，大自然是那样的令人摸不着头脑。
理解这些荒谬而复杂但很神奇的分子运动规律无疑是与药物研发关系密切且很有帮助的。最初研究受体和配体结合的“锁和钥匙”理论认为容纳小分子配体的受体的结合域构象是如冰刻般一动不动，不存在任何的结构重排现象。现在这个观点已经基本被否定。取而代之的是一种新的观点，认为分子间的结合不仅引起构象的变化，同时也存在相结合分子间随机的摇摆运动。
软体动物的乙酰胆碱结合蛋白（AChBP）作为人类烟碱乙酰胆碱受体（AChR）的配体结合域结构和功能的替代，是无数可以用来说明原子运动重要性的例子之一（图1）。在结合于AChBP的小分子AChR激动剂的晶体结构显示，在结合处一个关键的环状结构部分闭锁（图 1 a,c）。相反，当体积较大的乙酰胆碱受体拮抗剂，如蛇的α -神经毒素结合到AchBP时的晶体结构表明，该环状结构扩大了至少10?，使得活性位点更加凸显（图1 b,c）。在研究未结合配体的样品的构象状态时显示，同时存在开放和闭合构象。伯恩等人提出未结合配体的AChBP或AChR是高度动态的蛋白质，推测蛋白构象在有选择地结合激动剂或拮抗剂后才趋于稳定。这些结合口袋构象中的任何一个都可能具有药效作用，因此也可能在药理学上有重要意义。如果这个配体结合的一般模型是正确的，则对于基于结构的药物设计的影响将是深远的。因为相似的结合原则可能也适用于许多其它体系。
分子动力学模拟
虽然这些晶体结构的研究有力地证明了蛋白质受体结合配体时其结构发生重排的特性起到重要的作用，但是得到这些结果所需要的费用和大量的劳动力，致使很多人转向去研究可以预测蛋白质运动的计算技术。不幸的是，模拟复杂的
量子力学运动和大分子体系的化学反应所需的计算往往过于复杂和密集，即使是用最快的超级计算机也很难胜任。分子动力学（MD）模拟，最早在20世纪70年代末开始发展，试图通过基于牛顿物理学的，对原子运动的简单近似模拟，来克服这个限制，从而降低计算的复杂度。用于近似的一般过程如图2所示。首先，利用核磁共振（NMR），X-射线晶体衍射或同源性建模等物理和计算机技术得到所研究的分子体系的计算机模型。作用于体系统各个原子上的力，是用类似图3所示的近似方程估计的。大体上讲原子的受力来自与之由化学键结合的原子和未与之结合的原子两个方面。化学键和键角使用简单的虚拟弹簧模型来模拟，二面角（即一个键的旋转）使用正弦函数模型近似模拟之遮掩键位形和交错键位形间的能量差。没有相互连接的原子间的分子作用依赖于两种力，以兰纳德-乔恩斯电势能为模型的范德华力，和利用库仑定律建模的静电荷相互作用。了解这些力学规律将更有利于我们去理解近似方程是如何通过参数描述原子间相互作用的。
为了能够真实地再现运动分子的实际情况，近似方程的能量部分是参数化的，以适应量子力学的计算和实验（例如，光谱）数据。被参数化的项包括：定义用来描述化学键和键角的弹簧的理想强度和长度，确定用于计算静电相互作用的最佳局部原子电荷，确定合适的范德华原子半径等。总的来说，这些参数被称作为“力场”，因为其描述了支配分子动力学的各种原子间相互作用力的分配。有几个力场在分子动力学模拟中经常使用，包括AMBER，CHARMM和GROMOS。尽管这些力场的参数化方式不同，但是基本上会得到相似的模拟结果。
一旦作用于体系中各原子上的力被计算出来，这些原子的位置则可以根据牛顿运动定律进行移动。随后仿真模拟向前推进，大约只有1到2个飞秒。通常该过程将被重复数百万次。由于大量的计算要求，分子动力学模拟的进行，通常使用由几十个甚至几百个处理器并行组成的计算机集群或超级计算机。许多流行的仿真模拟软件包，往往使用和其所用力场相同的名称（例如AMBER，CHARMM和NAMD），且都兼容信息传递接口（Message Passing Interface，MPI），即一个用于传输从一台计算机到另一台计算机之间信息的系统，使得在一个软件应用程序上同时运行多个处理器的系统来更大程度上方便复杂任务的执行。
然而分子动力学模拟的效用不应被无限制地夸大，一些研究者已经正在把这些模拟结果与实验数据作比较，以验证计算技术方法的可靠程度。 NMR数据是特别有用的，如通过分子动力学模拟获得的许多受体和配体的构象可以被用来预测NMR所测量构像的解自旋程度，可实现实验数据和理论结果之间的直接比较。实际上，一些研究已经表明大分子动力学的计算结果和实验测量值之间具有良好的一致性。
分子动力学模拟：目前的局限性
除这些成功的例子之外，分子动力学模拟的应用仍然受限于两个主要的挑战和局限：第一，所使用的力场仍需要进一步细化；第二，对计算的高要求使得常规的模拟长度被局限在一微秒以内，这导致在许多情况下，不能够采集足够的构象状态模拟数据。下面是一个高计算能力需求的例子：对一个相对较小的分子体系（约25,000原子）进行一次微秒级的分子动力学模拟，如在24处理器的机器上运行则需要几个月的时间才能完成。
分子动力学模拟的除了高计算需求的挑战和局限外，所使用的力场也只是对原子间量子力学事实的近似模拟。虽然已经可以准确地预测许多主要的分子运动，但是使这些模拟适合于体系中的量子效应仍然是重要的。例如，当模拟体系包含
大量的金属原子时预测结果可能会不理想。
为了克服这一局限，一些研究人员向经典的分子动力学力场中引入了量子力学计算。酶的活性位点及其它有趣但在模拟中受限的领域根据量子力学定律进行模拟，而较大分子的运动体系则采用分子动力学进行模拟。虽然离使用量子力学来计算和描述整个体系的“理想”目标还有一定困难，但这种混合模拟技术已经被成功地用于研究多个分子体系。
除了化学键的断裂与形成和电子极化，体系中由化学键连接的原子间的电子流从一个原子核到另一个，是除了少数特例外通常被忽略的另一个量子力学效应。在经典的分子动力学模拟开始之前，每个原子都被分配了一个固定数值的一部分电荷。然而在现实中，围绕在原子周围的电子云总在不断根据它们的环境变化而改变，所以原子周围的电荷最好能够进行动态的变化和响应。尽管电子极化作用在模拟中的重要性已达成广泛的共识，但是被普遍接受的极化力场经过30多年的发展一直没有真正地实现，可以使用这些力场进行分子动力学模拟的体系仍是很少的。尽管如此，一些电极化分子动力学力场目前仍然正在开发之中，其实现会给未来的分子动力学模拟带来更高的精准度。
除了忽略量子力学效应，分子动力学研究也受限于传统模拟方法可以进行的时间尺度。要重现热力学性质和/或充分阐明与药物设计相关的结合位点配型，必须通过仿真模拟对该蛋白的所有可能构象状态进行探讨。不幸的是，许多生化过程，包括与药物结合有关的受体构象变化，发生的时间尺度比可以进行模拟的时间要长。一些重要的例子显示，分子动力学模拟目前只能够运行百万分之一秒，而事实上，大多数模拟在十亿分之一秒钟即可完成。
分子动力学模拟和药物研发
除目前所使用的力场和构象采样的局限之外，分子动力学模拟及其在蛋白质运动中的成果，在药物发现的研究中扮演了重要角色。正如一张跑步者的照片无法描述关于她的全程步幅，一个单一的蛋白质构象也无法完全描述蛋白质分子动力学。虽然通过核磁共振，X射线晶体衍射和同源建模生成的静态模型可以提供关于大分子结构的有价值信息，但分子识别和药物结合的过程是异常动态的。在溶液中，当像药物一样的小分子（例如，一个配体）接近其靶标位点（例如，受体），所接触的构象并不是单一的冻结结构，而是处于不断运动状态的大分子。
在某些少数情况下，蛋白质运动是有限的，并且所述配体可装配到一个相当静态稳定的构象中，就像一把钥匙配合到一个锁。更典型地，配体可能稳定结合于由其动态受体抽样得到的许多构象的集合，从中挑选可能是最适合结合的构象。结合后，配体可以进一步诱导未结合配体时通常不存在受体位的结构发生构象变化。总之，受体的运动显然在大多数小分子药物的结合中起到了至关重要的作用。分子动力学模拟已经可以提供几种技术来利用这些分子运动。
分子动力学模拟中自由能计算方法的发展
虽然分子对接程序更多的在运算速度而不是精度上取得最优化，更准确地讲，尽管计算是密集的，但用于预测结合亲和程度的技术是存在的。这些技术，包括热力学整合，单步程微扰和自由能微扰，在很大程度上是基于分子动力学模拟的。
自由能是化学反应状态的表征，这意味着在一个给定的例如药物与其受体结合的事件中，自由能变化仅由此事件发生之前的能量和发生后的能量来确定。虽然从初始到最终状态的变化路径可能影响受体与配体结合的动力学模拟，但其对
自由能却没有影响。存在多种反应路径：可能是配体是慢慢地向活性位点移动，而轻巧地与结合中心融合。可能是受体蛋白质完全展开其构象，然后再将配体分子折叠包裹。也可能是溶液中的配体被定向到一个特定的轨迹上运动，当移动到一个活性中心时停留若干秒钟。反应的过程步骤并不重要，自由能变化只与化学反应溶液的初始能量和配体受体结合后的能量差有关。
除一些明显的例子外，利用受体-配体系统足够长时间的动力学模拟来捕捉整个结合事件目前是不可行的。然而，使用于1984年提出的是一种称为“推导转化”的处理技术来计算药物的结合亲和力是可能的。这种转变思想与前面提到的轨迹模型没有很大的差别。在分子动力学模拟的过程中，通过配位原子产生的静电和范德华力足以避免不正确的人为因素。最终，配体不再与蛋白质受体或溶剂相互作用。对于所有实际目的，该配体已经消失了。转变过程是否全部真实存在并没有关系。自由能只表示状态函数，所以从初始到最终的状态变化路径，无论是真实的还是虚构的，都是无关紧要的。
在什么情况下虚拟的配体应以这种途径消除，目前尚不清楚。以何种方式结合的配体在分子动力学模拟运行中应被消除？对存在于溶液中的配体情况又是怎样的？为了解决这些问题，提出了如图4所示的基于热力学的化学转换处理技术。由于自由能是状态特征且自由能变化仅由初始能量和最终能量的差值来确定，一个体系从一个状态开始，沿着这个自由能的周期循环进行变换，当返回到同一个初始状态时总的自由能应该没有发生变化（也就是说，ΔGbind + ΔGprotein - ΔG - ΔGwater = 0）。我们注意到，ΔG在这个方程式中本身是零，根据图4下半部分中所示的两个构象状态变化显示，这意味着配体分子不与水溶液和任何构象的蛋白质受体发生相互作用。因此，ΔGbind + ΔGprotein - ΔGwater = 0，并且由此可得，ΔGbind = ΔGwater - ΔGprotein。这些方程表明，通过运行两个模拟，即一个模拟中与受体结合的配体没有作用，另一个模拟中溶液中的配体失效，有可能实现对药物结合能的估算，进而间接测量药物的效力。
类似的，当我们希望确定一个给定的化学变化是否会改善候选配体的效能，而对药物进行优化时，计算相对配体结合能是非常有用的。在这种情况下，利用给定配体的一个部分进行逐渐转变，明显优于直接研究整个全部结合过程。
结论
在这篇综述中，我们讨论了分子动力学模拟在药物研发中的许多方面的作用。由于配体结合和与其相关的重要大分子的运动是发生在百万分之一秒内的微观事件，要使用目前的实验技术对原子结合的能量学和力学有一个完整的了解是不可能实现的。分子动力学模拟在解释实验方法所不能达到的细节方面具有重要作用。
随着计算机运算能力和模拟算法设计的不断改进，计算机辅助药物设计的未来是充满有希望的；分子动力学模拟有可能在新型的临床治疗药物研究的发展中扮演着越来越重要的作用。
图1Lymnaea stagnalis (一种生物) 不同构象的乙酰胆碱结合蛋白。为了方便可视化，该蛋白质的一些部分已经被移除。 （a）具有烟碱样配体接合的封闭构象结合蛋白（PDB ID：1UW6），以蓝色显示。 （b）具有相同构象的烟碱样配体叠加的开放构象结合蛋白（PDB ID：1YI5），如图粉红色。 （c）以不同颜色的条带表示两种构象蛋白的差异。
图2 分子动力学模拟进行的流程示意图。首先，准备好受体-配体体系的计算机模型。即类似于图3中所示的用于估算作用在体系中每个原子上的分子力的方程式。各个原子位置的移动是根据牛顿运动定律计算决定的。随着模拟时间的推进，体系中原子的受力运动过程将被重复许多次。这个流程图是由凯·诺德兰特首先创建的版本改编而来的。
图 3 举例说明用方程式近似模拟支配分子运动的原子力。支配分子运动的原子力可被分为两类：一类是那些由化学键连接的分子间的相互作用，另一类是没有化学键连接的原子之间的相互作用。化学键和键角使用简单的弹簧模型，而二面角（即一个键的旋转）使用正弦函数模型，近似模拟遮掩键位形和交错键位形之间的能量差。没有相互连接的原子间的分子作用依赖于两种力，以兰纳德-乔恩斯电势能为模型的范德华力，和利用库仑定律建模的静电荷相互作用。
图 4 用于推导生物结合能量变化的热力循环。通常情况下，人们希望计算出分子间结合作用所产生的自由能变化，即ΔGbind，在图4顶部示出。然而，这基本上是不可实现的，因为分子动力学模拟能够运行的时间还没有长到足以展示整个结合事件。不过，我们可以利用分子动力学模拟进行一系列的能量变化转换，以推导出想要得到的结论。如图所示，ΔGprotein表示当与受体结合的配体作用无效的情况下产生的自由能变化。ΔG表示当一个不能与受体结合的，幽灵似的配体发生作用时产生的自由能变化，显然，由于幽灵配位体不能够与任何溶剂或受体原子产生相互作用，这种能量变化始终为零。而ΔGwater表示当在溶液中的未与受体结合的配位体作用无效而解离时产生的自由能变化。一个体系，沿着这个自由能的周期循环进行变换，则当返回到同一个初始状态时总的自由能没
有发生变化，因此，ΔGbind = ΔGwater - ΔGprotein。
Molecular Dynamics Simulations and Drug Discovery
Jacob D Durrant and J Andrew McCammon
Abstract This review discusses the many roles atomistic computer simulations of macromolecular (for example, protein) receptors and their associated small-molecule ligands can play in drug discovery, including the identification of cryptic or allosteric binding sites, the enhancement of traditional virtual-screening methodologies, and the direct prediction of small-molecule binding energies. The limitations of current simulation methodologies, including the high computational costs and approximations of molecular forces required, are also discussed. With constant improvements in both computer power and algorithm design, the future of computer-aided drug design is promising; molecular dynamics simulations are likely to play an increasingly important role.
Keywords molecular dynamics simulations, computer-aided drug discovery, cryptic binding sites, allosteric binding sites, virtual screening, free-energy prediction.
Introduction
Richard Feynman, recipient of the 1965 Nobel Prize in Physics, once famously stated: ‘If we were to name the most powerful assumption of all, which leads one on and on in an attempt to understand life, it is that all things are made of atoms, and that everything that living things do can be understood in terms of the jigglings and wigglings of atoms.’ Much of the biophysics of the last 50 years has been dedicated to better understanding the nature of this atomic jiggling and wiggling. The quantum-mechanical laws governing motions in the microscopic world are surprisingly foreign to those familiar with macroscopic dynamics. Motions are governed not by deterministic laws, but by probability functions; chemical bonds are formed not mechanically, but by shifting clouds of electrons that are simultaneously waves and particles. As Feynman eloquently put it, this is ‘nature as she is - absurd’.
Understanding these absurd molecular motions is undoubtedly germane to drug discovery. The initial ‘lock-and-key’ theory of ligand binding, in which a frozen, motionless receptor was thought to accommodate a small molecule without undergoing any conformational rearrangements, has been largely abandoned in favor of binding models that account not only for conformational changes, but also for the random jiggling of receptors and ligands.
The mollusk acetylcholine binding protein (AChBP), a structural and functional surrogate of the human nicotinic acetylcholine receptor (AChR) ligand-binding domain, is one of countless examples that illustrate the importance of accounting for these atomistic motions (Figure 1). In crystallographic structures of small-molecule AChR agonists bound to AChBP, a key loop (loop C) partially closes around the ligand (Figure 1a,c). In contrast, crystal structures of large AChR antagonists like snake α-neurotoxins bound to AChBP reveal that this same loop is displaced by as much as 10 ?, producing an active site that is far more open (Figure 1b,c). Bourne et al. proposed that the unbound AChBP and AChR are highly dynamic proteins that, in the absence of a ligand, sample many conformational states, both open and closed, that are selectively stabilized by the binding of agonists and antagonists. Any one of these binding-pocket conformations may be druggable and therefore pharmacologically relevant. If this general model of ligand binding is correct, the implications for structure-based drug design are profound, as the same principle of binding likely applies to many other systems as well.
Molecular dynamics simulations
While crystallographic studies like these convincingly demonstrate the important role protein flexibility plays in ligand binding, the expense and extensive labor required to generate them have led many to seek computational techniques that can predict protein motions. Unfortunately, the calculations required to describe the absurd quantum-mechanical motions and chemical reactions of large molecular systems are often too complex and computationally intensive for even the best supercomputers. Molecular dynamics (MD) simulations, first developed in the late 1970s, seek to overcome this limitation by using simple approximations based on Newtonian physics to simulate atomic motions, thus reducing the computational complexity. The general process of approximation used is outlined in Figure 2. First, a computer model of the molecular system is prepared from nuclear magnetic resonance (NMR), crystallographic, or homology-modeling data. The forces acting on each of the system atoms are then estimated from an equation like that shown in Figure 3. In brief, forces arising from interactions between bonded and non-bonded atoms contribute. Chemical bonds and atomic angles are modeled using simple virtual springs, and dihedral angles (that is, rotations about a bond) are modeled using a sinusoidal function that approximates the energy differences between eclipsed and staggered conformations. Non-bonded forces arise due to van der Waals interactions, modeled using the Lennard-Jones 6-12 potential, and charged (electrostatic) interactions, modeled using Coulomb’s law. For a more in-depth review describing how the equations describing these interactions are parameterized, see.
In order to reproduce the actual behavior of real molecules in motion, the energy terms described above are parameterized to fit quantum-mechanical calculations and experimental (for example, spectroscopic) data. This parameterization includes identifying the ideal stiffness and lengths of the springs that describe chemical bonding and atomic angles, determining the best partial atomic charges used for calculating electrostatic-interaction energies, identifying the proper van der Waals atomic radii, and so on. Collectively, these parameters are called a ‘force field’ because they describe the contributions of the various atomic forces that govern molecular dynamics. Several force fields are commonly used in molecular dynamics simulations, including AMBER, CHARMM, and GROMOS. These differ principally in the way they are parameterized but generally give similar results.
Once the forces acting on each of the system atoms have been calculated, the positions of these atoms are moved according to Newton’s laws of motion. The simulation time is then advanced, often by only 1 or 2 quadrillionths of a second, and the process is repeated, typically millions of times. Because so many calculations are required, molecular dynamics simulations are typically performed on computer clusters or supercomputers using dozens if not hundreds of processors in parallel. Many of the most popular simulation software packages, which often bear the same names as their default force fields (for example AMBER, CHARMM, and NAMD), are compatible with the Message Passing Interface (MPI), a system of computer-to-computer messaging that greatly facilitates the execution of complex tasks by one software application on multiple processors operating simultaneously.
While the utility of molecular dynamics simulations should not be overstated, a number of studies comparing these simulations with experimental data have been used to validate the computational technique. NMR data are particularly useful, as the many receptor and ligand conformations sampled by molecular dynamics simulations can be used to predict NMR
measurements like spin relaxation, permitting direct comparison between experimental and theoretical techniques. Indeed, a number of studies have shown good agreement between computational and experimental measurements of macromolecular dynamics.
Molecular dynamics simulations: current limitations
These successes aside, the utility of molecular dynamics simulations is still limited by two principal challenges: the force fields used require further refinement, and high computational demands prohibit routine simulations greater than a microsecond in length, leading in many cases to an inadequate sampling of conformational states. As an example of these high computational demands, consider that a one-microsecond simulation of a relatively small system (approximately 25,000 atoms) running on 24 processors takes several months to complete.
Aside from challenges related to the high computational demands of these simulations, the force fields used are also approximations of the quantum-mechanical reality that reigns in the atomic regime. While simulations can accurately predict many important molecular motions, these simulations are poorly suited to systems where quantum effects are important, for example, when transition metal atoms are involved in binding.
To overcome this challenge, some researchers have introduced quantum mechanical calculations into classic molecular-dynamics force fields; the motions and reactions of enzymatic active sites or other limited areas of interest are simulated according to the laws of quantum mechanics, and the motions of the larger system are approximated using molecular dynamics. While far from the computationally intractable ‘ideal’ of using quantum mechanics to describe the entire system, this hybrid technique has nevertheless been used successfully to study a number of systems.
Aside from bond breaking and formation, electronic polarization, caused by the flow of electrons from one atomic nucleus to another among groups of atoms that are chemically bonded, is another quantum-mechanical effect that, with few exceptions, is generally ignored. In classical molecular dynamics simulations, each atom is assigned a fixed partial charge before the simulation begins. In reality, however, the electron clouds surrounding atoms are constantly shifting according to their environments, so that the partial charges of atoms would be better represented as dynamic and responsive. Despite wide agreement on the importance of accounting for electronic polarization, after 30 years of development a generally accepted polarizable force field has not been forthcoming, and molecular dynamics simulations using those force fields that are available are rare. Nevertheless, a number of polarizable force fields are currently under development, and future implementations may lead to improved accuracy.
In addition to neglecting quantum-mechanical effects, molecular dynamics studies are also limited by the short time scales typically simulated. To reproduce thermodynamic properties and/or to fully elucidate all binding-pocket configurations relevant to drug design, all the possible conformational states of the protein must be explored by the simulation. Unfortunately, many biochemical processes, including receptor conformational shifts relevant to drug binding, occur on time scales that are much longer than those amenable to simulation. With some important exceptions, simulations are currently limited to at most millionths of a second; indeed, most simulations are measured in billionths of a second.
Molecular dynamics simulations and drug discovery
Weaknesses in current force fields and conformational sampling aside, molecular dynamics simulations and the insights they offer into protein motion often play important roles in drug discovery. Just as a single photograph of a runner tells little about her stride, a single protein conformation tells little about protein dynamics. The static models produced by NMR, X-ray crystallography, and homology modeling provide valuable insights into macromolecular structure, but molecular recognition and drug binding are very dynamic processes. When a small molecule like a drug (for example, a ligand) approaches its target (for example, a receptor) in solution, it encounters not a single, frozen structure, but rather a macromolecule in constant motion.
In some rare cases, protein motions are limited, and the ligand may fit into a fairly static binding pocket like a key fits into a lock. More typically, the ligand may bind and stabilize only a subset of the many conformations sampled by its dynamic receptor, causing the population of all possible conformations to shift towards those that can best accommodate binding. Upon binding, the ligand may further induce conformational changes that are not typically sampled when the ligand is absent. Regardless, receptor motions clearly play an essential role in the binding of most small-molecule drugs. Several techniques have been developed to exploit the information about these motions that molecular dynamics simulations can provide.
Advanced free-energy calculations using molecular dynamics simulations Though docking programs are optimized for speed rather than accuracy, more accurate, albeit computationally intensive, techniques for predicting binding affinity do exist. These techniques, which include thermodynamic integration, single-step perturbation, and free energy perturbation, are based in large part on molecular dynamics simulations.
Free energy is a state function, meaning that the free-energy difference associated with a given event like a drug binding to its receptor is determined only by the energy prior to that event and the energy following it; while the path taken from the initial to the final state may influence receptor-ligand kinetics, it has no bearing on the free energy. Perhaps the ligand diffuses slowly towards the active site and slips easily into the binding pocket. Perhaps the protein unfolds entirely and then refolds around the ligand. Perhaps the ligand in solution is beamed to a starship in orbit, only to rematerialize in the active site a few seconds later. The mechanics do not matter; the free energy depends only on the initial energy in solution and the final energy following the binding event.
With some notable exceptions, simulating a receptor-ligand system long enough to capture an entire binding event is not currently feasible. However, it is still possible to calculate a drug’s binding affinity using a technique called ‘alchemical transformation’, first described in 1984. This transformation is not too different from the starship example given above. During the course of a molecular dynamics simulation, the electrostatic and van der Waals forces produced by ligand atoms are turned down gradually enough to avoid undesirable artifacts. Eventually, the ligand is no longer able to interact with the protein or solvent. For all practical purposes, the ligand has disappeared. It does not matter that this transformation is not at all physical; free energy is a state function, so the path from the initial to the final state, whether real or imaginary, is irrelevant.
It is not clear, however, in what context the ligand should be figuratively annihilated in this way. Should a molecular dynamics simulation be run in which the bound ligand vanishes? What about the ligand in solution? To address these questions, alchemical transformations are selected based on the thermodynamic cycle shown in Figure 4. As free energy is a state function that
depends only on the energy of the initial and final states, a system that proceeds from one state around this free-energy cycle only to return to the same initial state should have no change in total free energy (that is, ΔGbind + ΔGprotein - ΔG - ΔGwater = 0). We note that ΔG in this equation is itself zero, since the ligand is entirely disappeared in both of the states shown in the bottom half of Figure 4, meaning the ligand cannot interact with the water solvent or the protein in either state. Thus, ΔGbind + ΔGprotein - ΔGwater = 0, and, consequently, ΔGbind = ΔGwater - ΔGprotein. These equations demonstrate that it is possible to estimate a drug’s free energy of binding, an indirect measurement of drug potency, by running two simulations, one in which the receptor-bound ligand disappears, and one in which the solvated ligand disappears.
A similar task, calculating relative ligand binding energies, is useful during drug optimization when one wishes to determine if a given chemical change will improve the potency of a candidate ligand. In this case, rather than annihilating the entire ligand, one section of the ligand is gradually transformed.
Conclusions
In this review, we have discussed the many roles that molecular dynamics simulations can play in drug discovery. As ligand binding and the important macromolecular motions associated with it are microscopic events that take place in mere millionths of a second, a complete understanding of the atomistic energetics and mechanics of binding is unattainable using current experimental techniques. Molecular dynamics simulations are useful for filling in the details where experimental methods cannot.
With constant improvements in both computer power and algorithm design, the future of computer-aided drug design is promising; molecular dynamics simulations are likely to play an increasingly important role in the development of novel pharmacological therapeutics.
Figure 1. The different conformations of the acetylcholine binding protein from Lymnaea stagnalis. Portions of the protein have been removed to facilitate visualization. (a) The protein in a closed conformation with nicotine bound (PDB ID: 1UW6), shown in blue. (b) The protein in an open conformation (PDB ID: 1YI5) with the same nicotine conformation superimposed on the structure, shown in pink. (c) Ribbon representations of the two conformations.
Figure 2. A schematic showing how a molecular dynamics simulation is performed. First, a computer model of the receptor-ligand system is prepared. An equation like that shown in Figure 3 is used to estimate the forces acting on each of the system atoms. The positions of the atoms are moved according to Newton’s laws of motion. The simulation time is advanced, and the process is repeated many times. This figure was adapted from a version originally created by Kai Nordlund.
Figure 3. An example of an equation used to approximate the atomic forces that govern molecular movement. The atomic forces that govern molecular movement can be divided into those caused by interactions between atoms that are chemically bonded to one another and those caused by interactions between atoms that are not bonded. Chemical bonds and atomic angles are modeled using simple springs, and dihedral angles (that is, rotations about a bond) are modeled using a sinusoidal function that approximates the energy differences between eclipsed and staggered conformations. Non-bonded forces arise due to van der Waals interactions, modeled using the Lennard-Jones potential, and charged (electrostatic) interactions, modeled using
Coulomb’s law.
Figure 4. The thermodynamic cycle used for selecting alchemical transformations. Typically, one wishes to calculate the free energy of binding, ΔG bind, shown across the top. However, it is generally impractical to run a molecular dynamics simulation long enough to capture an entire binding event. Instead, a series of alchemical transformations are performed using molecular dynamics simulations. ΔG protein is the change in free energy that occurs when a bound ligand is ‘annihilated’. ΔG is the change in free energy that occurs when an unbound ‘ghost’ ligand binds to the receptor; however, since a ghost ligand is not able to interact with any solvent or receptor atoms, this energy is always zero. Finally, ΔG water is the change in free energy that occurs when an unbound ligand in solution is ‘annihilated’. A system that proceeds from one state around this free-energy cycle only to return to the same initial state should have no change in total free energy; consequently, ΔG bind = ΔG water - ΔG protein.
分子动力学的模拟过程 分子动力学模拟作为一种应用广泛的模拟计算方法有其自身特定的模拟步骤,程序流程也相对固定。本节主要就分子动力学的模拟步骤和计算程序流程做一些简单介绍。 1. 分子动力学模拟步驟 分子动力学模拟是一种在微观尺度上进行的数值模拟方法。这种方法既可以得到一些使用传统方法,热力学分析法等无法获得的微观信息,又能够将实际实验研究中遇到的不利影响因素回避掉,从而达到实验研宄难以实现的控制条件。 分子动力学模拟的步骤为: (1)选取所要研究的系统并建立适当的模拟模型。 (2)设定模拟区域的边界条件,选取粒子间作用势模型。 (3)设定系统所有粒子的初始位置和初始速度。 (4)计算粒子间的相互作用力和势能,以及各个粒子的位置和速度。 (5)待体系达到平衡,统计获得体系的宏观特性。 分子动力学模拟的主要对象就是将实际物理模型抽象后的物理系统模型。因此,物理建模也是分子动力学模拟的一个重要的环节。而对于分子动力学模拟,主要还是势函数的选取,势函数是分子动力学模拟计算的核心。这是因为分子动力学模拟主要是计算分子间作用力,计算粒子的势能、位置及速度都离不开势函数的作用。系统中粒子初始位置的设定最好与实际模拟模型相符,这样可以使系统尽快达到平衡。另外,粒子的初始速度也最好与实际系统中分子的速度相当,这样可以减少计算机的模拟时间。 要想求解粒子的运动状态就必须把运动方程离散化,离散化的方法有经典Verlet算法、蛙跳算法(Leap-frog)、速度Veriet算法、Gear预估-校正法等。这些算法有其各自的优势,选取时可按照计算要求选择最合适的算法。 统计系统各物理量时,便又涉及到系统是选取了什么系综。只有知道了模拟系统采用的系综才能釆用相对应的统计方法更加准确,有效地进行统计计算,减少信息损失。 2. 分子动力学模拟程序流程 具体到分子动力学模拟程序的具体流程,主要包括: (1)设定和模拟相关的参数。 (2)模拟体系初始化。 (3)计算粒子间的作用力。 (4)求解运动方程。 (5)循环计算,待稳定后输出结果。 分子动力学模拟程序流程图如2.3所示。
分子动力学方法模拟基本步骤 1.第一步 即模型的设定，也就是势函数的选取。势函数的研究和物理系统上对物质的描述研究息息相关。最早是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是LJ势函数，还有EAM势函数，不同的物质状态描述用不同的势函数。 模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。 2 第二步 给定初始条件。运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如：verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。 一般意思上讲系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为模拟实践足够长时，系统就会忘掉初始条件。当然，合理的初始条件可以加快系统趋于平衡的时间和步伐，获得好的精度。 常用的初始条件可以选择为：令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。 第三步 趋于平衡计算。在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统的能量，并且这个状态本身也还不是一个平衡态。 为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，我们增加或者从系统中移出能量，直到持续给出确定的能量值。我们称这时的系统已经达到平衡。这段达到平衡的时间成为驰豫时间。 分子动力学中，时间步长的大小选择十分重要，决定了模拟所需要的时间。为了减小误差，步长要小，但小了系统模拟的驰豫时间就长了。因此根据经验选择适当的步长。如，对一个具有几百个氩气Ar分子的体系，lj势函数，发现取h为0.01量级，可以得到很好的相图。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s的量级呢。如果模拟1000步，系统达到平衡，驰豫时间只有10-11s。 第四步 宏观物理量的计算。实际计算宏观的物理量往往是在模拟的最后揭短进行的。它是沿相空间轨迹求平均来计算得到的（时间平均代替系综平均）
《装备制造技术》 2007年第 10期 收稿日期 :2007-08-21 作者简介 :申海兰 , 24岁 , 女 , 河北人 , 在读研究生 , 研究方向为微机电系统。 分子动力学模拟方法概述 申海兰 , 赵靖松 (西安电子科技大学机电工程学院 , 陕西西安 710071 摘要 :介绍了分子动力学模拟的基本原理及常用的原子间相互作用势 , 如Lennard-Jones 势 ; 论述了几种常用的有限差分算法 , 如 Verlet 算法 ; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。关键词 :分子动力学模拟 ; 原子间相互作用势 ; 有限差分算法 ; 系综中图分类号 :O3 文献标识码 :A 文章编号 :1672-545X(200710-0029-02 从统计物理学中衍生出来的分子动力学模拟方法 (molec- ular dynamics simulation , M DS , 实践证明是一种描述纳米科技 研究对象的有效方法 , 得到越来越广泛的重视。所谓分子动力学模拟 , 是指对于原子核和电子所构成的多体系统 , 用计算机模拟原子核的运动过程 , 从而计算系统的结构和性质 , 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动 [1]。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段 , 称之为“计算机实验” 手段 [2], 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。
根据模拟对象的不同 , 将它分为平衡态分子动力学模拟 (EM DS (和非平衡态分子动力学模拟 (NEM DS 。其中 , EM DS 是分子动力学模拟的基础 ; NEM DS 适用于非线性响应系统的模拟 [3]。下面主要介绍 EM DS 。 1分子动力学方法的基本原理 计算中根据以下基本假设 [4]: (1 所有粒子的运动都遵循经典牛顿力学规律。 (2 粒子之间的相互作用满足叠加原理。 显然这两条忽略了量子效应和多体作用 , 与真实物理系统存在一定差别 , 仍然属于近似计算。 假设 N 为模拟系统的原子数 , 第 i 个原子的质量为 m i , 位置坐标向量为 r i , 速度为 v i =r ? i , 加速度为 a i =r ?? i , 受到的作用力为 F i , 原子 i 与原子 j 之间距离为 r ij =r i -r j , 原子 j 对原子 i 的作用力为 f ij , 原子 i 和原子 j 相互作用势能为 ! (r ij , 系统总的势能为 U (r 1, r 2, K r N = N i =1! j ≠ i ! " (r ij , 所有的物理量都是随时 间变化的 , 即 A=A (t , 控制方程如下 : m i r ?? i =F i =j ≠ i
分子动力学模拟 分子动力学就是一门结合物理,数学与化学的综合技术。分子动力学就是一套分子模拟方法,该方法主要就是依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量与其她宏观性质。 这门技术的发展进程就是: 1980年:恒压条件下的动力学方法(Andersenの方法、Parrinello-Rahman法) 1983年:非平衡态动力学方法(Gillan and Dixon) 1984年:恒温条件下的动力学方法(能势‐フーバーの方法) 1985年:第一原理分子动力学法(→カー?パリネロ法) 1991年:巨正则系综的分子动力学方法(Cagin and Pettit)、 最新的巨正则系综,即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触,此时系统不仅同热源有能量交换,而且可以同粒子源有粒子的交换,最后达到平衡,这种系综称巨正则系综。 进行分子动力学模拟的第一步就是确定起始构型,一个能量较低的起始构型就是进行分子模拟的基础,一般分子的其实构型主要就是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度,这一速度就是根据玻尔兹曼分布随机生成,由于速度的分布符合玻尔兹曼统计,因此在这个阶段,体系的温度就是恒定的。另外,在随机生成各个原子的运动速度之后须进行调整,使得体系总体在各个方向上的动量之与为零,即保证体系没有平动位移。 由上一步确定的分子组建平衡相,在构建平衡相的时候会对构型、温度等参数加以监控。 进入生产相之后体系中的分子与分子中的原子开始根据初始速度运动,可以想象其间会发生吸引、排斥乃至碰撞,这时就根据牛顿力学与预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算,在这个过程中,体系总能量不变,但分子内部势能与动能不断相互转化,从而体系的温度也不断变化,在整个过程中,体系会遍历势能面上的各个点,计算的样本正就是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能,进而计算构型积分。 作用势的选择与动力学计算的关系极为密切,选择不同的作用势,体系的势能面会有不同的形状,动力学计算所得的分子运动与分子内部运动的轨迹也会不同,进而影响到抽样的结果与抽样结果的势能计算,在计算宏观体积与微观成分关系的时候主要采用刚球模型的二体势,计算系统能量,熵等关系时早期多采用Lennard-Jones、morse势等双体势模型,对于金属计算,主要采用morse势,但就是由于通过实验拟合的对势容易导致柯西关系,与实验不符,因此在后来的模拟中有人提出采用EAM等多体势模型,或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但就是对于二体势模型,多体势往往缺乏明确的表达式,参量很多,模拟收敛速度很慢,给应用带来很大困难,因此在一般应用中,通过第一性原理计算结果拟合势函数的L-J,morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想就是赋予分子体系初始运动状态之后,利用分子的自然运动在相空间中抽取样本进行统计计算,时间步长就就是抽样的间隔,因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞,体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力,因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但就是通常情况下,体系各自由度中运动周期最短的就是各个化学键的振动,而这种运动对计算某些宏观性质并不产生影响,因此就产生了屏蔽分子内部振动或其她无关运动的约束动力学,约束动力学可以有效地增长分子动力学模拟时间步长,提高搜索相空间的能
第3章 铁基块体非晶合金‐纳米晶转变的动力学模拟过程 3.1 Discover模块 3.1.1 原子力场的分配 在使用Discover模块建立基于力场的计算中，涉及几个步骤。主要有：选择力场、指定原子类型、计算或指定电荷、选择non‐bond cutoffs。 在这些步骤中，指定原子类型和计算电荷一般是自动执行的。然而，在某些情形下需要手动指定原子类型。原子定型使用预定义的规则对结构中的每个原子指定原子类型。在为特定的系统确定能量和力时，定型原子使工作者能使用正确的力场参数。通常，原子定型由Discover使用定型引擎的基本规则来自动执行，所以不需要手动原子定型。然而，在特殊情形下，人们不得不手动的定型原子，以确保它们被正确地设置。 图 3-1 1）计算并显示原子类型：点击Edit→Atom Selection，如图3‐1所示 图3-2 弹出对话框，如图3‐2所示 从右边的…的元素周期表中选择Fe，再点Select，此时所建晶胞中所有Fe
原子都将被选中，原子被红色线圈住即表示原子被选中。再编辑集合，点击Edit →Edit Sets，如图3‐3、3‐4所示。 图3-3 图3-4 弹出对话框见图3‐4，点击New...，给原子集合设定一个名字。这里设置为Fe，则3D视图中会显示“Fe”字样，再分配力场：在工具栏上点击Discover按 钮，从下拉列表中选择Setup，显示Discover Setup对话框，选择Typing选项卡，见图3‐5。 图3-5 在Forcefield types里选择相应原子力场，再点Assign（分配）按钮进行原子力场分配。注意原子力场中的价态要与Properties Project里的原子价态（Formalcharge）一致。
分子动力学模拟 一,软件： NAMD:https://www.doczj.com/doc/eb17721210.html,/Research/namd/免费注册之后进行免费下载, 只需要下载解压不需要安装 VMD:https://www.doczj.com/doc/eb17721210.html,/Research/vmd/免费,分子可视化和辅助分析软 件 二，分子动力学模拟需要的数据文件包括： (1)蛋白质的PDB文件，此文件只记录原子空间位置，能够从RCSB管理的PDB数据库(https://www.doczj.com/doc/eb17721210.html,/pdb/)下载。 (2)PSF文件，此文件负责储存蛋白质的结构信息，记录蛋白质原子之间的成键情况。用户需要根据自己要求生成该文件。 (3)力场参数文件。此文件是分子动力学模拟的核心。CHAYMM，X-PLOR，AMBER和GROMACS 是经常用到的四种力场。NAMD能够利用上述每一种力场执行分子动力学模拟。 (4)配置文件(configuration file)。此文件作用是告知NAMD分子动力学模拟的各种参数，例如PDB和PSF两个文件保存的位置，模拟结果储存在哪里，体系的温度是多少等等。此文件也是要用户根据需求自己生成。同一配置的电脑，蛋白质分子大小不同，模拟运行的时间也不同，通常大蛋白质需要较长的时间。 三．以蛋白质1L63为例给出操作说明。 在PDB数据库下载蛋白质1L63. 建立文件夹1L63，其中包括以下几个文件，其中.conf文件需要修改,下面第4步会讲到。 以下生成PSF文件： 1.单击VMD，file-New Molecule-打开Molecule File Browser对话框，单击Browse按钮，在文件浏览器中找到文件夹1L63,在此文件夹中选择1L63.pdb,单击Load按钮载入1L63.pdb 2.除去pdb文件中带有的水分子 单击Extension-TK Console，弹出VMD Tk Console窗口。 首先用cd命令改变当前目录到1L63文件夹下，然后输入下列命令： set L63[atomselect top protein] $L63writepdb L63p.pdb 这样，1L63文件夹下就生成了文件L63P.pdb。这一PDB文件仅包含蛋白质，不包含水分子。 3.生成psf文件。 注意，这里仅讲全自动的psf文件生成器，描述如下： 选择Extensions-Modeling-Automatic PSF Builder菜单项，点击左上角的Options,选择Add solvation box，和Add neutralizing ions,点击右下角的I’m feeling lucky按钮，
分子动力学模拟
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法）1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。
进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动，可以想象其间会发生吸引、排斥乃至碰撞，这时就根据牛顿力学和预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算，在这个过程中，体系总能量不变，但分子内部势能和动能不断相互转化，从而体系的温度也不断变化，在整个过程中，体系会遍历势能面上的各个点，计算的样本正是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能，进而计算构型积分。 作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、morse势等双体势模型，对于金属计算，主要采用morse 势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的L-J，morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想是赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞，体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，因此
分子动力学模拟基础知识 ? Molecular Dynamics Simulation o MD: Theoretical Background Newtonian Mechanics and Numerical Integration The Liouville Operator Formalism to Generating MD Integration Schemes o Case Study 1: An MD Code for the Lennard-Jones Fluid Introduction The Code, mdlj.c o Case Study 2: Static Properties of the Lennard-Jones Fluid (Case Study 4 in F&S) o Case Study 3: Dynamical Properties: The Self-Diffusion Coefficient ? Ensembles o Molecular Dynamics at Constant Temperature Velocity Scaling: Isokinetics and the Berendsen Thermostat Stochastic NVT Thermostats: Andersen, Langevin, and Dissipative Particle Dynamics The Nosé-Hoover Chain Molecular Dynamics at Constant Pressure: The Berendsen Barostat Molecular Dynamics Simulation We saw that the Metropolis Monte Carlo simulation technique generates a sequence of states with appropriate probabilities for computing ensemble averages (Eq. 1). Generating states probabilitistically is not the only way to explore phase space. The idea behind the Molecular Dynamics (MD) technique is that we can observe our dynamical system explore phase space by solving all particle equations of motion . We treat the particles as classical objects that, at least at this stage of the course, obey Newtonian mechanics. Not only does this in principle provide us with a properly weighted sequence of states over which we can compute ensemble averages, it additionally gives us time-resolved information, something that Metropolis Monte Carlo cannot provide. The ``ensemble averages'' computed in traditional MD simulations are in practice time averages : (99) The ergodic hypothesis partially requires that the measurement time, , i , in the system. The price we pay for this extra information is that we must at least access if not store particle velocities in addition to positions, and we must compute interparticle forces in addition to potential energy. We will introduce and explore MD in this section.
Gromacs中文教程 淮海一粟 分子动力学（MD）模拟分为三步：首先，要准备好模拟系统；然后，对准备好的系统进行模拟；最后，对模拟结果进行分析。虽然第二步是最耗费计算资源的，有时候需要计算几个月，但是最耗费体力的步骤在于模拟系统准备和结果分析。本教程涉及模拟系统准备、模拟和结果分析。 一、数据格式处理 准备好模拟系统是MD最重要的步骤之一。MD模拟原子尺度的动力学过程，可用于理解实验现象、验证理论假说，或者为一个待验证的新假说提供基础。然而，对于上述各种情形，都需要根据实际情况对模拟过程进行设计；这意味着模拟的时候必须十分小心。 丢失的残基、原子和非标准基团 本教程模拟的是蛋白质。首先需要找到蛋白质序列并选择其起始结构，见前述；然后就要检查这个结构是否包含所有的残基和原子，这些残基和原子有时候也是模拟所必需的。本教程假定不存在缺失，故略去。 另一个需要注意的问题是结构文件中可能包含非标准残基，被修饰过的残基或者配体，这些基团还没有力场参数。如果有这些基团，要么被除去，要么就需要补充力场参数，这牵涉到MD的高级技巧。本教程假定所有的蛋白质不含这类残基。 结构质量 对结构文件进行检查以了解结构文件的质量是一个很好的练习。例如，晶体结构解析过程中，对于谷氨酰胺和天冬酰胺有可能产生不正确的构象；对于组氨酸的质子化状态和侧链构象的解析也可能有问题。为了得到正确的结构，可以利用一些程序和服务器(如 WHATIF)。本教程假定所用的结构没有问题，我们只进行数据格式处理。 二、结构转换和拓扑化 一个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。由于.pdb 结构文件只含有原子坐标，我们首先必须建立拓扑文件，该文件描述了原子类型、电荷、成键情况等信息。拓扑文件对应着一种力场，选择何种力场对于拓扑文件的建立是一个值得仔细考虑的问题。这里我们用的是GROMOS96 53a6连接原子力场，该力场对于氨基酸侧链的自由能预测较好，并且与NMR试验结果较吻合。
分子动力学模拟 The Standardization Office was revised on the afternoon of December 13, 2020
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法）1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法） 1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。 由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动，可以想象其间会发生吸引、排斥乃至碰撞，这时就根据牛顿力学和预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算，在这个过程中，体系总能量不变，但分子内部势能和动能不断相互转化，从而体系的温度也不断变化，在整个过程中，体系会遍历势能面上的各个点，计算的样本正是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能，进而计算构型积分。 作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、morse势等双体势模型，对于金属计算，主要采用morse势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的L-J，morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想是赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞，体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但是通常情况下，体系各自由度中运动周期最短的是各个化学键的振动，而这种运动对计算某些宏观性质并不产生影响，因此就产生了屏蔽分子内部振动或其他无关运动的约束动力学，约束动力学可以有效地增长分子动力学模拟时间步长，提高搜索相空间的能力。
《材料计算设计基础》 学号： 流水号： 姓名： 完成日期：
分子动力学模拟及其在材料中的研究进展 摘要：本文综述了分子动力学模拟技术的发展，介绍了分子动力学的分类、运动方程的求解、初始条件和边界条件的选取、平衡系综及其控制、感兴趣量的提取以及分子动力学模拟在材料中的研究进展。 关键词：分子动力学模拟平衡态系综金属材料感兴趣量径向分布函数 引言 科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在金属材料研究中显得非常有吸引力。 分子动力学MD (Molecular Dynamics)模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。MD模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法”, 它所计算的是时间平均,而MC进行的是系综平均。然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的比较计算能给出几乎相同的结果。 经典的分子动力学方法是Alder等于1957年提出并首先在“硬球”液体模型下应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变”。后来人们称这种相变为Alder相变。Rahman
物理计算与设计报告书 院（系）名称： 学生姓名： 专业名称： 班级： 时间： 金属铝分子动力学模拟
摘要：分子动力学模拟，是指对于原子核和电子所构成的多体系统，用计算机模拟 原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在全部其它 原子核和电子所提供的经验势场作用下按牛顿定律运动。我们用c语言编写程序，VMD 动画演示得到原子在拉伸过程中的变化。在控制温度不变的情况下，得到了金属铝分子 的动力学模拟过程。通过不断拉伸，趋衡铝分子，计算其势能，力,速度，观察每次拉 伸过程中以及拉伸后铝原子的排列，得到金属铝的运动细节，从而更加利于我们了解铝 的性质。 结论：原子两端的拉力与原子势能的变化曲线基本一致。原子间断层以滑层方式断 裂。 关键词：铝分子，分子动力学，c语言，势能 1 引言 人们很早就知道材料的力学性能随尺度发生变化尺度减小, 材料中缺陷存在的几率降低, 材料的强度提高同时尺度的变化可能导致材料内在变形竞争机制的改变, 例如多晶材料晶粒粒径在微米级以上时, 强度主要受位错强化机制控制, 而粒径进入纳米级后, 材料的变形主要来源于晶界滑移等机制原子尺度下, 微观效应占主导地位, 材料的理化、力学性能表现出与宏观不同、甚至相反的特性。Brenner发现金属单晶晶须拉伸强度与晶须直径呈反比,Fleck在微米级细铜丝的扭转试验中观察到尺寸效应纳米电机系统(NEMS)的出现同迫切要求了解纳米尺度下材料的力学行为, 当前从实验上较难获得详细的信息, 而分子动力学模拟可以提供相关细节. 分子动力学通过直接模拟原子的运动过程, 使我们能够详细了解模拟对象的演化发展历史分子动力学模拟的一个关键在于原子势函数的选取原子势早期一般采用简单的对势, 但对势无法正确描述弹性常数, 其结果不理想世纪年代提出的镶嵌原子法、有效介质理论更客观地反映了原子间多体作用的本质, 可得到较合理的结果.认为体系总能量为
分子动力学模拟方法的基本原理与应用 摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如Lennard-Jones势; 论述了几种常用的有限差分算法, 如Verlet算法; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。 关键词: 分子动力学模拟; 原子间相互作用势; 有限差分算法; 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 从统计物理学中衍生出来的分子动力学模拟方法(Molecular Dynamics Simulation, MDS) , 实践证明是一种描述纳米科技研究对象的有效方法, 得到越来越广泛的重视。所谓分子动力学模拟, 是指对于原子核和电子所构成的多体系统, 用计算机模拟原子核的运动过程, 从而计算系统的结构和性质, 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段, 称之为“计算机实验”手段, 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。 科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在材料研究中显得非常有吸引力。 分子动力学模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈
分子动力学模拟： 对于原子核和电子组成的多体体系，求解运动方程（哈密顿，牛顿，拉格朗日），用经典和量子化方法求解粒子的运动状态。 MC方法：系综（抽样）平均法分子动力学：时间平均 一优点：遇到的不利影响因素回避掉,从而达到实验研宄难以实现的控制条件。核心算法：粒子的运动状态就必须把运动方程离散化，离散化的方法有经典Verlet算法、蛙跳算法(Leap-frog)、速度Veriet算法、Gear预估-校正法等。缺点：元胞体积和形状不变，不含有自由电子，对金属体系计算不理想。 注意：一般而言，MD模拟时间足够长，初始条件不会影响计算结果，但是会加大构型平衡的计算时间。 二步骤： 1.选取所要研究的系统并建立适当的模拟模型。 2.设定区域的边界条件，选取粒子间相互作用势模型；要注意观察PBC边界条 件的使用，以及计算格子和建模的晶格子之间的关系。体系是单胞沿不同方向重复叠合而组成。但模拟时只保留基本单元，由平移对称矩阵计算得到其他原子的空间坐标。最小近邻的截断半径。 3.设定系统所有粒子的初始位置和初始速度； 4.计算粒子间相互作用力和势能，以及各个粒子的位置和速度；最好与实际模 型相符，以减少达到平衡的时间。势场参数调整，最小近邻的截断半径。 对势：LJ势（惰性气体，过渡金属，柔韧材料），Born-lande势（离子晶体），Morse势，Johnson势（金属）发展到三体势，缺点是导致Cauchy关系，即不能描述晶体的弹性性质。 多体势：80年代以后，EAM势等（晶体对势+核嵌入电子云嵌入能），多用于金属。 5.待体系达到平衡后，构型积分获得体系的宏观性质。选取合适的系综，控制
vasp的分子动力学模拟 VASP 2010-01-15 02:26:36 阅读57 评论0 字号：大中小 vasp做分子动力学的好处，由于vasp是近些年开发的比较成熟的软件，在做电子scf速度方面有较好的优势。 缺点：可选系综太少。 尽管如此，对于大多数有关分子动力学的任务还是可以胜任的。 主要使用的系综是NVT 和NVE。 下面我将对主要参数进行介绍！ 一般做分子动力学的时候都需要较多原子，一般都超过100个。 当原子数多的时候，k点实际就需要较少了。有的时候用一个k点就行，不过这都需要严格的测试。通常超过200个原子的时候，用一个k点，即Gamma点就可以了。 INCAR： EDIFF 一般来说，用1E-4 或者1E-5都可以，这个参数只是对第一个离子步的自洽影响大一些，对于长时间的分子动力学的模拟，精度小一点也无所谓，但不能太小。 IBRION=0 分子动力学模拟 IALGO=48 一般用48，对于原子数较多，这个优化方式较好。 NSW=1000 多少个时间步长。 POTIM=3 时间步长,单位fs, 通常1到3. ISIF=2 计算外界的压力.
NBLOCK= 1 多少个时间步长，写一次CONTCAR，CHG和CHGCAR，PCDAT. KBLOCK=50 NBLOCK*KBLOCK 个步长写一次XDATCAR. ISMEAR=-1 费米迪拉克分布. SIGMA =0.05 单位:电子伏 NELMIN=8 一般用6到8, 最小的电子scf数.太少的话,收敛的不好. LREAL=A APACO=10 径向分布函数距离, 单位是埃. NPACO=200 径向分布函数插的点数. LCHARG=F 尽量不写电荷密度,否则CHG文件太大. TEBEG=300 初始温度. TEEND=300 终态温度。不设的话，等于TEBEG. SMASS -3 NVE ensemble;-1 用来做模拟退火。大于0 NVT 系综。 【转】vasp的分子动力学模拟 ★★★★★★★★ 小木虫(金币+1):奖励一下，谢谢提供资源 uuv2010(金币+1): 您是否可以做一个专题，详细讲讲怎么做？比如第一步需要干什么，第二步需要干什么，结果怎么分析……如果能做一个这样完整的专题就太好了，不知道您是否有兴趣？2011-07-13 18:20:12 uuv2010(金币+1): 多谢提供资源2011-07-16 17:39:55 uuv2010(金币+5, 1ST强帖+1): 多谢您的详细讲解！感谢就此专题与大家分享！2011-08-12 18:25:12 vasp做分子动力学的好处，由于vasp是近些年开发的比较成熟的软件，在做电子scf速度方面有较好的优势。 缺点：可选系综太少。 尽管如此，对于大多数有关分子动力学的任务还是可以胜任的。 主要使用的系综是NVT 和NVE。 下面我将对主要参数进行介绍！
1 分子动力模拟计算的基本原理 分子动力计算的基本原理，即为利用牛顿运动定律。在分省储存空间。其缺点是位置与速度不同步。这意味着在位置一子动力模拟中，体系原子的一系列位移是通过对牛顿运动方程定时，不可能同时计算动能对总动能的贡献。 的积分得到的，结果是一条运动轨迹，它表明了系统内原子的位置与速度是如何随时间而发生变化。 先由系统中各分子位置计算系统的势能，按照经典力学，系统中任一原子i 所受的力为势能的梯度： 将牛顿运动定律方程式对时间积分，可预测i 原子经过时间t 后的速度与位置： 式中， 及 分别是粒子i 的位置与速度，上标“0”为各物理量的初始值[1]。 2 牛顿运动方程的数值解法 为了得到原子的运动轨迹，必须解式（3）的牛顿运动方程，可采用有限差分法。有限差分法的基本思想就是将积分分 成很多小步，每一小步的时间固定为 。常用的有以下几种算法：① Verlet 算法；② Velocity-Verlet 算法；③ leap-frog 算法（蛙跳算法）；④ Beeman 算法；⑤ Gear 算法。leap-frog 算法和Gear 算法由于使用简便，准确性及稳定性高，节省储存空间等作者：photon 优点，已被广泛采用。 2.1 leap-frog算法 Leap-frog 算法速度与位置的数学式为： 为了执行leap-frog 算法，必须首先由t-0.5 时刻的速度与t 时刻的加速度计算出 ，然后由方程 计算出位置 。时间为t 时的速度可由式（6）算 2.2 Gear算法[1] Gear 所提出的一种利用数值解的方法，称为校正预测法（predictor-corrector method ）。时间t+ 时的位置、速度等可由时间t 的泰勒展开式预测得到： 式中的 的1次、2次、3次微分。式（7）所产生的速度、加速度不是由牛顿运动方程解得的，所以并非完全正确。可由所预测的位置 计算所受的力及正确的加速度 。设正确的加速度与预测的加速度之间的误差为： 式中， 均为常数。以上为Gear 的一次预测校正法，也可将此计算推展至更高次的校正。 3 势函数 势函数表明了原子间的相互作用。针对不同的计算物质，不同的模拟目的，势函数有不同的形式。计算结果的可靠性与势函数密切相关。在分子动力学发展初期，主要采用对势。随着模拟体系的复杂性，逐渐出现了多体势，以弥补对势的不足。 3.1 对势 主要是Lennard-Jones 势（L-J 势），又叫12－6势能，它的数学表达式是： 式中，r 为原子对间的距离， 、 是势能参数。在 L-J 势能中， 项是排斥项， 项是吸引项。当 r 很大时，L-J 势能趋近于零，表示当原子对相距很远时，彼此之间已经没有非键 分子动力学模拟方法概述 周晓平 田壮壮 忽晓伟 （郑州大学 西亚斯国际学院 河南 郑州 451100） 摘 要： 主要介绍分子动力学模拟的基本原理，阐述分子动力学方法的运动方程、数值解法、势函数、边界条件、适用系综以及体系相关性质的计算。最后指出分子动力学模拟方法的优势和发展方向。 关键词： 分子动力学；势函数；边界条件 中图分类号：O414 文献标识码：A 文章编号：1671－7597（2012）1210040－02 由牛顿第二定律可得i 原子的加速度为： 可得各量的校正式为：