分子动力学模拟方法及其在激光加工中的应用
摘要：介绍了分子动力学模拟的基本原理以及常用的概念，论述了几种常用的有限差分算法，分析和总结了分子动力学模拟的功能、特点和应用，列举了在激光加工中分子动力学模拟方法的一些应用，并在材料科中的应用情况进行了简要概述．
关键词：分子动力学模拟有限差分法飞秒激光
0 前言
随着计算机模拟技术的发展，使实验上尚无法获得或很难获得的大量重要信息的获取成为可能，虽不能完全代替实验但为科研工作者们提供了重要的参考、指导实验、验证某些理论假设，降低试验的盲目性、成本低廉广等，其中特别是分子动力学模拟在各个学科中都有着广泛而重要的应用。分子动力学模拟(moleculardynamics simulation，MD)【1】。所谓分子动力学模拟，是指对于原子核和电子所构成的多体系统。用计算机模拟原子核的运动过程．从而计算系统的结构和性质。其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下接牛顿定律运动【2】，它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段，称之为“计算机实验”手段，在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。是在评估和预测材料结构和性质方面模拟原子和分子的一种物质微观领域的一种重要模拟方法，通过计算机对原子核和电子所构成的多体体系中的微观粒子之间相互作用和运动进行模拟，在此期间把每一原子核视为在全部其他的原子核和电子所构成的经验势场的作用下按照牛顿定律进行运动，进而得到体系中粒子的运动轨迹，再按照统计物理的方法计算得出物质的结构和性质等宏观性能。简而言之即是应用力场及根据牛顿运动力学原理所发展的一种计算机模拟方法。
分子动力学模拟是一种非常有效的计算机技术已成为重要的科学研究的方法之一。
1 分子动力学方法的基本原理
计算中根据以下基本假设【3】：
(1) 所有粒子的运动都遵循经典牛顿力学规律。(2)粒子之间的相互作用满足叠加原理。显然这两条忽略了量子效应和多体作用，与真实物理系统存在一定差别，仍然属于近似计算。从分子动力学的定义可以看出要想理解什么是分子动力学模拟，就必须首先清楚地理解力场、牛顿运动方程及其数值解法等基本概念。同时，在分子动力学模拟领域中，系综、周期性边界条件、积分步长等也是经常提及的术语名词，对它们的正确理解也影响着对分子动力学的深入理解。
1.1 力场
力场就是势能面的表达式，它是分子动力学模拟的基础，是分子的势能与原子间距的函数，针对特定的目的，力场分为许多不同形式，具有不同的适用范围和局限性，计算结果的可靠性与选用的力场有密切关系。在各种形式的力场中，Lennard-Jone(LJ)势能是目前较为常用的势能．其势能表达式为：
U(r)为对应于r值下的分子的势能；r为原子间距；ε，δ为势能参数。
在众多科学家的努力之下，力场已由最初的单元予分子系统发展到多原子分子、聚合物分子甚至生物分子系统。力场的复杂性、精确性、适用范围都有了很大的进步。在诸多力场中，每个力场都有着各自的优缺点及其适用条件。因此，在模拟时应该对当时模拟的条件、系统的特征等诸多因素加以分析选取适合的力场，才能保证模拟的速度和准确性。
1.2 牛顿运动方程及其数值解法
在分子动力计算中必须先解以下牛顿运动方程：
根据计算结果再算出粒子的速度与位置，从而确定粒子运动的轨迹。这是分子动力学模拟计算的基本思路。关于牛顿运动方程的解法有很多，一般采用常用Verlet所发展的数值解法，其中最早的Verlet方法是将粒子的位置以泰勒式展开，经过计算得出结果，由于该法容易导致误差本文不再详细介绍。之后Verlet为解决这个问题发展出了跳蛙方法(1eap frog method)，此方法计算速度与位置的数学式是：
则时间为t的速度由一下公式算出：
可以看出，该算法只需与两个已知条件，节省了计算机的存储空间，具有较高的准确性和稳定性。现今，该方法已广泛地应用于分子动力学模拟中。
1.3 一些重要术语
系综系综(ensemble)是指具有相同条件系统(system)的集合．如正则系综(canonical ensemble)是指具有相同分子数目N、相同体积V与相同温度T的系统的集合，符号为(N，V，T)，其他还有等粒子等温定压系综(N，T，P)、等粒子等容等系统的能量系统(N，V，E)等多种系综。
系综是统计力学中非常莺要的概念，系统的一切统计特性基本都是以系综为起点推导得到的。实际应用时，要注意选择适当的系综，如(N，T，P)常用于研究材质的相变化等。
周期性边界条件分子动力学计算通常是选取一定数目的分子，将其置于一
个立方的盒子中，该盒子即为模拟系统，周围是与它具有相同的粒子排列和运动的盒子。在粒子的运动过程中，计算系统中若有一个或几个粒子跑出盒子，则必有一个或几个粒子由其他盒子跑进该计算系统以维持模拟系统中的粒子数为定植从而保证该模拟系统的密度恒定，才能符合实际状况。这种为保证体系密度恒定而设定的条件称为周期性边界条件。
积分步长积分步长即为分子动力学计算公式中的8t(integration time step)，它的选取决定了模拟的时间和准确性，积分步长越小准确性越高但越费时，相反积分步长越长计算速度越快但会降低计算的准确性，所以节省计算时问又不失去其精准性是选取适当的积分步长的原则。一般取系统最快运动周期的十分之一。
2 分子动力学模拟有限差分法
为了得到原子的运动，可以采用各种有限差分法来求解运动方程。常用的有以下几种算法：
Verlet算法【4】：Verlet提出的Verier算法在分子动力学中运用最为广泛．它运用原子在t时刻的位置r(t)和加速度a(t)及t -δt时刻的位置，计算出t +δt时刻的位置
这里为简单计，省略了下标i。显然，速度并没有出现在Verlet算法中，计算速度有很多种方法，一个简单的方法是用t+δt时刻与t-δt时刻的位置差除以2δt。Verlet算法执行简明，存储要求适度，但它的一个缺点是位置r(t+5t)要通过小项与非常大的两项2r(t)与r(t-δt)的差的相加得到，这容易造成精度损失。Verier算法中没有显式速度项，在下一步的位置没得到之前，难以得到速度项。另外，它不是一个自启动算法；新位置必须由t时刻与前一时刻t-δ的位置得到，在t=0时刻，只有一组位置，所以必须通过其他方法得到t-δt的位置。获得t-δt时刻的位置方法之一是应用近似式
Velocity-Verlet算法【5】：Swope提出的velodty-verlet算法可以同时给出位置、速度与加速度，并且不牺牲精度。这种算法的优点是给出了显式速度项，并且计算量适中，目前应用比较广泛。
Leap-frog算法【6】：Hockney提出的Leap—frog算法，是Verlet算法的一种变化。Leap-frog算法相比Verlet算法有两个优点：(1)包括显速度项，(2)计算量稍小。它也有明显的缺陷：位置与速度不是同步的。这意味着在位置一定时，不可能同时计算动能对总能量的贡献。
Beeman算法【7】：Beeman提出的Beeman算法也和Verier算法有关。Beeman算法运用了更精确的速度表达式。因为动能是直接由速度计算得到的，所以它更好地保持了能量守恒，然而由于它的表示式比Verlet算法复杂得多，所以计算量较大。
Gear算法【8】：Gear提出了基于预测一校正积分方法的Gear算法，这种方法可分为二步；首先，根据Taylor展开式，预测新的位置，速度与加速度．然后，根据新的位置计算得到的力计算加速度a(t+δt)，这加速度再与由Taylor级数展开式预测的加速度a(t+δt)进行比较，两者之差在校正步里用来校正位置与速度项，Gear的预测一校正算法需要的存储量为3(O+1)N，O是应用的最高阶微分数，v是原子数目，本例的存储量为15M而Verlet算法的存储量为9N。更重要
