分子动力学和模特卡罗计算物理简介 分子动力学原理简介 分子动力学程序计算物理简介计算机模拟与计算物理 在当今科学研究中，计算机模拟已经起到 举足轻重的作用。传统的物理研究可以分为理 论物理和实验物理两大领域。20世纪中叶以来 ，随着计算机技术的飞速发展，计算机在各个 领域得到广泛的应用。物理学家也将计算机模 拟应用到物理学的研究中，逐渐形成了一个新 兴的物理学分支—计算物理。 计算物理已经成为和理论物理与实验物理 相平等的一个新的物理研究领域，成为联系理 论和实验的重要的桥梁。计算凝聚态物理和材料设计三种基本数 值模拟手段： 第一性原理计算（First principles calculation） 经典和量子力学分子动力学模拟（molecular dynamics simulation, MD） 模特卡罗方法（Monte Carlo Simulation, MC）分子动力学简介分子动力学 所谓经典分子动力学模拟，是指对于由原子 核和电子所构成的多体系统，用计算机模拟原子 核的运动过程，并从而计算系统的结构和性质。 通常，分子动力学模拟基于Born-Oppenheimer 近 似（即绝热近似），将其中每一个原子核视为在 其它原子核和电子所提供的势场中按牛顿定律运 动。在分子动力学模拟中，系统中原子的一系 列位型是通过求解牛顿运动方程得到的，即通 过下面的方程组得到各个原子在不同时刻的运 动细节。分子动力学模拟不仅能得到原子运动的细节 ，还能像作实验一样进行各种观察。对于平衡系 统，可以用分子动力学模拟作适当的时间平均来 计算物理量的统计平均。对于非平衡系统，也可 以通过分子动力学模拟来在观察时间内对物理现 象进行直接模拟。特别是一些实验上无法观察的 微观物理细节，分子动力学模拟过程中都可以很 容易的观察到。这种优点使得分子动力学在物理 、化学、材料科学等领域得到广泛的应用。 分子动力学特点相对于第一性原理计算，分子动力学模拟 在采用经验势时，虽然计算精度差一些，但是 由于程序简单，可以模拟含大量原子的系统以 及第一性原理难以或无法计算的动力学以及有 限温度的性质，因而有其独特的优点。 当然，目前结合第一性原理和分子动力学 方法，发展出来的第一性原理分子动力学可以 结合两者的优点，具有很强大的功能。但通常 计算量很大，还不太适用于大规模的系统以及 长时间的模拟。分子动力学模拟现状和发展方向 早期的分子动力学模拟研究受到计算机速 度和内存的限制，模拟的空间尺度和时间尺度 都很小。随着微电子技术的进步，现代计算机 的计算速度按照莫尔定律高速发展，同时价格 也越来越便宜。现代处理器的并行设计和计算 机集群技术的发展也使得大规模和长时间的分 子动力学模拟成为可能。2006 年，在IBM 的蓝 色基因（BlueGene/L）超级计算机上，已经成 功的进行了三千二百亿个原子的分子动力学模 拟。中国科研人员利用“天河一号“创分子模拟计算世界纪录 新华网天津６月１１日电（记者毛振华）中国科学院过程工程研究所科研人员日 前利用全球最快计算机“天河一号”的ＧＰＵ超强计算能力，运行了一项分子 动力学模拟项目，成功创下全球分子模拟计算的世界纪录。 中国科学院过程工程研究所研究员葛蔚在接受新华社记者采访时说，研究所的科 研人员采用“天河一号”上的７１６８颗“英伟达ＧＰＵ”，开展了一项规模 巨大的分子动力学模拟，以了解太阳能电池和半导体行业中常用的晶体硅的微 观行为。为此，他们编写的程序大约有１００００行，其中最核心的计算由约 ２０００行ＣＵＤＡ代码执行，最终取得了每秒１．８７千万亿次单精度浮点 运算的可持续性能。 整个模拟总共用了３个小时左右，达到了深入统计分析所需的时间尺度。通过反 复验证，科研人员确认模拟结果与真实材料相符。目前，科研人员仍在紧张分 析模拟结果和改进程序，最终的计算性能有望进一步提高。 据“天河一号”ＧＰＵ生产方英伟达公司对外发布的消息称，此次模拟的运行性 能是之前最高性能的分子模拟的五倍，而模拟规模也在其两倍以上。这次模拟 勾画出了约为１１００亿个原子的微观行为，而之前此类模拟的纪录是对４９ ０亿个原子达到每秒３６９万亿次浮点运算性能。 业内专家表示，计算机模拟对于新材料的研究与生产具有非常重要的意义，可最 大限度地揭示细节，但成本却比实验要低得多。 据了解，“天河一号”位于国家超级计算天津中心，曾于去年荣膺ＴＯＰ５００ 组织颁发的全球最快超级计算机称号。目前被广泛应用于勘探石油、医药开发 、天气建模、数据备份等领域。当前大多数的分子动力学模拟研究基于经验 势，通常包括两体互作用势，多体互作用势，各 种经验力场等。模拟的物理系统在材料科学方面 包括各种金属、半导体、聚合物和陶瓷等。在生 物物理、生物化学以及药物研究方面，分子动力 学更是被广泛的应用来模拟蛋白质、DNA 和各种 药物。目前使用经验势的分子动力学模拟的系统 规模在一般的计算机上也可以达到数万个原子， 模拟的时间在纳秒数量级。分子动力学的发展方向可以概括为以下几 个方面：更高的精度、更快的速度、更大的模 拟系统和更长的模拟时间。 首先是在精度方面，通过与第一性原理， 以及半经验的量子力学方法相结合，可以得到 不错的精度。另一方面是发展更为复杂和精细 的力场模型，得到与实验更为吻合的经验势。在加快模拟速度方面，针对分子动力学的并 行计算方法的提出可以充分利用现代高速并行计 算机的计算能力，还有一些专门针对分子动力学 而设计的硬件加速卡，比如日本的GRAPE，以 及利用图形显示加速卡进行加速（GPU计算）。 目前GPU计算十分流行，对于一些简单的分 子动力学计算，GPU的计算速度是CPU的几十到 一百倍左右。不少分子动力学软件已经直接支持 GPU计算。在扩大模拟系统的规模和延长模拟的时间方 面，许多针对大规模和长时间模拟的近似方法被 提出，比如多尺度方法将连续体模型与分子动力 学结合加速模拟大体系的演化，温度加速算法则 通过分子动力学与动力学蒙特卡罗的结合来加速 在低温和常温下的一些稀有事件模拟。MD的基本要素和流程 统计物理学量、边界条件、互作用势函数、 动力学积分和热力学控制构成分子动力学模 拟的基本要素。 分子动力学流程：最小影像约定 对于宏观材料，具有摩尔量级的原子数， 要模拟这么多原子是不现实的。如果不是特别 要去研究样品的表面，那么可以采用周期性边 界条件，如下图所示：此时，模拟的原子被认为是放在一个盒子 里，一个盒子里的坐标为r的原子用来表示位 于原子r’的一系列影像原子： 盒子里的原子i不仅仅和盒子里的另一个原子j有 互作用，而且和其它所有的影响原子i’,j’ …… 都有互作用。利用周期性边界条件，我们 消除了表面的影响，但这意味着我们要计算无穷 多原子之间的相互作用。所以在实际模拟中，我们都是对是势函数做 截断，这样只要考虑有限原子间的互作用了。 这种截断对分子动力学模拟的编程实现非常重 要，也是现代大规模并行模拟的基础 当我们选取的盒子足够大的时候，还可以引 入一个更为广泛采用的简化计算—最小影像约 定。即如果互作用随着距离衰减很快，盒子中 的第i个原子和j的距离Rij可以取为第i个原子与 第j个原子以及其所有影像原子（j’,j’’）之间的 最小距离：rij为真实的空间距离，Rij为采用周期性边界 条件后分子动力学程序中实际使用的距离。 需要注意的是，当我们使用最小影像约定 的时候，我们已经默认了这样一个条件：我们 的盒子足够的大，使得一个原子仅仅和自己所 在盒子或相邻的盒子里的一个影像原子才有互 作用，对于一个长方体形状的盒子，即要求最 短的一条边大于截断距离的两倍：非周期性边界条件 非周期性边界条件其实是最自然的处理方 式，随着纳米技术的兴起，在纳米尺度下材料 所含有的原子数不再是巨大的天文数字（已经 完全在现代计算机可以处理的范围之内），而 且纳米尺度下的表面和纬度效应非常重要，对 于低维纳米材料，一般采用非周期边界条件， 至少某一维度上采用非周期性边界条件。互作用势函数 分子动力学模拟中每个原子受到的力主要 由原子间的互作用势函数U来决定，势能的梯 度即为受力。所以选取合适的势函数直接关系 到分子动力学模拟结果的可靠性。 第一性原理计算可以比较精确地获得原子 间的互作用力，但由于计算量的问题，第一性 原理分子动力学所能处理的体系不可能太大， 模拟时间也不能很长。因此我们需要构造一些 经验的势函数形式。用经验参数拟合的势函数 具有简单高效的特点，被广泛应用于分子动力 学模拟。LJ势主要用于描述Van der Waals类型的相 互作用，最常用的LJ势是12-6形式： Lennard-Jones势 一般6次项用来描述没有化学键连接的两个 原子的吸引势，而12次项模拟原子间满壳层电子 的排斥力。LJ势是一个两体势，形式非常简单。LJ势能曲线多体势--金属：嵌入原子势(EAM) 多体势于20 世纪80 年代开始出现。其基本 的物理思想是把电子分为两部分。一部分为 内层满壳层电子，这部分电子的作用在势函 数上表现为两两之间的排斥势；另一部分则 为价电子部分，这部分和原子类型以及周围 的成键环境密切相关，正是这部分体现出晶 体中原子的“多体”互作用。嵌入原子势（EAM）是第一个被广泛采用 的经验多体势函数，用来描述含d 电子的金 属。EAM 势函数形式可以通过紧束缚模型 的低阶近似而得到，具有如下形式除了EAM势之外，还有很多类似的多体势函 数，比如Finnis-Sinclair势，Glue势，有效介质 势等等。多体势—半导体：Stillinger-Weber势 含s-p 电子的共价键的体系（比如半导体） ，对经验多体势函数是更大挑战。因为共价 键具有比金属键强得多的方向性。 SW势是第一个针对半导体的经验多体势函 数，包含了两体和三体项：Tersoff-Brenner势 Tersoff在1988年根据紧束缚的键序理论，提 出了对第四主族半导体（C,Si,Ge)的著名的 Tersoff势：其中fR和fA分别代表排斥和吸引互作用，bij 反映出Tersoff势的多体性。g(θ)反应了键角 改变的势函数。fC(r)是一个具有截断距离的 衰减函数，R+D就是截断距离 后来Brenner将其形式简化并推广到碳氢系 统中，得到了广泛的应用，通常称为 Tersoff-Brenner势。力场 在分子动力学模拟领域，力场是一大类被广 泛采用的经验势，特别是在有机分子和生物 分子的模拟中有着极为广泛应用。力场模型 是通过化学键的途径，把势函数分为若干项 。通常的力场包含：非键互作用（Van der Waals 作用、静电互作用）、键的伸缩能、 键角弯曲能、二面角扭曲能、交叉项能量（ 上述互作用之间的耦合能）等等。力场通常是通过大量的实验数据拟合出来的 。常用的力场可以分为三大类：传统力场、 第二代力场和通用力场。这里我们仅仅给出 简单介绍。传统力场 AMBER力场：广泛使用，适合生物大分子。 CHARMM力场：对小分子体系到溶解化的大 分子体系都有很好的拟合。 MM2/MM3力场：使用非常广泛，主要针对有 机小分子。 CVFF力场：可用于无机体系的力场。第二代力场 CFF力场：包含CFF91,PCFF等，适用于很 多体系，包括有机小分子，生物大分子等。 COMPASS力场：Accelery公司开发的，适 合高分子体系。 MMF94力场：目前最准确的力场模型之一通用力场 ESFF力场：Accelery公司开发，适用于多种 有机和无机分子。 UFF力场：可以适用于几乎所有元素 Dreiding力场：适用于有机分子、主族元素动力学积分 我们如果选定了边界条件，也知道了原子 间相互作用势函数，那么接下来就需要求解牛 顿方程组： 该方程为二阶常微分方程组，可以通过对 时间的有限差分法求解，即将连续的时间离散 化为一个时间序列，离散的时间间隔为Δt。时间步长 在实际分子动力学中，时间不长的选取需要 丰富的经验。步长过大，则数值误差对结果 又明显影响，而步长太小，又会浪费很多时 间。 步长的选取和积分方法有关，但总的来说取 决于我们所模拟的系统。一般来说，系统中 最高振动频率越高，则时间步长应该取的越 小；反之则可以取大些。Verlet算法 t时刻的加速度a(t)由牛顿方程决定。Verlet算 法具有Δt4级别的误差，不显含速度变量。速度 可以由下面公式求得，具有Δt2的精度：Verlet算法的优点是执行明确简单，存储要 求适度。在实际中很少有长时间的能量漂移 现象，也就是说具有很好的稳定性。因此得 到广泛的应用。 缺点是不显含速度，要求速度只有Δt2的精 度，低于位置的精度。Leap-Frog算法 Leap-frog算法可以从Verlet算法导出，它在 半个积分时间步的时候得到速度，并利用 这个速度计算新的位置： 该算法的优点是显含速度项。缺点是计算 的位置和速度不是同步的，不能同时计算 势能和动能。Velocity-Verlet算法 速度Verlet算法是Verlet算法基础上的一个显 含速度的形式。它和Verlet算法具有相同的 计算量，又不损失位置和速度项的精度。该算法中，先要计算新的位置和力之后，才 可以计算新的速度。 速度Verlet算法同时给出位置和速度，而且 计算量和Verlet算法相同，但速度计算精度 更高。 速度Verlet算法具有极好的数值稳定性。自 该算法提出后，被广泛采用，目前已经成为 分子动力学运动方程积分的标准积分方法。Beeman算法 Beeman算法在位置上其实和Verlet算法是一 样的，但速度具有更高的精度。但表达式比 Verlet算法复杂，计算量也相应较大。Rugge-Kutta算法 Rugge-Kutta算法是工程和科学计算中应用 最为广泛的数值求解常微分方程组的方法。 该算法有很多不同精度的格式，也有一些可 以估算积分误差和自适应步长的方法，可以 根据不同的需要选取。 比较常见的是四-五阶和七-八阶的格式。 Rugge-Kutta方法具有很好的数值稳定性。 具有很高的短时积分精度，缺点是计算量和 存储量都偏大。统计物理学量 不同于蒙特卡罗方法，分子动力学通过求解 牛顿方程来得到系统在相空间的动力学演化 过程，因此是一种确定性的模拟方法，即当 给定了一个初始位型之后，分子动力学将模 拟真实的分子运动轨迹。分子动力学模拟系 统在相空间的演化过程给我们提供了关于系 统的大量样本，我们可以用统计物理学的方 法来得到系统的相关性质。对于一个时间无限长的模拟过程，我们可 以认为系统的相空间是被完全抽样了的。在这 种极限情况下，我们可以使用时间平均来代替 系综平均，从而得到系统的热力学量。这样通 过统计物理学的应用，我们可以用分子动力学 模拟来“测量”系统的各种热力学性质，当然 也可以用来模拟热力学相变，比如固液相变和 液气相变等。在实际的模拟过程中，由于模拟 时间的限制，系统可能并没有达到平衡，所以 在模拟时要注意确保模拟的时间足够长。系统的一个物理量A的系综统计平均值<A>可 以用时间统计平均值来代替： 其方差的涨落为：势能、动能和温度 势能U即势函数U的值 动能K和经典统计物理学中一样： 系统的总能为两者之和 E=U+K 根据能量均分定理定义一个动力学温度T Nf为系统总的自由度。对于纳米尺度下，少数粒 子系统的温度定义是一个尚有争议的问题。但由 于这么计算T比较简单，仍被广泛使用。均方位移(MSD) 均方位移是时间的函数，< >是对所有原子的 平均。一般地，对平衡态体系，MSD正比于 时间t。扩散系数D 扩散系数可以由均方位移求的，其中d是 系统的维度（三维d=3）。通常对于固体，扩 散系数D是一个非常小趋近于0的值。而对于 液体和气体则是一个与温度有关的有限值。径向分布函数 粒子分布的对关联函数g(r)定义为： 径向分布函数为： 径向分布函数意义为给定某个粒子的坐标， 在距离这个粒子为r的地方找到其它粒子的概 率。可以用来描述物质的有序性。关联函数与响应函数及振动谱 如果A=B，就是自关联函数。上面提到的对关联 函数g(r)就是粒子密度在空间上的关联。关联函数与热力学响应函数密切相关，比如热 传导率与热流密度的自关联函数的Green Kubo 关系： 一个比较常用的关联函数就是速度在时间上的 自关联函数（N为粒子数）：C(t)的傅立叶变换就是振动谱（功率谱）：模拟系综以及热力学控制 实际的分子动力学模拟总是在设定的外界 条件下进行的，需要进行热力学控制，即需要 和外界进行热交换及力学平衡。平衡态的模拟 ，通常可以用不同的系综进行分类。分子动力 学中常用的有微正则（NVE）系综、正则（ NVT）系综、等温等压（NPT）系综。NVE 系综是分子动力学的自然系综，只需 要调整系统的能量就可以比较容易地实现。对 分子动力学而言则需要对温度T 进行调节，这 一般是通过热浴来实现。NPT 系综则需要对温 度T 进行调节的同时对体积V 进行调节以保证 压强P 不变，一般是通过热浴和压力浴来实现 的。 温度和压强是非常重要的热力学参量，它 们的控制也是分子动力学模拟中需要特别注意 的地方。速度标度热浴 改变温度最简单和最直观的方法是直接对 速度进行标度。即如果t 时刻系统的温度为T(t) ，而我们期望的温度是T0，那么我们可以通过 改变速度v 来调节系统的温度。最简单的标度方 法是每一步都对速度乘上一个标度因子：Anderson热浴 Anderson在1980年提出了通过随机选取的粒子 速度进行速度标度来控制温度： 其中n是随机选取的粒子标号，vr是满足温度为 T0的Maxwell-Boltzmann分布的随机速度。Berendsen热浴 Berendsen在1984年提出通过一个温度为 T0的虚拟热浴与系统接触来控制温度。该热浴 每一步对速度都按照速度标度方法进行标度， 但标度因子与热浴和系统的温度差有关：这里的参数τ是驰豫时间，具有时间量纲， 反应了系统和热浴的耦合强度。τ越小则两者耦 合越强。 一般模拟中τ应取得比系统中最高 频率周期要长，通常可以取为最高频率周期的10 倍。这是因为原子振动方式的热传导过程中，一 次能量交换可以认为至少要通过一个高频振动周 期才能完成。 Brendsen热浴允许系统温度在期望温度上下 波动。Nose-Hoover热浴 Nose-Hoover热浴通过引入一个热浴的自由 度ζ到整个系统的运动方程之中来模拟热浴温 度对一个具体粒子的运动的影响。这个自由度 体现为“摩擦力”的形式，与系统和热浴的温 度差有关。此时粒子的运动方程为：这里ζ具有阻尼系数的含义，而τ和 Berendsen热浴中一样是驰豫时间。Nose-Hoover 热浴也允许系统温度波动，而且没有速度标度 带来对运动积分的不连续性影响，是一种比较 缓和的热浴方法，因而对扩散和输运系数等的 计算影响很小。压力裕 Anderson压力浴 Berendsen压力浴 Parrinello-Rahman压力浴分子动力学程序 AMBER CPMD(DFT) CHARMM NAMD GROMACS GROMOS GULP LAMMPS DL_POLY TINKER …………石墨片在温度下的起伏 原子团簇B8 原子团簇B10模特卡罗模拟 Monte Carlo Simulation蒙特卡洛 模特卡罗是著名的赌城，而模特卡罗模拟是 一种基于“随机”的计算方法。这一方法源 于美国在第一次世界大战进研制原子弹的“ 曼哈顿计划”。但模特卡罗的思想很早被 人们发现和利用。Buffon投针实验(1777) 假设地板上画着一组间距为1的平行线。把 一根长度为1的针扔到地上。设N次抛针，M次 相交，则相交概率的期望值为2/pi，即： 2N M 历史上一些实际的投针实验：模特卡罗求π 将n根飞镖随机投向一正方形的 靶子，计算落入此正方形的内切 圆中的飞镖数目k。假定飞镖击 中方形靶子任一点的概率相等。 2 244 kr nr 那么落入内切圆的概率为： 通过计算机随机撒点，可以求出π： n k4计算机算法 为简单起见，只以上图的右上1/4象限为样本 k ← 0; for i ← 1 to n do { x ← uniform(0, 1); y ← uniform(0, 1); // 随机产生点(x,y) if (x2 + y2 ≤ 1) then k++; //圆内 } return 4k/n; 实验结果： π=3.141592654 n = 1000万: 3.140740, 3.142568 (2位精确) n = 1亿: 3.141691, 3.141363 (3位精确) n = 10亿: 3.141527, 3.141507 (4位精确)计算定积分 1 x y 1 设f: [0, 1] → [0, 1]是一个连续函数， 则由曲线y=f(x), x轴, y轴和直线x=1 围成的面积由下述积分给出： 1 0 ()S f x dx  向单位面积的正方形内投镖n次，落入阴影部分的镖 的数目为k，则： /1 kS S k nn   模特卡罗方法计算积分的方法比较简单，程 序实现也很简单。模特卡罗特别适合于计算 高维的积分，比如4维以上的定积分。 模特卡罗方法在许多领域都有着广泛的应用 ：比如统计物理、核物理、金融。