计算物理方法(Sec3).ppt

1、第三部分,Monte Carlo 方法,第一节 用Monte Carlo 方法模拟凝聚态物理系统的基本思想,凝聚态系统的随机特性：1. 粒子系统量子态的量子随机性2. 大量粒子的热统计随机性- 适合用计算机模拟,第二节 随机游动及应用,随机游动是一种基于运用 0，1 区间的均匀分布随机数序列来进行的计算。早在1906年Pearson就提了“随机游动”的问题。以后随着其理论的逐步完善，随机游动模型在物理学、生物学和社会科学中都得到广泛的应用。许多教科书中都可以找到它在诸如气体分子扩散、液体中悬浮物的布朗运动、量子力学中薛定锷方程的求解、高分子长链的特件研究、求解偏微分方程和数学积分的近似计算等中

2、的成功应用。我们在介绍它的应用之前，有必要首先介绍一下随机游动模型。,醉汉的一维行走问题初始： 电杆位置 x=0，步长 l ，每一步的取向是随机的，右行几率为 p，左行几率为q=1-p。问题：醉汉在行走N步以后，离电杆的距离为x的概率 。有了 后，可以计算：,可用概率理论解析地分析,虽然这里用了很简单的解析方法得到上式，但是一般情况下，能精确求解游动问题的技术却不是这样简单。有两种重要的方法可以用于游动问题，它们是查点法和蒙特卡洛 方法。查点法：对给定的行走总步数N及总位移x，要求把游动时可能的每一步的坐标和几率都确定下来。,从上面的分析可以看出：查点法只有在总步数较小时才可以应用，N比较大时

3、用起来就比较困难了。对比查点法，蒙特卡洛方法就可以克服在游动中的这个困难，具有可操作性。蒙特卡洛方法可以对许多步的游动过程进行抽样。我们以随机游动的蒙特卡洛方法在求解泊松型微分方程中的应用作为例子。若该泊松方程及其边界条件为,为求解区域D的边界，s为边界 上的点。这里我们采用等步长h的正方形格点划分的差分法。在区域D内的任意正则内点o(其相邻的节点都在区域D内)的函数值可以用周围四个邻近点1，2，3，4上的函数值来表示。这个表达式有如下差分方程表示,前面所述类型的随机游动或链具有如下特征；它在游走中任一阶段的行为都不被先前游动过程的历史所限制，即区域内的点可以被多次访问，这种随机游动过程叫做马

4、尔科夫过程。又因为游动最终会终止在边界上，故而上述的这类游动也称为马尔科夫链。马尔科夫涟正是这样生成相继各状态的，它使得后一个状态在前一个状态的邻近。由此可以知道相继各状态之间的确存在着关联。马尔科夫链是分子动力学中由运动方程生成的轨道在概率方面的对应物。对统计力学系统进行蒙特卡洛模拟计算将在本章第4节中介绍。另外还有一种非马尔科夫过程。自规避随机游动过程就是属于这一类。在这个过程中任何一步的游动概率都要考虑前面游动的历史，因而游动将有可能在碰到边界前就被强行终止掉。随机游动对一些更抽象的问题也是非常有用的。,两个重要的概念,随机行走的概念-行走方向的几率-按该几率实现行走 路径平均、路径求和

5、、路径积分的概念-某一路径给出所求物理量的一个值-不同路径给出不同值-对所有路径给出的值求平均,第三节 量子Monte Carlo 方法,量子力学中的波函数是直接与几率密度相关的量、我们有分布密度函数的关系式,其中， 为归一化常数，因此波函数也被称为几率幅度。人们很自然地想到可以利用蒙特卡洛方法来求解量子力学问题。用于求解量子系统的薛定諤方程的蒙特卡洛模拟方法通称为量子蒙持卡洛方法。在实际应用中主要有变分蒙特卡洛方法(VMC)，路径积分蒙特卡洛方法(PIMC)和格林函数蒙特卡洛方法(GFMC)等。在本节我们将介绍路径积分量子蒙特卡洛方法和变分量子蒙特卡洛方法，并对格林函数Monte Carlo

6、方法作入门的了解。,一、量子力学回顾 量子力学的基本方程是薛定锷方程：,为微观体系的哈密顿算符。对微观粒子，其哈密顿算符可以写为,为势函数算符。求解哈密顿算符所对应的能量本征态的波函数和能量本征值是量子力学的基本内容。若知道初始态的波函数为 ，波动方程则有唯一的波函数解及以后时刻的几率密度 。从费曼的观点来看，一个粒子,在某个时刻 t ，某空间位置 x 的波函数应当是来自所有的初始态位置“传播”到该时空点的幅度。即,上式中的 称为“传播子”。它表示在初姑时刻 ，空间位置 点的波函数值对下一时刻 ，在x点上的波函数值的贡献强度。该传播子可以表示为,如果 为与时间无关的哈密顿算符的本征态波函数，它

7、满足的薛定锷方程为,于是：,其中：,于是：,假定该等式在延拓到t为虚值时仍成立，令 ，则有,当 足够大时，特别是在 时( 是基态能量， 为第一激发态的能量)，上式的右边主要是来自能量最小的基态能量 的贡献。如果我们取 并忽略其他的贡献项，则有,于是：,利用归一化的要求，基态波函数绝对值的平方可用传播子表示为,我们现在必须计算传播子。将 时间间隔分为 个等时间间隔 的小区间，则此间隔为 ，并且 ，( ) 。根据完备座标表象的关系式,有：,取连续极限得到,其中常数 ， 为沿路径的经典作用量。上式表示传播子是由连接初态 和末态 的所有路径，通过相因子 所做的贡献。其中 是系统的拉氏量。上式中 是所有

8、各种可能的分 段直线段构成的路径 ( )之和的总作用量。同样，如果我们假定将 延拓到虚数范围时，上述等式仍然成立。令 ，则上式中的作用量可以推出为,于是：,类似前面的推导，上式中指数中有一个路径积分，它的积分是沿路径 ，即我们把路径积分的空间起终点 和 分别放在 上，则该积分为,因而对于每一条路径，就有一个能量。上式于是有如下形式；,由于 ，并对 进行积分，此时须加进一个 函数在被积函数中，则上式可等价写为,上面的公式显示出量子力学中的费曼路径积分在欧氏时空的表示，揭示出量子理论与统计力学之间的深刻联系。这时的路径积分与配分函数两者在数学上是相同的，因而我们可以用计算经典统计力学配分函数的做法

9、来计算路径积分问题。,二、路径积分量子蒙特卡洛方法下面我们就用路径积分蒙特卡洛方法求解薛定锷方程的基态能量和基态波函数的数值。从上面两个公式可以使我们联想到玻尔兹曼分布。变量 的位形分布密度函数正好是将玻尔兹曼分布中的 换成 。 可以被视为函数 在位形 (每个位形对应一条路径) 在此分布下的平均位。其分布的数学表示为,这里存在的一个关键问题是：上面公式中给出的 具体形式计算起来并不方便。在计算归一化常数 时，包含了一个积分。这个计算实际上是一个高维的多重积分的计算。费曼路径积分量子化的欧氏积分表示中的积分计算也仍然主要是个蒙特卡洛计算问题，对它们的积分计算可以离散化为对路径的求和。但是采用一般

10、随机抽取位形点的办法，效率是很低的。尤其是在此高级空间中做均匀抽样时，由于指数项的缘故，大量的点会落到对求和贡献非常小的区域。 此时如果我们采用马尔科夫随机游动的重要抽样方法Metropolis方法，将是十分有效的。利用Metropolis方法，按照类似玻尔兹曼分布的分布函数来抽取若干位形 ，便可以计算出基态波函数 的估计值，然后对该估计值求平均便得到 的值。,这种方法在求解一维基态波函数时优越性并不明显。但是在更复杂的量子力学计算中，采用路径职分方法就显示出极大的优越性。这半要是由于在传统的场论计算中，势函数的作用是用在真空上的微扰方法来处理的；而在路径积分中，是将势函数插入到作用量积分中去

11、求数值解，事实上是在做精确计算的尝试。前一种方法对电弱作用的计算很有效，但对于有强相互作用的问题，其使用价值不大。在强相互作用中，矩阵元不能够以强耦合常数展开为收敛的级数。另一个优点是该方法将时空离散化为格点，这将带来数值计算上的方便。此外，采用Metropolis游走方法来选择具有代表性的态是非常有用的。该方法不仅可以以简洁的数组方式给出场的描述，还能够对积分加上截断，以保证在将格点上的离散时空延拓到连续时空时微扰理论的重整化。,三、Metropolis 随机抽样方法,在随机游动的Monte Carlo方法中，有一种最常用方法称为Metropolis方法，它是重要抽样法的一个特殊情况；采用此

12、方法可以产生任意分布的随机数，包拈无法归一化的分布密度函数。Metropolis方法是通过某种方式的”随机游动”来实现的。只要这个随机游动过程拉照一定规则来进行，那么在进行大量的游动，并达到平衡之后所产生的点的分布就满足所要求的分布 。Metropolis方法所采用的游动规则是选择一个从x 点游动到x 点的“过渡几率” ，使得它在游动中所走过的点 的分布收敛到系统达到平衡时的分布 。要达到这样的重要抽样的目的，就需要对过渡几率 的选择加上适当的限制：,对于相空间中点集的一切互补的对偶集 ，存在着 和 ，使得 。这是相空间区域的连通性和遍历性的陈述。 由于概率正定性的要求，对于一切 ，。 概率归

13、一化要求，对所有的 ， 。 由于要求极限分布为平衡分布，所以对所有的x， 可以证明，只要游动所选的”过渡几率”满足如下的细致平衡条件就可以达到平衡时分布为 这样的目的：,第五项中的细致平衡条件实际只是一个充分条件，并不是一个必要条件。该条件并不能唯一确定过渡几率 。所以，过渡几率 的选择有很大的自由度。不同的选取即是不同的方法；在Metropolis方法中一般采用一个简单的选择过渡几率的方法：,设游动到了 点，以下的具体操作步骤为：（1）首先选取一个试探位置，假定该点位置为：，其中 为在间隔 内均匀分布的随机数。（2）计算 的数值。（3）如果不等式 满足，那么进行这一游动，并取 ，返回（1）开

14、始对游动到 点的试探。（4）如果 ，产生一个 区间内的随机数 。,（5）如果此时 ，那么还接受这步游动，并取这步游动所到达的点为 ，然后返回到步骤（1），开始下一步到达 点的游动。（6）如果此时 ，就拒绝游动，仍留在 的位置不变。（7）返回到步骤（1），重新开始对游动到 点的又一次试探。,必须指出：采用这样的游动过程，只有在产生了大量的点 后才能得到收敛到满足分布 的集。这里有一个明显的重要问题，就是如何选择 的大小，才能提高游动的效率?如果 选得太大，那么绝大部分试探的步子都将会被舍弃，就很难达到平衡分布：反之，如果 取得太小，那么绝大部分试探步子都会被接受，这同样难以达到所要求的平衡分布。

15、根据实际应用的经验，选取 的一个粗略标准应当是：选择适当 大小的原则是要在游动的试探过程中，有13到12的试探步子将被接受。按照这样的标准选择得到的 ，就可以大大提高游动的效率。,另一个在Metropolis方法中的问题是：进行这样的随机游动，从哪一点出发才可以比较快地达到平衡分布呢?原则上讲，从任何一个初始位置出发均可达到平衡分布，但是为了尽快地达到平衡分布，我们最好是要选择一个合适的初始位置这个初始位置应当是在游动范围内所要求的几率分布密度 最大的区域。,四、例子-一维简谐振子假定有一个质量为m的粒子，在一维简单简谐势,中运动。取 为长度单位， 为时间 中的单位。有,首先，选择任意的、连接

16、 个时间间隔、且的一条路径，计算上式的能量，然后再选一系列路径，每条路径与前条路径最多只有在一个时刻（例如 ），有不相同的空间点（见图）。采用Metropolis方法来确定满足上面要求的新径迹。其中将随机定下的坐标 改变到 的过渡几率为 ， 为分别包括在 时刻坐标为 和 的两条径迹的能量差，可由上式算出。这样的随机游动抽样得到的径迹也许会与前一个径迹相同。每当新径迹选出后，就利用前式计算被积函数 的估计值，并累加到求和之中。最终该求和所得的值与抽样路径的总数相除所得到平均值就得到 的数值结果。按上述方法，游动足够多的步数后，找们就可以得到x点上 的值。,五、变分量子Monte Carlo方法通

17、过薛定谔方程求解基态本征能量和本征波函数,选择试探波函数 计算试探能量,其中， 可看作局域能量。,可看作空间点出现的几率。由哈密顿，有,采用Metropolis随机游动方法产生满足 分布的 个点 ，则,不断改变试探波函数的值，并计算试探能量的平均值，直到取得最小值，这时得到的试探波函数和能量平均值就是基态波函数和基态能量本征值。,下面我们以一个一维的量子体系的变分蒙特卡洛模拟步骤作为示范：（1）选择一个物理上合理的近似基态波函数 作为试探波函数；（2）采用Metropolis方法，按照分布密度函数 随机抽取 个点 用上述公式计算能量平均值 。（3）改变试探波函数的值，使得 的值在区间 内随机变

18、化一个小量即 ，重复（2）中能量平均值的计算得到 。,（4）计算能量平均值的改变值 ，如果 则接受这一 变化，否则，便拒绝这个改变，回到第（3）步，重新选择试探波函数的改变值。（5）返回到第二步，反复循环直到能量平均值不再有明显的改变为止。,六、格林函数Monte Carlo方法,考虑坐标空间的薛定谔方程,进行能量零点的移动,使得最低的能量本征值大于零,于是，对于 可定义格林函数,而这个格林函数是非负的。这样，薛定谔方程可写成积分方程的形式,我们可以积分这个方程得到一个函数序列,我们可以设定这个序列的初始函数,从概率论的基本原则可以证明： 如果 （a）从波函数 随机抽样产生出 的一系列数值 ；

19、（b） 从每一个根据作为从 到 的过渡密度的随机产生出 的一系列数值，则 这一系列数值的密度就是 。,第二步之所以可行是因为 。用这个方法我们可以生成一系列坐标使其满足分布 。 注意到有关这一分布函数的所有知识都包含在这一座标集合 中。如，以下的数值积分可化成求和的形式：注意到最关键的一步是以 为条件从 生成 的一系列数值。,一般来说，尽管 不知道，但可以通过局部几率无规行走的方法得到所需的数值集合。我们来分析一下从 开始如何产生序列 。用本征函数展开，有于是，上面的迭代成为,上面最后一个方程说明了三个要点:1) 数值集合 的分布密度渐进地趋向于 ， 与初始的集合 无关。2）这一渐进行为与 无

20、关。3）在 中的分布在渐进分布中的贡献具有正比于 的系数。对Monte Carlo 计算来说，如果某个集合 的一个结构对所计算的结果具有 的贡献， 那么对 的积分方程的计算将比对初始密度 的方程的计算具有更小的统计误差。因此，我们将计算如下的递推公式：,可以根据Metropolis方法进行重要抽样。格林函数的计算一般通过运动方程方法和一定的近似进行。,第四节 统计物理的Monte Carlo计算,假定我们对一个处于热平衡的恒温T 的体系感兴趣。对该热力学问题我们作如下表述：设有个包含 N 个粒子的恒温的平衡态系统，我们要计算系统的可观测量A，即该物理量的平均值,上面公弌中 x 表示在相空间中的

21、态矢（例如，其坐标为各个粒子的空间位置、动量和自旋等）， 它给出该状态在相空间中点的坐标。显然上面公式计算都是涉及很高维数积分问题。在统计力学的实际问题中，只有像理想气体、简谐振子系统二维Ising模型等少数类型的问题可以解析地严格积分求解。用蒙特卡洛方法来积分时，只是在把相空间离散化时可能会引起误差，采用蒙特卡洛方法时本身存在的统计误差以及由于计算机有限字长所引起的数值舍入误差，此外一般不存在任何其他的近似误差。,商业软件简介商业软件的三个要素：1）晶体结构或原子结构的输入2） 核心计算部分3 ） 结果的输出和分析,1） 输入晶格结构及原胞的形状和大小原子的种类和数量各个原子的空间位置带电情况表面和界面情况自旋取向设定2） 核心计算计算方案的选择参数选择3） 结果的分析和输出,