4.1基本原理
分子动力学模拟（moleculardynamicssimulations）是一种用来计算一个经典多体体系的平衡和传递性质的方法。这里的经典意味着组成粒子的核心运动遵守经典力学定律，这对于许多体系是一个很好的近似。只有当处理到一些较轻的原子或分子(He，H2，D2)的平动、转动或振动频率满足h>kBT的振动时，才需要考虑量子效应。
考虑含有N个分子的运动体系，系统的能量为系统中分子的动能与系统总势能的和。其总势能为分子中各原子位置的函数(，，，)12nUrrr。通常势能可分为分子间（或分子内原子间）的范德华作用（vdw）与分子内部势能（internal，int）两大部分：
vdwintUUU（4-1）范德华作用一般可将其近似为各原子对间的范德华作用的加和：
依据经典力学，系统中任一原子i所受之力为势能的梯度：
由牛顿运动定律可得i原子的加速度为：
将牛顿运动定律方程式对时间积分，可预测i原子经过时间t后的速度与位置：
i式中，v为速度，上标“0”为各物理量的初始值。
分子动力学计算的基本原理，就是利用牛顿运动定律。先由系统中各分子的位置计算系统的势能，再计算系统中各原子所受的力及加速度，然后令上式中的tt，则可得到经过t后各分子的位置及速度。重复以上的步骤，由新的位置计算系统的势能，计算各原子所受的力及加速度，预测再经过t后各分子的位置及速度……如此反复循环，可得到各时间下系统中分子运动的位置、速度及加速度等资料。一般将各时间下的分子位置称为运动轨迹(trajectory)。
4.2牛顿运动方程式的数值解法
在分子动力学中必需求解(4-5)式的牛顿运动方程以计算速度与位置，有限差分方法(finitedifferencemethods)常被用来解具有连续势能函数的分子动力学轨迹，也就是速度与位置。其基本思路就是把积分划分为许多小段，在时间的间隔上为固定的t，在t时刻，粒子受到的总作用力是其和其它粒子相互作用的向量之和。从粒子受到的力，我们可以确定粒子的加速度，结合t时刻的位置和速度得到t+t时刻新的位置和速度。这里假定在时间步长t范围内粒子受到的力是不变的。再计算在新的位置粒子受到的力，继而得到t+2t时刻的新位置和速度，这样周而复始的计算就可以得到粒子运动的轨迹。
利用有限差分方法来积分运动方程的算法有好多种，其中有许多在分子动力学的计算中被广泛地应用。在所有的算法中都假定位置和动力学性质（如速度、加速度等）可以通过Taylor展开式表示：
其中，v(t)是速度（位置对时间的一阶微分）；a(t)是加速度（二阶微分），b(t)
是三阶微分，等等。Verlet算法[193]是在分子动力学模拟中最广泛使用的积分运动方程的算法。Verlet算法利用t时刻的位置和加速度，以及前一时刻的位置r(tt)，来计算新时刻tt的位置r(tt)，我们给出这些性质和速度在t时刻的相互关系。
把这两个方程相加：
r(tt)2r(t)r(tt)t2a(t)（4-11）在Verlet积分算法中不显含速度，速度可以用很多方法来计算，如一个最简单的方法：
v(t)r(tt)r(tt)/2t（4-12）另外，速度还可以利用半个步长ttVerlet算法的应用是很直接的，并且所需的内存也不是很多，包括两组位置〔r(t)和r(tt)〕以及加速度a(t)。Verlet算法的缺陷之一是位置r(tt)的确定利用了较小项〔t2a(t)〕和两个较大项[2r(t)和r(tt)]的差值，这将导致计算精确性的降低。
而且，在方程中，Verlet算法缺乏精确的速度计算项使它获得速度很困难，并且速度实际上是在下一步位置确定后才能得到。另外，Verlet算法不是自我开始的（self-start），新的位置是通过当前时刻位置r(t)和前一时刻位置r(tt)获得的，在t=0时刻，很显然只有一组位置，因而利用其它一些方法获得tt的位置是很必要的。一种获得r(tt)的方法是利用Taylor展开，如方程（4-6）在第一项截断，则r(t)r(0)tv(0)。
随后的跳蛙法（leap-frog）[193]是对Verlet算法进行了一些改进：
应用跳蛙算法，首先要利用tt这样就“跳过”位置给出了tt21时刻的速度，位置则跳过速度给出了tt时刻的新位置，从而为tt23时刻的速度做了准备，等等。相对于传统的Verlet算法，跳蛙算法有两个优点：它显含了速度的计算，而且不需要计算两个较大数的差值。然而，它的明显缺点就是位置和速度并不是同时确定的，这就意味着在确定位置的同时计算动能对总能量的贡献是不可能的，而位置的确定就同时确定了势能对总能量的贡献。
速度Verlet(VelocityVerlet)算法[194]同时给出了位置、速度和加速度，而没有任何精确度的损失：
Beeman算法对于速度给出了更加精密的表达式，因为动能是直接从速度计算得到的，所以它通常给出较好的能量守恒。然而它的表达式比Verlet更加复杂，在计算过程中就需要更多的计算机时。
4.3周期性边界条件
正确模拟边界或边界作用对于模拟方法是极为重要的，因为边界条件可以利用相对较小的分子来计算“宏观”性质。边界条件的重要性可以利用下面简单的性质来加以说明。假设我们有1个1L的容器，在室温条件下充满水，这个立方体容器中含有水分子的个数为3.31025个，水分子与器壁的相互作用可以延伸到10个分子直径的大小，水分子的直径大概为2.8，与边界碰撞的水分子的个数为21019。所以，在150万个水分子中只有1个受到容器壁的影响。在MC和MD模拟中，粒子的个数要远远小于10251026，并且通常要小于1000个。在一个不超过1000个分子的体系中，它们都不包含在和容器壁相互作用的粒子中。显然，用容器中1000个水分子来获得“宏观（bulk）”性质是不太正确的。可供选择的方法是去掉（dispensewith）容器。
执行分子动力计算通常选取一定数目的分子，将其置于一个立方体的盒子里。
但执行的计算必需与实际的体系相符，通常实验的系统密度为必须满足的条件。设立方体盒子的边长为L，则其体积为V=L3。若分子的质量为m，则N个分子系统的密度为：
计算系统的密度应等于实验所测定的密度，以此作为调整盒长的依据。以水分子系统为例，假设执行含有1000个水分子的动力学模拟计算，水的密度为1g/cm3。则L=3.1010-7cm=31（4-23）为使计算中系统的密度维持恒定，通常采用周期性边界条件（periodboundarycondition）。周期性边界条件使利用较少粒子进行宏观性质的模拟成为可能。在周期性边界条件下，粒子所受到的力和宏观液体中粒子受到的力是一样的。考虑立方体盒子中的粒子，它在各种方向重复以给出周期性排列。图4-1给出了一个二维盒子，其中，每个盒子周围有8个邻近盒子，这样在三维中每个盒子周围就有26个近邻。
镜像盒子中粒子的坐标可以简单地用加上或减去盒长的整数倍计算。在图4-1中给出，如果模拟过程中一个粒子离开盒子，则从相对应的边上有一个粒子进入盒子，因而在整个模拟过程中体系的总粒子数保持不变，即密度不变，符合实际的要求。尽管在计算机模拟中已广泛地使用周期性边界条件，但是它确实还存在着一些不足。周期性晶胞的一个很显着的不足是它不能获得大于晶胞长度的波动，比如接近于液-汽临界点。
4.4截断半径与最近镜像
MD模拟中最耗时的部分是能量最小化中非键相互作用和（或）力的计算。在力场模拟中键长伸缩，键角弯曲和二面角扭曲的个数正比于原子的个数，而非键相互作用项的个数随着原子个数平方的增加而增加，因此它的数量级为N2。原则上，非键相互作用在体系中应当对每两对原子间都进行计算，然而在模型中这一点还没有完全实现。Lennard-Jones势能函数随着距离的衰减很快，在2.5时的Lennard-Jones势大概是1.0时的1%。解决非键相互作用最普遍的方法就是利用截断半径和最小镜像方案。在最小镜像方案中，每个原子最多只能“看到”体系中每个原子的一个镜像（通过周期性边界条件无限重复）。如图4-2所表示，只需要计算距离最近原子或镜像之间的能量和（或）力就可以。比如，计算分子1与3的作用力，是取与分子1最近距离的镜像分子3。在所有镜像系统中，分子1和3距离最近的是模拟系统中分子1与D盒中的分子3，而非模拟系统中的分子1与分子3。同样的，计算分子3与分子1的作用是取模拟系统中的分子3与E盒中的分子1。
由于在计算中利用最近镜像的概念，因此就需要采用截断半径（cut-offradius）的方法来计算远程相互作用力。当利用截断半径时，如果原子对之间的距离大于截断半径，则它们之间的相互作用为零，当然同时也要考虑最近镜像。如果利用周期性边界条件，则截断半径不应太大，否则它将看到自身的镜像，那么同一个原子将被计算两次。这就暗示，在模拟立方体盒子中的液态体系时，截断半径不能大于立方盒子的一半，长方体晶胞中截断半径不能大于长方体最小边长的一半。对于分子的模拟，截断半径的上限受到分子大小的影响。模拟中非键相互作用如果只考虑Lennard-Jones势能函数，截断半径就不应过大，用2.5的截断半径就会带来一定的误差。然而，当考虑长程静电相互作用时，截断半径必须要大，事实上证明此时利用多大的截断半径都是有误差的，一般的原子所选取的静电相互作用的截断半径约为10左右，但这只是一个简单的估计，对于不同的原子是不能一概而论的。
一般的非键势能函数在所选取的截断半径处势能值并不为零。因此，分子动力学计算中就会出现能量不连续的情况。处理这种问题的方法，通常将势能函数乘上一个开关函数(switchingfunction)以弥补这一缺点。设截断半径为Rc，开关函数的起始点为Rs，则该函数S(r)之形式为：
