Saturday, February 22, 2014

mingda1986 量子Monte Carlo模拟 对波函数投骰子次数足够多,也可以保证系统的态趋向于真正的基态。两者的共同点,就是辛钦定理。

对波函数投骰子次数足够多,也可以保证系统的态趋向于真正的基态。两者的共同点,就是辛钦定理。










量子蒙特卡罗模拟法

量子Monte Carlo法是目前计算多粒子凝聚态系统的最精确的方法。一个有N个相互作用的电子的系统,需要3N维的Schrodinger方程才能求解。这看起来根 本做不到,因此才有了独立电子近似的能带论等近似理论。然而,以随机行走为基础的量子Monte Carlo模拟,恰是直接求解多体Schrodinger方程的方法,它通过对多体试验波函数进行抽样,进而优化,获得系统的性质。

要想深入的掌握几种常用的量子Monte Carlo模拟法,需要一定的数学储备。这里只用积分,也可以对它有所理解。

需要指出的是,量子Monte Carlo模拟法并没有很多物理图像,而是由很多技术细节组成,使得计算效率增加,精度提高等等。

一、赌博和随机模拟为什么赌博的原理能成为解决量子多体物理的基础?举个例子就知道了。
在 拉斯维加斯,赌大小是最简单的玩法。想必大家都知道,简述如下:同时掷三个骰子,如果点数之和<=9算小,>=10为大。这样看来,胜负概 率相等,赌场没太有赚头(见注1),所以又加了一条,三个骰子数值相等时,庄家赢。就是这小小的1/36的概率差,保证了押注次数足够多时,庄家必胜。投 骰子次数足够多,可以保证庄家获胜;对波函数投骰子次数足够多,也可以保证系统的态趋向于真正的基态。两者的共同点,就是辛钦定理。

二、Monte Carlo法Monte Carlo法就是用统计来模拟事件的方法。一个简单的例子是求圆周率pi,如图1. 如果想知道圆周率pi,可以在二维平面上边长为R的方块内,均匀撒点。然后计算点的坐标到正方形中心的距离。如果距离小于R,则点在圆内,反之点在圆外。 圆周率可以从落在圆内的点的数目和总的点数目的比得到。那么点恰好落在圆周的时候,算是圆内还是圆外呢?其实没有影响。撒点足够多时,恰好落在圆上的点相 对于播撒的总点数趋近于零测度。



Monte Carlo法在物理上有很多应用(更别提其他领域了),比如激光器发射的光子呈Poisson分布,可以被随机模拟;再比如电子显微镜里的电子束打进样本,发生散射,其散射截面也是概率性的,从而也可以被模拟。但是这里,为了研究量子多体系统,Monte Carlo的威力是可以计算多重积分。

除了积分表里的积分及其变种,其他的积分都需要数值求解。一个重要的里程碑是Simpson公式,定积分可以近似成:



这 样,一维积分就可以分成很多小段,每个小段都可以用Simpson公式写出。然而,在高维积分情形下,Simpson公式的误差正比于M^(-4/d), 其中M为小段的数目,d为维度。这样看来,对于高维积分,其误差下降非常缓慢,运算效率很低。这时就体现出了用Monte Carlo法计算多维积分的威力。至于为什么选择计算积分做例子,是因为它恰恰和多体理论的量子Monte Carlo法紧密相关。方法如下:
1) 假设要计算一个N重积分I=∫g(R)dR,其中g(R)是被积函数,R是3N维矢量
2) 令P(R) 是任意概率密度函数,其值为正,定义新函数f(R) = g(R)/P(R),则可以把积分重写为I=∫f(R)P(R)dR.
3) 此时,对概率密度P(R) 进行采样,获得采样点R1, R2….RM。换句话说,R1, R2…是按照概率密度P(R) 分布的M个点(其实每个R都是矢量,但是叫它们为“点”无伤大雅)。
4)这个情形下,积分I就可以近似的写成I ≈ {f(R1)+f(R2)+….+f(RM)}/M即积分的值近似等于所有采样的算数平均值。表面看起来,这样胡乱撒点,很难保证计算精度,其实不然。想 想求圆周率时,均匀撒了足够多的点,最后可以以任意精确度逼近pi,这里按照P(R)概率密度函数撒了M个点,然后当M足够大时,这些求和的平均值,就会

趋于真实的积分值。
这个结论,就是中心极限定理。据此,其计算统计误差σ,与采样点的数目M的平方根成反比,σ∝1/sqrt(M)。这样,只要采样足够多,误差可以控制在任意精度之内。以上求积分的过程,就成为了量子Monte Carlo里最简单的变分Monte Carlo的基础。

三、变分Monte Carlo 法

能够完全描述N个相互作用电子的系统的是3N维的波函数。假设系统处于任意的试验波函数ΨT(R), 则变分原理告诉我们,处于这个试验波函数下的系统的能量,一定不可能小于基态能量:
1)

左 边式子积分号里被积函数,就相当于上一节Monte Carlo积分法里的被积函数积分g(R)。2)对它做变换。联想到波函数的模方恰是概率密度,从而令P( R)= |ΨT(R)|2,此P(R)就可以当作刚才Monte Carlo积分里面的概率密度,则积分被重写为


这已经上一节从1)到2)的变换很像了。不难看出,这里相当于上一节f(R)的是 ,这里叫它EL­(R),局域能量。因为它有能量的单位,但是又和R相关。
3)此时对波函数进行采样,得到很多符合概率密度P(R) 的采样点 R1, R2….RM。。
4)从而,变分能量也可以写成E变分 ≈ { E­L(R1)+ E­L(R2)+….+ E­L(RM) } / M这样就得到了关于试验波函数的能量平均值。
5) 从E变分到基态能量: 尽管求得了变分能量值,但该能量并非基态能量。接下来该怎么获得基态能量呢?奥秘在波函数里。波函数除了是N个电子的空间坐标R(R为3N为矢量)的函 数,其实还带有一些待定参数。这些参数会被优化,使得变分能量尽可能的小,更加接近基态能量。如果采样采的好,并且试验波函数选取的好,那么变分能量的最 小值,就会非常接近基态能量。如果用这个方法求分子的结合能,可以把 这 样,波函数带有的参数优化之后,整个变分Monte Carlo过程就结束了。

变分Monte Carlo最早用于求解量子液体He,作者McMillian发表于《物理评论》。对于这项重要的工作,他的导师,因三极管和超导理论获得两次诺贝尔物理 奖的Bardeen,竟然谦虚的没有放上自己的名字。然而,还有一些细节没有提到。比如如何获得采样,以及试验波函数如何选取。这两个问题对于正确的解答至关重要。

四、如何获得呈现P(R)分布的采样点
Metropolis 算法可以实现采样。Metropolis仅仅是原始作者里的姓氏最靠前的一个人,而不是主要贡献者,但是还是以第一作者的位置把名声流传下来。不同的采样 不是假设两组采样各自独立,而是一个采样可以从上一个采样生成。假设已经有一组采样R,即第一个电子的坐标为r1, 第二个为r2 … 第N个为rN,1)把现在的电子坐标R按照任意的概率密度函数T(R’ßR) 移动到新的电子坐标R’。2)设定接受/拒绝新的坐标R’的机制,R’有一定概率被拒绝作为新的采样,也有一定概率接受。这个概率与T和P都有关。下 图就是以N=6个粒子为例子,简述从现有采样点R过渡到R’的过程,由于采样是对每个电子的3维坐标采样,这里也仅仅简化为1维坐标。每一次并不是全部的 粒子都运动,好比N个电子的坐标构成一个聚合物,如果同时聚合物的每个单体都运动,那么就会耗费很多能量;聚合物的运动仅仅是一部分的运动,和整体平移的 运动。然而如果每次不是所有单体都运动,那么采样后与采样前的坐标相差不大,这样的采样效率很低。从而,更好的办法是先把整体都平移,如b),再选取一部 分粒子,仅对这部分粒子坐标(相当于聚合物的单体)做按照概率密度T的随机行走,注意T是与坐标有关的函数(如c图)。当窗口内的坐标行走之后,再利用
回复
        Metroplis算法,判定是否选取行走之后的坐标,为新的坐标R’。




        有一点值得注意,如果给出算法的表达式(不复杂),可以发现,这样的采样,在足够长的采样步数(时间),即趋于平衡,之后,服从P(R)的概率分布,而且与T无关。


        五、试验波函数的选取显然,在变分法中,波函数的选取有很大的任意性。只要保证基本对称性,选取任何波函数都可以。然而,波函数的好坏,决定了运行的效率,更决定了能量的下限到底能够多么趋近与基态能量。这时就体现了选择波函数的重要性。

        电子作为费米子,如果只考虑独立电子的交换效应,那么波函数就成了Slater行列式——两个自旋相同的电子不能占据空间的同一位置。假设电子之间的相互作用也被考虑进去,新的试验波函数就可以从Slater 行列式得出。

        假 设电子i和电子j的距离为rij, 势能为u(rij), 总的势能就可以对所有的势能求和,∑u(rij)。由于势能高的位置出现的概率较低,从而波函数写成势能的指数形式:Exp{ -∑u(rij)},把这样的一个式子与Slater 行列式相乘,就得到Slator-Jastrow的试验波函数。

        六、投影Monte-Carlo模拟法
        尽 管变分Monte Carlo法已经可以很逼近真正的基态,但是仍旧有些不足。一个最大的缺陷就是,试验波函数的选取很大的决定了结果的好坏;变分后的能量与真实基态能量有 一定差异,因此它不是预测相图的好方法。这时,在变分Monte Carlo的基础上,另一种量子Monte Carlo法——投影Monte-Carlo模拟,可以提供更高的精度。投影Monte-Carlo法,把试验波函数做时间演化,足够长时间之后,试验波函数的投影出来的就是真正的基态波函数,证明如下:

        1. 假设Hamiltonian的本征方程为:Hφi=Eiφi,它满足封闭关系: ∑|φi><φi|=1.

        2.假设试验波函数是Ψ(0), 把它按照封闭关系分解:| Ψ(0)>= ∑|φi><φi|Ψ(0)>.

        3. 在虚时间内,薛定谔方程相当于一个扩散方程,从而时间演化写成| Ψ(t)>=exp[-t(H-E)] | Ψ(0)>= exp[-t(H-E)] ∑|φi><φi|Ψ(0)>= ∑ exp[-t(Ei-E)] |φi><φi|Ψ(0)>.

        4. 注意到所有可能的试验波函数的本征值E都大于等于基态E0, 而从变分Monte-Carlo得到的试验的本征值E非常接近E0, 当t趋于无穷大时,由于指数项exp[-t(Ei-E)]的存在,除了基态以外的其他态都衰减为零。从而,投影剩下的态为真正基态.,即

        5. | Ψ(t→∞)> = exp[-t(E0-E)] |φ0><φ0|Ψ(0)>这里选用虚时间的作用也不难看出:尽管理论上,实时间里的无穷频率的振荡项也趋于零,但是在模拟采样过程中 根本无法实现。这样,就证明了一个中心结论: 虚时间下,对于波函数在无穷远的或者将来的时间演化趋于系统的真正基态。由此,如果已经获得了从变分Monte Carlo的波函数,把它做时间演化,那么新的态将更加接近基态,从而减少系统误差。

        七、如何对试验波函数做时间演化
        第一步,首先做变分Monte-Carlo,获得已经初步优化过的试验波函数。在下图的一维势场的例子,甚至不必优化试验波函数,直接令试验波函数为均匀分布即可。

        第二步,为了使用计算机计算,把试验波函数离散化。如果是均匀分布的试验波函数,那么只需要写成一系列Delta函数的和:Ψ(0)=∑δ(x-xi)对于每个Delta函数,都各自独立的进行随机行走(扩散本身就是一种随机行走),把每个delta函数叫做一个步行者。
        第 三步,对这个试验波函数进行时间演化如果没有外势场,那么虚时间里的薛定谔方程就相当于扩散方程,从而每个Delta函数形式的步行者,会演化为一个 Gaussian波包,即扩散方程的解。假设步行者的中心从xi扩散至x, 这个过程相当于把Delta函数展宽,其格林函数写成是:


        此 时时间演化后的波函数就会变为Ψ(t)=∑G(x-xi)在有外势场的情况下,可以证明,在时间步长t比较小的时候,传播子只需要多乘以一项P(xi, x)=exp[-(V(xi)+V(x)-2E))t],就可以做总的时间演化。这样,在势场的作用下,获得的新的波函数,就可以写成 Ψ(t)=∑G(x-xi)P(xi, x)


        然而,实际应用中,每个步行者可以死亡或者生殖;其概率根据P来决定。从P的形式可以看出,如果势场比较高,那么P值比较小,步行者就有较高的死亡 率。如果势场比较低,P值较大,从而步行者有更多的机会存活。经过稍微改进的生存/死亡算法,最终就可以得到经过生死劫的步行者,也就是一系列的格林函 数,它们的和,就最终会收敛到真正的基态波函数上。

        七、结语限于篇幅,很多重要的信息没有来得及介绍。比如费米子的符号问题以及解决办法,路径积分Monte-Carlo如何推广到有限温度,甚至应用于等离子体,以及用它获得相图,或者求解量子固体氢和液态氦的实际例子。这些将适时增补。


        注1:即使概率相同,比如翻硬币,那么本金大的也有更大的概率赢。这是由于涨落随着投掷次数增加而增加。假设赌场本金为无穷大,长久看来,也是赌场必胜。

        收起回复
              回复
                    除了知道蒙特卡罗法外后面的完全不懂。。。标记先
                    回复
                          回复
                                我想问问这个方法能处理多大尺度的体系,计算效率如何?比如相对密度泛函方法要慢多少
                                收起回复
                                      • sugar0763: 我扫过一眼用量子蒙特卡洛算量子散射的文章,大概是10年前的,周期表前两行的原子能有几个吧
                                        2012-8-2 05:04 回复
                                      • admministrater: 明知故问?Monte Carlo可以处理较小的生物分子嘛。
                                        2012-8-2 07:05 回复
                                      • sugar0763: 回复 @admministrater :你说的那个是不超过经典力学范围的,量子蒙特卡洛就不是这么回事了
                                        2012-8-2 07:32 回复
                                      • admministrater: 回复@sugar0763 :哦,抱歉我没仔细看帖…
                                        2012-8-2 08:14 回复
                                      • mingda1986: 这个方法的计算scale是N^3+x*N^4,其中x~0.001,N=电子数目。但是效率比dft慢很多,计算规模也要大很多。为了算得有意义的结果经常要求助超级计算机
                                        2012-8-2 11:05 回复
                                        • mingda1986:应用也很多,包括地球物理,那些高压的相,用dft都不准,所以需要求助于它。还有等离子体,这里没有介绍有限温度的路径积分monte carlo.
                                          2012-8-2 11:07 回复
                                        • ボケ少女:回复 @mingda1986 :DFT没法处理大量的弱相互作用 QMC的方法我大概看了一些,分很多种, 但是具体的应用还不太了解,比如计算什么样的性质?能处理周期性么?超级计算机要用多大规模(多少core/flip?),能处理多大的体系要花多少时间?(大概是技术细节,这些东西paper里应该不会写的)
                                          2012-8-2 11:49 回复
                                        • ボケ少女:回复 @admministrater :QMC我不太熟悉的 不过有点感兴趣
                                          2012-8-2 11:49 回复
                                        • admministrater:回复@ボケ少女 :我当MC了…我错了。QMC我仅限于听过名词而已。
                                          2012-8-2 12:07 回复
                                        • mingda1986:回复 @ボケ少女 :科学本来就由细节组成。这些“细节”当然在paper里会提到。我不懂你说的什么“计算什么样的性质”,因为优化了波函数,就什么性质都有了。计算规模如你所说,和哪种qmc有关,如果是variational,尺度比如比dft慢10倍,如果是diffusion,慢100倍甚至更多,但是精度更高。
                                          2012-8-4 08:44 回复
                                      • 还有17条回复,点击查看
                                      回复
                                            我晕 高等数学不会啊
                                            收起回复
                                                  看不懂。。。
                                                  回复
                                                        正在加载中,请稍候...
                                                        正在发豆:54万拜年领T豆
                                                        还没签到呢!
                                                        我在贴吧
                                                        在线奖励:
                                                        • 0'
                                                          10'
                                                          00:00
                                                        • 30'
                                                        • 50'
                                                        • 90'


                                                        吧友热玩游戏排行






                                                        No comments:

                                                        Post a Comment