p(可=mA.(r,)d (2.2-7) 注意，写出此式时我们已隐含假定t是足够长的，该时均值与初始 条件无关。实际上，一般来说这是一个不真实的微妙很设。然而我 们将不顾此种微妙，而直接假定，一旦指定N,V及E,各个时 均值与初始坐标和动量无关。这样的话如果对许多不同的初始状态 求平均，将不会改变P:(r)的结果。也就是说我们考察的是假想状 态，在此状态下我们在同一N,V及E下，具有不同的初始坐标 和动量进行大量分子动力学模拟 JA(,ro,poc】 超始赛件 初始条件数 (2.2-8) 现在来考察极限情况，即对所有与赋值N,V及E相容的初 始条件进行平均，则可用积分替代射所有初始条件的加和 ∑fiN(0),p(0)] 初蟑条修 初始条件数 f[rN(0),pN(0)]drNdpN 可(N,V,E) (2.2-9) 式中f表示初始坐标(0)，p(0)的任一函数，且 n(N,V,E)=drNdpN(已忽略一定值系数●)。积分符上的 下标E表明积分局限在能量为恒定值E的层面上。此种“相空间” 积分对应了在前一部分所讨论的系综平均的经典极限●。本文用 ●如考素一量子为学体系，则(N,V,E)表示对一给定N,V及E的量子 态数目。在经真的范倒中，具有N个可分辨，无结构粒子的维体系的量子杰数目为 口(N,V,E)。对于N个不可分辨的粒子，应将后者再除以系数NI。 ●此妞我们考察一个轻典等同的所谓撒正则系综，即具有固定N,V及E的体 系的系综。在墩正则系综中相空闯积分的经典表达式可以由一种量子力学表达式导出， 在后者中所含的各量子态的求和很大程度上与由量子力学表达式推导经典的恒定N, V,T(正则)系综的方法相间。 16
<…>表示系综平均以便与由一横线表示的时间平均相区别。如果 将时间平均与对所有不同的初始态的平均的次序切换，则可得 pn可=▣0a(r0,poEr2.2-10) 然而在上式中的系综平均与时间t′无关。之所以如此是因为 在体系的初始相空间坐标及其在随后时刻t的坐标之间存在一一 对应关系[41,2]。因此对所有初始相空间的坐标的平均值等同于对 随时间演化的相空间坐标求平均值。由此可以省去式(2.2-10)中 的时间平均，并得到 P:(r)=〈a(r)》NwE (2.2-11) 此式表明如果想计算为粒子体系的坐标和动量的函数平均值，可以 计算其时间值(MD法)或系综平均(MC法)。应予指出上述各 节只是意味着式（2.2-11)是可接受的，而非一种证明。实际上由 于式(2.2-11)一般来说并不正确，要想证明看来似乎是不太可能 的。然而，在下文中将仅假定“各态历经假说”[如式（2.2-11) 所示]适用于在计算机模拟中的所有体系。不过读者应明白体系中 的各个举例，实际上并非各态历经的，例如玻璃态和亚稳态，甚至 有的原理上也并非如此，如近似谐振固体。 17
3 Monte Carlo模拟 在本章中叙述Monte Carlo方法的基本原理，尤其着重于固定 粒子数N的给定体积V及温度T下的体系 3.1 Monte Carlo方法 在上一章中已介绍了（经典)统计力学的基本概念。下一步目 标是指明Monte Carlo方法的切人点。我们从配分函数的经典表达 式（2.2-5)开始 Q=cexp[-H(rV,pN)/E8T]dpNdrN (3.1-1) 式中表示全部N个粒子的坐标，为对应的动量。H (,P)是体系的哈花尔顿函数。.它表示一个孤立体系的总能量 是组分拉子的坐标和动量的函数：H=K+U,其中K是体系的 动能，U为势能，最后，c是一比例常数，其选择应使式(2.1- 14)中的对应量子态的加和在-·0的极限时趋向经典配分函数。 例如对于有N个相同原子的体系，c=1/(五3N!)。对应于式 (2.2-1)的经典公式为 A(pN,)expl-BH(pN,rN)]dpdrN (A)= (3.1-2) exp[-H(pN,rN)]dpNdrN 式中3=1/T。在此式中，已将可观测量A看成坐标及动 量的函数。由于K是动量的二次函数，可以对动量得到解析积分。 因此与动重有关的函数的平均值，往往易于求得。难题是计算各种 函数的平均值。只有在一些例外的情祝中，对粒子坐标的多维积分 18
才能解析计算，在其他的情形中，必须采用数值方法。 至此已确定了必须解的数值问题的性质，接着来看可能解决的 方案。看来最直接的途径或许是由数值面积积分，如用辛普森法 则，计算式(3.1-2)中的(A》。然而很容易看出，即使是独立坐 标数dN(d是体系的维数)很小的O(100),这种算法也是完全 不适用的。假定我们打算用在dN维构型空间的网点上的被积函数 求积。如沿各个坐标轴取m个等距离点，则总的必须计算的被积 函数的点数为mw。对于除最小体系之外的所有的体系，此数目 会变成天文数字一般大。例如，在二维空间中的100个粒子，且 =5,则将在1010个点上计算被积函数。这种数量级的计算，在 已知的世界上是无法进行的。由于所得到的答案可能受到大的统计 误差的影响，这种计算量还算是值得庆幸的。当然，在对应于网目 尺寸的距离上各个函数是平滑的情况下，数值求积是最适用的。不 过对于绝大多数分子的势能，式(3.1-2)中的玻尔兹曼因子随粒 子坐标而急速变化，因此准确的求积需要小的网格距离（即大的 m值)，此外，当计算稠密液体的被积函数时，对占绝大多数的点 来说，玻尔兹曼因子的值小到可忽略。例如对于处于凝固点的100 个硬球流体而言，每100个构型中仅有一个值不为零。 上述例子清楚表明，需要一种较好的数值方法来计算热平均 值。一种方法是Monte Carlo法，或更为准确地说是由Metropolis 等[}于1953年提出的Monte Carlo重要性抽样方法。此方法在稠密 分子体系的数值模拟中的应用是本章的主题。 3.1.1重要性抽样 在讨论重要性抽样之前，首先考察最简单的Monte Carlo方法， 即随机抽样。假定要数值估算一维积分1 I =f(z)dz (3.1-3) 代之以常用的面积积分，即在横坐标的预定值处计算被积函数，我 们可另辟新径。注意到式(3.1-3)可另写为 I=(b-a)(f(x)> (3.1-4) 19
式中〈f(x)〉表示f(x)在区间[a,b]上的无权重平均 值。在强制性的Monte Carlo中，此平均值是由在区间[a,b]中 于大量（如数目为L)的随机分布的x值处计算f(x)来确定 的。显然，当L→∞，此方法一定可得到正确的1值。然而如同 采用常用的面积积分法一样，在计算如式（3.1-2)所给出的平均 值时，此种方法是无所作为的，因为大多数计算花费在玻尔兹曼因 子可忽略不计的点上。显然很倾向子在玻尔兹曼因子大的许多点上 抽样，在其他地方则少抽样。这就是重要性抽样的基本观点。 应该如何在整个构型空间分配抽样？要了解这一点，让我们先 考虑一个简单的一维例子。假定想用Monte Carlo抽样来计算式 (3.1-3)的定积分，按某种非负的概率密度边(x）在区间[a, b](为了便利，假定a÷0及b=1)非均匀分布点上抽样。于是， 可将式(3.1-3)改写为 (3.1-5) 假定已知w(x)是另一个茵数u(x(非负，非递诚)的导数，且u(0)=0 及u(1)=1(此边界条件意味着(x)是整过的)。于是I可写成 1-} (3.1-6) 在式(3.1-6)中x(w)表明如果将4作为积分变量，则x必定 为u的函数。下一步即为在区间[0,1]产生L个非均匀分布的 植。子是可得到！的估叶值 i≈是为 f(u)l w[x(u)] (3.1-7) 按此方式写出I,我们获得了什么？答案主要取决于所选择的 w(x)。为了说明这一点，估算【红中的方差，其中表示采 用L个随机抽样点由式(3.1-7)所获得的1的估算值 i会2胡-]别 (3.1-8) 式中角括号表示真值，也即是当L→∞的极限时所应取的值。 20