行。这与如量子力学编码相比就非常短少。因此一个模拟系统往往 写出许多适用于特定用途的不同程序。其结果是没有一种标准的 Monte Carlo或分子动力学模拟程序。然而绝大多数MC或MD程 序的核心部分（如不完全相同）是很相似的。下面将构筑这样一个 核心，它将是初步的，以清晰代之以效率。但它将演示Mo咖te Carlo 方法是如何工作的。 3.2.1算法 我们将讨论的这类Monte Carlo或分子动力学程序的主要目的 是计算经典多体体系的平衡物性。此后，我们将此类程序称之为 MC或MD程序，尽管Monte Carlo方法有许多其他的应用。现在 让我们考察一简单的Monte Carlo程序。 在前一节中Metropolis方法是以一种马尔科夫过程引入的，在 此方法中随机行走是按访问一特殊点的概率正比于玻尔兹曼因 子exp[一U(r)打]来构筑的。有许多方式来构筑随机行走。在 Metropolis等[6]所介绍的方法中，提出了下述步臻： 1)随机选择一个粒乎，并计算其能量U(r)。 2)给该粒子一随机位移，r'=r+A,并计算其新能量。 3)按下式的概率接麦从至rN的移动 P∝(o+n)=min{1,exp[-[U(r)-U(rN)](3.2-1) 此基本Metropolis方法示于算法3-1及算法3-2中。 算法3-1盖本Mietrepolis算法 RROGRAM me 基本的Metropolis算法 do icycl =1,ncyel 执行ncycl次MC循环 call mcmove 移动一个分子 if (mod (icycl,nsamp).eq.0) call sample 样本平均 enddo end 算法注解 26
1)子程序ncmove。对被选中的分子进行尝试移动（见算法3-2）。 2)子程序sample.。每nsamp次循环进行一次抽样。 算法32移动分子娄试 SUBROUTINE mcmove 移动分子尝试 o=int (ranf ()npart)+1 随机选取一分子 call ener(x(o),eno） 老构型的能量 xn=x (o)+(ranf ()-0.5)*delx 给分子一个随机移动 call ener (xn.eno) 新构型的能量 if(ranf().lt.exp(-beta￥(enn-eno)） 接受准则(3.21) x (o)=xn 接受：用xn替换X（o) return end 算法注解 1)子程序ener计算给定位置分子的能量。 2)注意，如果构型被拒受，则保持原构型。 3)ranf()是均匀分布在[0,1]上的随机数， 3.2.2技术细节 在本节中我们讨论许多在设计一个有效率的模拟程序时有重要 实际意义的计算技巧。应予强调的是绝大多数技巧虽然是十分有用 的，但却非独有的并且缺乏深刻的物理含义。然而这并不意味着使 用这些计算工具既无风险又不乏细徽差别。理论上，节省时间的一 些方法并不一定会影响摸拟的结果。但在一些情形中，节省时间的 技巧对模拟的结果确有可观的影响。对于许多不同的避免直接计算 相距甚远的粒子间的分子相互作用的方法，这一点尤为真切。此 外，我们通常可以估算这种省时方法的所不希望的副作用，并对其 加以校正。 (1)边界条件 原子和分子体系的Monte Carlo及分子动力学的目的是提供一 个宏观样本的物性信息。诚然，在当今的计算机中通常还只能处理 几百至几百万这样数目的自由度。绝大多数模拟检测几百至儿千个 27
粒子的体系的结构及热力学物性。显然这个数目与热力学极限相去 甚远。更为精确地说，对于这类体系不能假定可靠的边界条件（例 如自由的或刚性的，或周期的)对体系的物性影响甚微。实际上在 具有自由边界的三维N个粒子的体系中处于界面的分子数正比于 N1B。例如在一个具有1000个原子的简单立方晶体中，约49% 的原子位于界面，若是10个原子，此分率减少至仅为6%。 为了模拟主体相，选择能反映围绕着我们的N个粒子的无限 的主体相是至关紧要的。这往往可以借助于使用周期性边界条件来 实现。将含有N个粒子的体积当作是具有与其相同的单元的无限 周期点阵的原始单元（见图3-2)，一给定粒子（臂如）则与在此 无限周期体系中的全部其他粒子相互作用。也就是说，全部其他粒 子在完全相同的单元中，而且全部粒子（包括它自己的周期性映 像)在全部其他单元中。举例来说，如假定所有的分子间相互作用 是成对加和的，则N个粒子在任何一个周期盒子中的总势能为 Um=受u1+L) 式中L是周期盒子（为便利假定为立方的）的直径，是三 个整数的任一向量，面加和上的一撒表明当=0时，=j的项应 予排除。在此非常一般化的式子中，周期性边界条件并不十分有 图3-2周期性边界条件的图示表征 28
用，因为若要模拟主体性质，必须将有限加和改写成无限加和●。 实际上我们大多数处理短程相互作用，此时往往允许对所有的分子 间相互作用在超出一定的截断距离。进行裁断。在实践中如何做 到此点将随后讨论。 虽然已表明边界条件是模拟均相主体体系的特别有效的方法， 但人们还应当明白这种边界条件的应用会导致并不存在于真实的宏 观主体体系中的虚假的相关性。模型体系的周期性的后果是，只有 那些与周期点阵波长相容的起伏是允许的，即与周期盒子相适应的 最长的波长为A=L。如果预想波长涨落是重要的（如在连续相变 附近)，人们则应该事先考虑采用周期性边界条件导致的各种问题。 另一种由于采用周期性边界条件的非物理影响是已发现稠密原子流 体的径向分布函数多()不是各向同性的。 最后，应予指出周期性边界条件的常见的错误概念，即周期性 盒子的边界有某种特殊含义。其实不然。·原始单元的周期性点阵的 原点可选在任何地方，·这种选择不会影响所考察的模型体系的任何 物性。相反，应该固定的是周期单元的形状和方位。 (2)相互作用的截断 现考虑执行具有短程相互作用的体系的模拟。短程意味着对于 一给定的粒子i的全部势能是由较某一截断距离r。为近的相邻粒 子之间的相互作用所决定。用选择足够大的ǐ。的方法，可以使得 由于忽略了与较大距离的粒子的相互作用而造成的误差任意小。如 果我们采用周期性边界条件，。小于L2(周期盒子的直径的一 半)具有特殊的意义，因为此时需要考虑的是给定的粒子：仅与任 一其他粒于j的最近的周期映像的相互作用（见图3-2中的虚线盒 子)。若分子间作用势能在r≥r。时不严格为零，分子间相互作用 的截断会引起U的系统误差。假定分子间相互作用迅速衰减，便 可以对U增加一尾部贡献来校正此系统误差 鲁事实上，在最初的Lennard-.Jones粒子三维MC摸拟中，Wood和Parker4讨论 ‘了与此处讨论的常用方法有关的这种无限加和。 29
0w=是(r)+当2u(,)4xr2dr (3.2-2) 式中。表示截断后的势能函数，P为平均数密度。在写这一 表达式时，隐含假设当r≥r。时径向分布函数g(r)=1。显然， 只有当尾部校正小时，最近邻的周期映像规定才适用。由式(3.2 2)还可看出，只有当势能函数u(r)比r3(在三维方向)衰减 得更快时，对势能的尾部校正才是无限小。如分子间的长程作用以 色散力为主，则满足此条件。然而对于非常重要的库仑力及偶极相 互作用情形，尾部校正发散，因此最近邻映像规定不能用于这类体 系。此时必须直接考虑与所有的周期映像的相互作用。附录B1中 表述了处理此作用的几种方法。 几种因素使得势能截断成为一件具有技巧性的工作。首先必须 认识到虽然势能函数的绝对值嘘粒子间的距离而减小，但对于 足够大的x,近邻的粒子数目随而快速增加。实际上在一给定原 子距离为y的粒子数目按4~1逐渐增大，此处d是体系的维数。 例如，让我们计算简单的三维Lennard-Jones流体的成对势能的截 断效应，这一较通用的成对势能可表示为 (r).=4[(）严-(] (3.2-3) 任一给定原子：的平均势能（在三维方向）为 =re(r)u(r)dr 式中(r)表示离一给定原子i距离为r的平均数密度。所含的 系数12校正了分子间的相互作用的双重计数。如在距离r。处截断 势能，则忽咯的尾部势舶贡忧1为 （）4r2o()(,dr (3.2-4) 式中由于体系的全部原子都相同，我们省去了下标i。为简化 ua的计算，假定r≥r。时，密度p(r）等于平均数密度p。如果 u(r)为Lennard-Jones势能，则uai为 30