第四讲 运动方程的积分算法

第四讲运动方程的积分算法黄敏生积分运动方程的注意点•何为一个好的运动方程积分算法？但这并不是很贴切，因为花在积分运动方程上的时间分率（相对于计算相互作用来说）很小，至少对原子或简单分子体系是这样。计算速度？？准确度对较大的时间步长来说更重要，因为所能使用的时间步长越长，单位时间内力的计算量越少。因此，这意味着采用允许用长时间步长的成熟算法是有利的。（然而，保证准确度，不允许发散）准确度？？积分运动方程的注意点要想获得允许使用长时间步的算法，必须将信息存储在粒子坐标的较高阶导数中。结果是这需要更多的内存。对于一个通常的模拟来说，这并不是一个严重的缺点。因为除非处理很大的体系，存储这些导数所需的内存与即便是在一个通常的工作站上所能获取的总量相比来说仍是很小。节省内存？？能量守恒是一个重要的判据，但实际上需要区分两种能量守恒，即短时间的和长时间的。复杂的高阶算法通常在短时间内（如在几个时间步内）有很好的能量守恒性，然而该方法通常会有所不期望的特征，即时间较长时总能量漂移。（不一定最好）能量守恒？？积分运动方程的注意点最好有一个算法能同时在短时间和长时间准确地预测所有粒子的轨迹。不存在MD模拟所研究的所有体系，体系的轨迹穿过相空间（即对于用由所有粒子坐标和动量所跨过的6N维空间）时敏感地依赖于初始条件。这意味着两个初始靠得很近的轨迹随时间的演绎将会显著分开。李雅普诺夫不稳定性MD获得的轨迹在某种意义上与真实的轨迹相接近。MD的目标并不是精确地预测一个已知初始条件的体系轨迹将会发生什么（卫星轨道预测）。而对统计预测感兴趣。在MD中，统计预测是足够精确的。标准Verlet算法•Verlet提出的Verlet算法在分子动力学中应用最为广泛。由Taylor公式展开有：位移+标准Verlet算法由Taylor公式展开有：速度在典型分子动力学模拟中，只有原子的初始位置和初始速度是给定的，而在verlet算法中，计算下一步的原子的位置需要前两步的信息。那么如何采用verlet算法计算第二个时间步的原子位置呢？不出现在算法中。-•解决方法一是，在第二个时间步中，把该步力看成常量，且使用普通运动学方程:标准Verlet算法•解决方法二是trrtriii)0()0()(•流程标准Verlet算法标准Verlet算法加速度项Leap-frog算法•对标准Verlet算法进行改进而得到的蛙跳算法(Leapfrogmethods)。相比标准Verlet算法，它有两个优点。•一、包含显式速度项。•二、计算量稍小。•缺点原子的位置与速度计算不同步，这就意味着在确定位置时，不能同时计算体系的动能，给模拟过程带来不便。Leap-frog算法速度Verlet算法速度Verlet算法Velocity-Verlet算法不仅可以获得相同精度的原子位置和速度量，给出了显式的速度项，而且在每步积分中只需要存储一个时刻的状态变量，模拟稳定性好，允许采用较大的时间步长，计算量适中，因而在分子动力学方程的积分算法中得到了最广泛的应用。Gear预测校正法预测校正法是分子动力学模拟中的常用算法之一，其基本思想是Taylor展开，这种算法包含三个部分。第一步，利用泰勒展开预测下一时刻的位置及其一阶、二阶、三阶导数：式中v、a、b、c分别是位置矢量r的一阶、二阶、三阶和四阶导数。Gear预测校正法第二步根据新的原子位子rp，计算受力以及修正加速度ac(t+t)。定义预测误差第三步根据加速度的预测误差对各预测量进行修正：•预测－校正法允许的时间步长比其它算法长两倍以上•每个积分步内要计算两次体系势能，以得到原子间相互作用力。•该算法的稳定差，能量波动较大，较verlet算法占用更多的内存。Gear预测校正法时间步长Toolong-errorsresultfromapproximationsJustright-errorsacceptable,maximumspeedTooshort-computationneedlesslyslow时间步长过长的时间步时间步长过大，原子的作用力急剧改变。误差逐渐累计，导致结果发散.两个Argon(氩)原子在两个不同时间步长ts的模拟。图中画出的是计算模拟值与理论值的差。RMS(ReservoirModelingSystem保守系统的模拟)EnergydeviationCircles:VerletSquares:Gear4thorderTriangles:Gear5thorderDiamonds:Gear6thorder(log/logscale)误差与时间步长的关系时间步长•时间步长与研究对象、系统温度、所采用的数值积分算法及势能函数有关。不存在一个通用的时间步长值。•一般情况下，体系的温度越低，允许采用的时间步长越大；而模拟较高温度时必须采取较小的时间步长。•一般认为，时间步长应小于原子振动周期的十分之一，而通常原子振动周期的数量级为0.1皮秒(10-12s)，即时间步长应选择在飞秒级(10-15s)。宏观物理量的计算•分子动力学计算最终得到的是系统各个时刻的相空间轨道(thephase-spacetrajectory)，包括任意时刻所有原子的坐标和速度，这些都是微观原子层次的物理量。•一般来说，在探讨其力学行为时，我们需要研究一些宏观力学概念的物理量，比如温度、能量、压强、应力状态等，这就需要我们对分子动力学计算出来的粒子数据进行分析。•统计力学是连接微观层次的物理量和宏观概念物理量的桥梁，利用统计力学原理我们可以从系统中单个粒子的运动学状态得到整个粒子系统的一些性质。•系统的物理性质是系统中粒子坐标和速度的函数，对于任意一个时刻宏观概念物理量A，定义为宏观物理量的计算统计平均•势能部分可以按势函数计算，对于不同的原子势函数表达式有不同的计算公式。•动能计算公式为宏观物理量的计算能量温度温度T直接与粒子动能相关，即著名的均匀分布公式，每个自由度赋予kBT/2的能量,N个粒子的总自由度为3N，故动能为计算不同热力学状态的总能E和温度，可以得到内能一温度曲线E(T)。这对于监测相变的发生非常有用，相变发生时，该曲线会有跳跃。•应变是表示物体变形大小的测度，可以有不同的定义方式。•应变是一个相对量，反映了物体相对初始构型的变形程度。•应变又是纯几何量，与空间尺度变化无关。•因此宏观、纳观尺度下应变可以采用相同的定义方式原子应变原子应力宏观应力反映了单位面积上作用力的大小，是关于面积的强度量。原子应力离散原子系统的原子应力是关于体积的强度量。实际上，原子应力只是形式地沿用了应力的概念，具有与宏观应力完全不同的特征。原子应力表明一个原子与周围原子相互作用的强弱程度原子应力动量势能二阶对称算子Hamilton自由能可以看出原子应力具有“能量密度”的量纲。包括了原子动量流和原子间作用力的贡献。是原子的一种力学“活性能”，反应了原子产生运动的潜在能力。原子应力越高，则原子越容易发生位错运动。原子应力Voronoi几何构形的数学描述原子应力Lammps原子应力的定义原子对势相互作用原子键作用键角作用二面角作用面弯曲作用约束作用Lammps中原子应力没有除于体积，是能量的量纲。（e.v）•宏观应力实际上是原子应力在一定程度上的统计平均。•如果要计算体系中某个区域（由region定义，可以是整个模拟盒）所围成的“块”的应力，在Lammps只需将该区域里的所有原子的单原子应力值加起来，再除以这个区域的体积即可，无须进行单个原子体积的计算。原子应力原子应力与宏观应力的关系Lammps中原子应力谢谢