分子动力学入门第二章

1、第二章：分子模拟的基本部分2.1模拟的物理体系模拟的最主要组成部分就是所研究物理体系的模型。对于分子动力学来讲就是选择势能函数：V（r，.r）该函数是有关原子核位置的函数，它表示当所有原子的位置组成一特定构型时体系的势能。该函数是原子的平动和转动的不动点，通常的位置是指原子间的相对位置而不是绝对位置（内坐标表示，而不是笛卡尔坐标）。原子所受到的力就是势能相对于位移的梯度：F=-V（r，.r）（1）。这种形式暗示存在一种有关总能量E保守的定律，E=K+V，这里的K值得是瞬间动能。最简单的势能函数V的写法是表示成成对相互作用的和：V（r，.r）=（2）第二个求和下的j>i的目的是考虑没对原子

2、只能求和一次。在以前大多数势能函数都是有成对的相互作用构成的，但是现在情况不在是那样啦。现在已经知道towbody 近似对一些相关系统非常不合适，例如金属和半导体。许多种many-body势能函数在凝聚态模拟中普遍得到运用，这会在第四章简单的做一介绍。寻找精确的势能函数也是目前重要的一个研究领域。在第四章会介绍一些目前有关这方面的研究情况。现在我们来看看目前最普遍运用的相互作用模型：LennardJones的成对势能函数。2.1.1 LennradJones势能函数LennradJones的公式：（r）=4（）-（）（3）该函数表示一对原子间的势能，而总势能是有（2）决定。该势能函数在很大r处

3、具有一个“attractive tail”（相互吸引），r能达到的最小为1.122，在很短距离能强烈排斥，在r=处相互作用为0，随着r的减小渐渐增大。1/r,在短程起主导作用，模拟当两原子间靠的的非常近时的原子间的排斥作用。它的物理源于与泡利原则：当原子周围的电子云出现重合时，体系的能量骤然增加。指数12的选择是有一个实际的依据：等式（3）是非常容易计算的。实际上，站在物理学层面上，指数函数是更加合适的1/r在远距离起主导作用，构建相互吸引部分。该部分会给体系一个内聚力。1/r的引力范德华色散力，源于偶极子的相互作用（波动的偶极子相互作用），范德华引力是非常弱的，但是它主导了封闭体系的bond

4、ing 特性。这些证据就是LJ势能函数模拟的非常好的原因。参数和是选择拟合材料的物理特性。另一方面，LJ势能函数不能充分模拟开放体系的，因为存在局部强的键（在共价键体系）或者存在自由的“电子海洋”（金属）。在这些体系中two-body模型本身就很不适合。有关这些体系的势能函数会在第四章讲解。但是，无论LJ模拟实际材料到底怎样，目前LJ12-6势能函数仍是模拟体系的极其重要部分。有很多文章利用LJ在不同几何构造中（固体，液体，表面，成簇，二维系统等）研究了原子相互作用的行为。在一些基本问题上，可以说LJ是标准的势能函数。模拟通过LJ体系可以帮助我们（现在仍然是）理解凝聚态物理许多领域的基本问题，

5、正由于这样LJ函数的重要性是不可估计的。（通常习惯当用LJ势能函数进行模拟时，和，许多示例代码都遵循这个惯例）。2.1.2 截断和远程相关公式(3)中的r的范围是无限的。但是在实际运用过程中，习惯上要选择一个截断半径Rc,当原子间的距离超过Rc时，默认原子间的相互作用为0.这样会简化程序和减少计算资源。但是简单的势能截断会引起新的问题：无论什么时候一个粒子对“经过”截断距离，能量都有一个小的变化。在模拟中大量的这样的事件可能破坏能量守恒。为了避免这样的问题，我们必须要强制使势能在截断半径以外消失： V（r）=势能截断当然会影响物理量的大小。关于截断的影响我们可以近似估计其大小，即把Rc之后的体

6、系看成是一个均匀体系（密度恒定）。对于大体系（沿着每个方向都有周期性边界条件），这就相当于添加一个常数的纠正。例如势能的尾部（引力引起的）对于凝聚能量和总压力的贡献是很小的。对于有些具有自由表面的几何构型，截断效应是很难估计的，这里原因有对称性低，以及像表面能量相当大。普遍用于LJ的截断半径为2.5和3.2。值得一提的是这种LJ截断模型计算去来的值可以作为一般的two-body体系的一个参考。在很多情况下，没必要考虑截断效应的纠正，因为截断模型本身也是研究的的主题。对于金属和半导体的是能函数的经常是包含截断半径设计的。因此这些势能函数是不含真正范德华力的。由于该力相对于金属键和共价键非常弱，因

7、此通常问题不大，但是当几何构造有两个以上独立体时就理当别论啦。2.2 周期性边界条件在模拟体系的边界我们应该做些什么？一种可能是就是什么都不做：在体系的简单的末端，靠近边界的原子相对于内部原子来说具有很少的邻居。换句话说，围绕原子的是表面。除非我们真的是想模拟一簇原子，这种条件是不真实的。无论模拟的体系多大，相比于宏观一小块物体（是10的数量级），模拟体系的原子数量N是微不足道的，表面原子数与总原子数的比值是非常大的（相比现实）。因此表面效应是非常重要的。一个解决这问题的途径就是运用边界条件。在运用边界条件时，粒子是被封闭在盒子里的，我们科想象，这个盒子被无限复制后在笛卡尔坐标的三个方向上展开

8、，充满整个空间。换句话说，当我们的粒子在盒子里的位置为r，我们假设该粒子实际上代表了一系列粒子，他们的位置为： r+la+mb+nc（l，m，n=-）l，m，n是整数，a，b，c是盒子各边的向量。所有这些想象的粒子一起移动，实际上在计算机程序中他们中只有一个。现在的关键点就是，每个粒子i应该被认为不仅与盒子里j粒子相互作用，也应该与邻近盒子里他们的像相互作用。这就是，相互作用可以“通过”盒子的边界。实际上，我们可以容易了解到（a）我们实际上我们从我们体系中消除力表面效应。（b）盒子边界的位置没有什么影响（这就是说，a translation of the box with respect to

9、 the particles leaves the force unchanged）显然，成对相互作用数量增加非常多时，PBC对结果的影响就会很明显。在实际情况中，这是不真实的因为势能函数具有短程相互作用。接下来讨论的最小映像准则会使问题进一步简化，利用PBC会使程序的复杂度。2.2.1最小映像法则我们假设我们所用的势能函数是在有限的范围内：当分离的距离等于或大于截断距离Rc时，两个粒子就不再相互作用啦。我们还假设我们用的盒子在笛卡尔坐标的方向上每条边都大于2Rc。当这些条件都满足时，显然在所有的盒子里的粒子i和一系列所有的周期映像粒子j形成的相互作用中之多只有一个。为了说明这一点，我们假设i

10、与映像j中的j1和j2相互作用，两个映像必须被由一个盒子到另一个盒子的平移向量分开，根据假设长度至少为2Rc。为了和j1和j2两个都相互作用，i应该到两个粒子的距离都应该小于Rc，由于映像粒子之间的距离大于2Rc，所以这是不可能的。因此在这种条件下，我们可以运用最小映像准则：在所用的可能的映像粒子j中，选择最近的，舍弃所有其他的。实际上只有最近的才有相互作用，而其他的是没有相互作用的。这种规则被普遍运用着，大大的简化了MD的程序。当然，我们在PBC有效应的盒子在每个方向的长度至少为2Rc。2.2.2表面几何图形PBC的目的就是为了消除表面效应。然而我们也许也对有表面的条件也感兴趣。显然所有的表

11、面物理就是这类问题。对于表面模拟，经常采用的模型是厚板模型：厚片材料，限定具有两个自由表面。这个模型的获得可以通过保留正交平面沿一个方向移动PBC（经常是沿着Z轴移动）。因而，厚板可以被认为是把xy平面无限复制。如果厚板足够厚，内部可以达到和大体积材料一样。两个表面（上面和下面）可以被看做两个独立的表面。在这种假设下，系统的行为可以近似于具有一个单平面的半个无限系统虽然条件与实际表面实验相接近，但是在模拟中是没办法实现的。还有一些其他条件，研究者发现最好在材料的一端“冰冻”一些层原子固定在完美大晶体的位置而让其他的是自由的。这种条件只有当在同样的分离距离这种“冰冻”效应要比其他自由端的效应小。

12、这种条件在一下例子中非常合适，当表面发生大量的扰动：表面重排，和局部融化。当厚板的冰冻的端被选择，必须关注被冰冻的一些列的层厚度至少为Rc，这是为了确保，没有运动的粒子可以通过冰冻面。我们也可以只沿一个方向创造PBC，这样可以获得一个线性构型，具有沿一个线性轴方向的PBC或者同时移动。后一种情况与原子簇相关。自从出现了分子动力学和蒙特卡罗，原子簇的模拟得到了好的解决。2.3 时间积分算法分子动力学的程序的核心就是时间积分算法，需要随着他的运动轨迹积分相互作用粒子的运动方程。时间积分算法是基于有限差分法，时间被分割成有限的格子，时间步t是连续格子上两点间的距离。知道在时间t时（细节依赖于算法）的

13、位置以及时间依赖的相关量，那么积分会给出在时间t+t时的同样的量的值。通过不断的迭代，就可以追踪系统时间的变化。当然，这是一种近似计算，存在相关的误差。尤其，我们可以区分相关的误差：截断误差，指的是有限差分方法与真实解决办法的间的差异。有限差分通常依赖于在泰勒展开式在一定程度上的截断，因此这是误差命名的原因。误差不依赖于implementation：而依赖于算法。舍入误差，与实现算法的a particular implementation 的误差。例如计算机计算中利用的有效数字。两个误差都可以通过减少而得到减少。对于很大时，截断误差占主导，但是随着的减少快速降低。例如，verlet算法的每步时

14、间积分截断误差与成正比。舍入误差随降低很慢，在无限小时占主导作用。利用64位精确计算可以使舍入误差保持更低（今天大多数的工作平台，利用Fortran语言，双精度）。MD计算的两种流行的算法是Verlet 算法和矫正预测算法。这将会在下面介绍，更多的时间积分方法可以参考下面的文献：下面的文献对相关算法进行了深入的分析：2.3.1 Verlet 算法在分子动力学模拟中，普遍运用的时间积分算法可能就是Verlet 算法。基本想法就是三阶r（t）的泰勒展开，一次在时间上向前，一次在时间上向后。v是速度，a是加速度，b是r相关t的三阶导数，因此：r（t+）=r（t）+v（t）+（1/2）a（t）+（1/

15、6）b(t) +O()r（t-）=r（t）-v（t）+（1/2）a（t）-（1/6）b(t) +O()两式相加： r（t+）=2r（t）+a（t）+O()-r（t-）这是verlet算法的基本形式。既然我们是对牛顿方程的积分，a（t）是力除以质量，因此力也是r（t）的函数： a(t)=-(1/m)V(r(t)我们可以马上看到，当算法以阶进行系统演化，即使三阶导数不明确，出现截断误差。同时这种算法容易实现，准确，稳定，这就解释了分子动力学中使用这种算法普遍性。这种Verlet算法的一个问题是速度不能直接生成。尽管速度对时间演化不是必须的，但是有时对他们的了解是必须的。然而，他们需要去计算动能K，

16、K的计算对测试总能量E=K+V的守恒是必须的。这是一个非常重要的测试来验证MD在正确进行。我们可以通过位置计算速度： v（t）=然而，这种表达式的误差与相关而不是。 为了克服这种困难，许多Verlet算法已经发展了，提高了同样轨迹的精确度，只是存储的变量与时间是不同的。蛙跳算法，虽然没报道，在处理速度方面更好些。一个更好的实施基本的算法叫做速度Verlet策略，位置，速度，加速度在时间t+从在时间t同样的数量中获得： r（t+）=r（t）+v（t）+（1/2）a（t） v（t+/2）=v(t)+(1/2)a(t) a(t+)=-(1/m)V(r(t+) v(t+)=v(t+/2)+(1/2)a(t+)尽管我们需要9N的内存储存3N个位置，速度，加速度，但是我们从不需要同时存储三个量中的任何一个在不同的时间。2.3.2 预测矫正算法预测矫正算法是另一个普遍用于运动方程的积分的方法。这种方法在分子动力学中普遍运用归功于Gear，具体包括以下三个步骤：