（理论物理专业论文）4he液体基态能量的蒙特卡罗方法研究.pdf

摘要 摘要 量子多体问题是目前理论研究的一个热点问题，低温时人们认 为二体相互作用是多体问题中的基本作用形式。4 胁液体是一种典 型的玻色多体系统。蒙特卡罗数值模拟是一种国际上通用的用来研 究多体相互作用的权威方法。一般地讲，可以分为经典与量子两种 蒙特卡罗方法，主要是依赖所研究的多体系统的运动是附合量子行 为还是经典力学的理论而定。对于量子多体系统而言，又可以分为 变分蒙卡方法、格林函数蒙卡方法以及d f m c ( d i f f u s i o n m o n t e c a r l o ) 。在不同的问题中人门可以采用不同的方法进行计算。 本文利用变分蒙特卡罗方法对系统的基态能量进行了比较详细地 研究。通过调整不同的二体相互作用函数，发现正确的二体相互作 用所满足的关系。本文中我们对三种不同二体相互作用势时4 胁液 体的基态能量的上限进行了计算估计。结果表明：m s l 2 g 6 与实验结 果最为接近，l e n n a r d j o n e s6 - 1 2p o t e n t i a l 和h f d h e 2 也可得到 与实验值相近的结果。 关键词；变分蒙特卡罗方法；4 胁；勒纳德一琼斯势；m sj 2 g 6 ；h f d h e 2 l i i a b s t r a c t t h eq u a n t u mm a n yb o d yp r o b l e mi so n eo ft h eh o tr e s e a r c hp o i n t i nt h ea l lo ft h ew o r l d w h e nt h et e m p e r a t u r ei sv e r yl o w , t h et w ob o d y i n t e r a c t i o np l a ya n i m p o r t a n tr o l e i nt h es y s t e m t h em o n t ec a r l o n u m e r i c a ls i m u l a t i o ni so n ep o w e r f u lm e t h o dt os t u d yt h em a n yb o d y i n t e r a c t i o ns y s t e m g e n e r a l l y , i tc a nb ec l a s s i f i e da sc l a s s i c a lm o n t e c a r l oa n dq u a n t u mm o n t ec a r l o ，d e p e n d e do nt h el a wo ft h es y s t e m g o v e r n e db yt h ec l a s s i c a lm e c h a n i c so rt h eq u a n t u mm e c h a n i c s t ot h e q u a n t u mm a n yb o d ys y s t e m ，t h e r ea r ev m c ( v a r i a t i o n a lm o n r e c a r l o ) ， d f m c ( d i f f u s i o nm o n t e c a r l o ) a n dg f m c ( g r e e nf u n c t i o nm o n t e c a r l o ) t ot h eg i v e np r o b l e m ，o n ec a ns e l e c tt h ep r o p e rm e t h o dt od e a lw i t hi t t h el i q u i dh e l i u m4i so n et y p i c a lb o s em a n yb o d ys y s t e m i nt h i s p a p e r ，w ec a r e f u l l yi n v e s t i g a t e dt h eg r o u n ds t a t ee n e r g yo f t h es y s t e m b ya d j u s t i n gt h ed e t a i lf o r mo ft h et w ob o d yi n t e r a c t i o n ，w ef i n dt h e g o o df o r m u l af o rt h i st w ob o d yi n t e r a c t i o n w ei n v e s t i g a t e dt h ee f f e c t o ft h r e ed i f f e r e n tt w o b o d yp o t e n t i a l so nt h ee q u i l i b r i u me n e r g yo f l i q u i dh e l i u m4 o u rr e s u l ts h o w st h a to n ep o t e n t i a l ，m s 12 g 6 ，g a v e b e t t e ra g r e e m e n tw i t ht h ee x p e r i m e n t ，b u th f d h e 2a n dl e n n a r d j o n e s 6 12g a v ea na c c e p t a b l ea g r e e m e n t k e y w o r d s ： 、，a r i a t i o n a lm o n t ec a r l om e t h o d ( v m c ) ；4 h e ； l e n n a r d - j o n e s ( l j ) 6 - 1 2p o t e n t i a l ；m s12 g 6 ；h f d h e 2 i v 引言 引 言 蒙特卡罗( m o n t ec a r l o ) 方法又叫做随机模拟( r a n d o ms i m u l a t i o n ) 方法，有的文献中也称作随机抽样( r a n d o ms a m p l i n g ) 技术或统计试验 ( s t a t i s t i c a lt e s t i n g ) 方法“。 蒙特卡罗方法最早可追溯到十九世纪的蒲丰( b u f f o n ) ，但其名称的确定 。- 6 3 和系统发展却开始于二十世纪四十年代中期。由于电子计算机的发明，使 得用随机实验方法解决实际问题成为可能。人们可以利用电子计算机实现大 量的随机抽样试验。第二次世界大战期间，蒙特卡罗方法在原子弹研制工作 中的应用堪称当时最具有代表性的工作。在这一研究工作中，人们对中予随 机扩散的行为进行了模拟试验。即在电子计算机上模拟中子的行为，并进行 随机抽样，通过对所获得的大量有关中子行为的数据，计算并推断出所要求 的参数。然而，即使是在这项研究工作的早期，m e t r o p o l i s 、u l a m 、v o nn e u m a n n 等人”- 7 1 也并不是对中子行为进行完全的直接模拟，因为这样做所得到的结果 是不很理想的。他们采用了某些降低方差的技巧，例如“重要性抽样”、“俄 罗斯轮盘( r u s s i a nr o u e t t e ) 和分裂( s p l i t t i n g ) 方法”1 。但是，这些技 巧毕竟是在中子行为直接模拟基础上采用的。这种粒子输运问题的早期研究 工作一直延续到二十世纪五十年代初的h a r r i s 和k a h n 等人o + ”1 。 蒙特卡罗方法的基本思想可以概括如下：为求解数学、物理、工程技术 以及生产管理等方面的问题，首先建立一个概率模型或随机过程，使它的参 数等于问题的解；然后通过对模型或过程的观察或抽样试验来计算所求参数 的统训特征，最后给出所求解的近似值。而解的精确度可用估计值的标准误 差来表示。 假设所要求的量x 是随机变量善的数学期望e f 孝) ，那么近似确定z 的方法 是对害进行次重复抽样，产生相互独立的善值的序列备、岛、晶，并计 算其算术平均值： 知= 万i 善n 茧 ( o 1 ) 由柯尔莫哥罗夫( ko 月moropob ) 加强大数定理 1 1 有 4 月0 液体基态能量的蒙特卡罗方法研究 p 【溉知2 。j 2 1 ， ( o 2 ) 因此，当充分大时，下式 靠“e ( 孝) = x ( o 3 ) 成立的概率等于1 ，亦即可以用5 乍为所求量x 的估计值。 样本方差 蠢= 万11 f 万1 缶生n 舌2 一e 2 ( 鲁) ( 0 4 ) 作为理论方差的估计值。 利用蒙特卡罗方法可以解决很多问题，应用范围非常广泛，但总体来说， 根据该问题是否涉及随机过程，可以分为以下两类： 第类是确定性的数学问题。用蒙特卡罗方法求解这类问题的方法是： 首先建立一个与所求解问题有关的概率模型，使所求的解就是我们所建立模 型的概率分布或数学期望；然后对这个模型进行随机抽样观察，即产生随机 变量：最后用其算术平均值作为所求解的近似估计值。计算多重积分、求逆 矩阵、解线性代数方程组、解积分方程、解某些偏微分方程边值问题和计算 微分算子的特征值“”1 等都属于这一类。 第二类是随机性问题。例如前文所提到的中子在介质中的扩散等问题“7 ”3 就属于随机性问题，这是因为中子在介质内部不仅受到某些确定性的影响， 而且更多的是受到随机性的影响。对于这类问题，虽然有时可表示为多重积 分或某些函数方程，并进而可考虑用随机抽样方法求解。然而一般情况下都 不采用这种间接模拟方法，而是采用直接模拟方法，即根据实际物理情况的 概率原理，用电子计算机进行随机抽样试验。原子核物理问题、运筹学中的 库存问题“、随机服务系统中的排队问题“ 。”1 、动物的生态竞争和传染病的 蔓延等部属于这一类。本文所要讨论的4 胁系统就属于随机性问题。 蒙特卡罗方法是一种具有独特风格的数值计算方法，它的优点以及它与 一般数值方法的不同点可归纳为如下三个方面“： 1 ) 蒙特卡罗方法计算积分是通过大量的简单的重复抽样实现的，与一般 数值积分方法如g a u s s 求积法等大不相同，方法和程序都很简单。 2 ) 蒙特卡罗方法的收敛是概率意义下的收敛，与一般数值方法的一致收 r一上、 敛不同。由于蒙特卡罗方法的收敛速度比较慢，其主阶仅为0 in 2 l 。因此， 不能用来解决精度要求很高的问题。但其收敛速度与问题维数无关，而这正 是其它数值方法所无法比拟的，可以认为是蒙特卡罗方法最独具特色的地方。 3 ) 蒙特卡罗方法在解题时受到问题的条件限制影响较小，适应性强。 蒙特卡罗方法有许多分支，其中包括变分蒙特卡罗方法( v a r i a t i o n a l m o n t ec a r l o ) 、扩散蒙特卡罗方法( d i f f u s i o nm o n t ec a r l o ) 、格林函数蒙特 卡罗方法( g r e e nf u n c t i o nm o n t ec a r l o ) 等。本文主要介绍变分蒙特卡罗方 法，只涉及一般的函数，而另外两种方法要用到格林函数，本文暂不涉及。 然而必须指出的是：格林函数蒙特卡罗方法汁算出的结果和实验值最为接近， 而变分蒙特卡罗方法计算出的结果与实验值相差较大。本文的目的不在比较 这几种方法的优劣，而是比较几种不同的势函数，故采用比较简单的变分蒙 特卡罗方法。 4 h e 液体基态能量的蒙特卡罗方法研究 第一章变分蒙特卡罗方法 1 1 算法简介 多体问题啥密顿量的一般形式为： 疗= 一嘉；v 善r ( “) ( 2 ) 包括一个动能项和一个两体相互作用的势能项。两体相互作用与粒子间 的距离以及自旋自由度、同位旋等有关。多体问题处理起来非常复杂，它的 一个显著特征是在几乎所有人们感兴趣的多体系统中，两体相互作用势在非 常小的距离表现出极强的排斥力，因而一般的微扰方法很难直接运用。 既然微扰法不能用，只能采用其它的方法来求解( 2 1 ) 。可以近似粗略 地考虑两个排斥核的相互作用，找一个合适的试探波函数，采用变分法来计 算基态能量的上限。其中最简单的是j a s t r o w 形式。“，它由两部分组成，一 部分是单粒子轨道波函数( s i n g l ep a rc i c l eo r b i t a ls ) ，另一部分是相互作 用因子( c o r r e l a t i o nf a c t o r ) ，是所有每对粒子间相互作用的乘积。 其中一些物理上常用的形式有： 1 ) 有限玻色于系统( 液滴) 蚧= 丌p n r ) 丌，( ，；，) ( 2 2 ) 第一项是所有围绕质心的粒子波函数的乘积，第二项是相互作用因子。 2 ) 无限玻色子系统( 液体或气体) y 广兀( _ ) ( 2 - 3 ) 这种问题中，试探波函数使所有处于同一个格子中的有限数目的粒子相 互作用因子的乘积，然后根据周期性边界条件延拓至全空间。 3 ) 晶体 = 兀妒“- 4 ) n ，( _ ) ( 2 4 ) 4 是晶体中的格点，单粒子波函数将每个粒予束缚在格点周围。同样地， 这个波函数也是处于同一格子中的有限数目的粒子波函数的乘积，并周期性 延拓至全空间。 第一章变分蒙特卡罗方法 费米子系统也可以采用同样的方式描述，唯一的不同是要求波函数具有 反对称的形式。 更一般的试探波函数形式还包括三体相互作用、相互作用因子( 依赖于 动量) 、自旋和同位旋等的影响。为简便起见，以下仅考虑最简单的形式，基 态能量的上限由下面的多维积分给出： 时 昂：筚掣蜂辈 亿；， 昂。节丽丽 g 5 凡代表3 个自由度的集合。当不考虑自旋且哈密顿量是一个局域算子 式( 2 5 ) 可以改写为近似于蒙特卡罗积分的形式。定义局域能为： e l ( r ) 2 南膏( r ) ， ( 2 a ) 并引入多维概率分布函数： 小) = 赫 偿， 其中p ( r ) 恒为正值，并归一化。于是变分能量可改写为： e v = p ( r ) 丘( r ) 搬 ( 2 8 ) 这个积分是标准的蒙特卡罗积分的形式。p ( 凡) 是重要抽样函数，e ( r ) 是被 抽样的函数。利用变分蒙特卡罗方法求解能量的任务就变得简单明了：用 m e t r o p o l i s 算法产生一组随机变量蜀、尺2 、a 、粕，并计算能量的估计值： 它的方差为 b = l f e l ( r , ) ( ：，) ( 2 1o ) 注意：对e ：进行的抽样与膏2 的平均值并不相同，除非哈密顿量矗是一个局 域算子! 日e 液体基态能量的蒙特卡罗方法研究 局域能的形式并不是唯一的。哈密顿量中势能项的贡献就是势能，但是 关于动能的计算包含有微商运算，直接计算将导致所谓的三体项( t h r e eb o d y t e r m s ) 。通过分步积分技巧有可能得到比式( 2 6 ) 更为简单的形式。但是不 管由什么形式算出的变分能量都应该相同，仅仅是变分值与方差不同而已。 利用直接的方法进行计算，也就是说将l a p l a c e 算符作用在试探波函数的右 边，在极限情况下，采用波函数的精确解为试探波函数有可能得到零方差的 结果。这可能是一种比较好的方法，但正如前文所提到的，一般地会产生三 体项，使计算变得很困难。 1 2 一个简单的例子：一个处于高斯势阱中的粒子 考虑一个处于高斯势阱中的核子。 矿c ，= e x p 一( ； 2 】 c 2 一， h 2 其中= - - 4 5 m e v ，6 ：2 0 9 f h l ，素2 4 1 5 0 m 。v f mz 。基态能量可由径向薛 定谔方程数值积分得到，结果为：一1 0 8 1 3 m e v 。 一个简单而又较准确的试探波函数可选取谐振子的基态波函数。 局域能为 2 纠一) 功 e ( r ) = 嘉( 3 a 2 爿r 2 ) + y ( ，) ( 2 1 3 ) 可以很容易地算出变分能量为： 即，= 和q 播 i 亿 可以用两种方法来计算能量值。第一种是采用直接抽样的方法，波函数 服从正态分布，程序名为 一2 m ，( 详见附录) ；第二种是采用m e t r o p o l i s 算 法，产生一组随机数序列，进而得到我们想要的结果，程序名为t 4 m ( 详见 第一章变分蒙特卡罗方法 附录) 。 表1 ：高斯势阱中粒子基态能量两种算法比较 a n a l v t i cg a u s s i a ns a m p l i n g m e t r o p o l i ss a m p l i n g a ( 霞) j ，0 0 01 0 , o o o1 0 0 , 0 0 0j ，0 0 01 0 , 0 0 0 1 0 仉0 # 0 ( 毛) 0 46 8 9 41 2 0 3 7 0 3 士0 2 76 8 4 士0 0 96 8 8 8 士0 0 2 77 7 5 0 2 86 9 2 - 士0 0 96 8 4 2 - 士0 0 2 7 0 59 1 9 01 3 5 6 8 9 9 * 0 2 2- 9 1 5 士0 0 79 2 2 4 0 0 2 3- 9 2 5 土0 2 29 4 2 士0 0 79 1 8 8 0 0 2 3 0 61 0 3 0 11 3 1 21 0 0 2 士0 1 51 0 3 3 士0 0 51 0 ，2 9 8 土o 0 1 6- 1 0 2 0 士0 1 51 0 2 4 i - 0 0 51 0 3 0 1 0 0 1 6 o 7，1 0 0 6 91 2 15】o 2 l 土o 1 699 9 * 0 0 41 0 0 4 1 士0 0 1 49 9 6 土0 111 0 0 3 0 0 4一1 0 1 3 0 * 0 0 l5 o 88 5 2 51 3 0 28 0 7 土02 28 5 7 士0 0 88 5 4 0 0 ，0 2 48 5 4 土0 2 58 6 7 + 0 0 7- 8 4 7 2 0 0 2 4 表1 中是程序运行的结果比较。第- - y , 是变分参数日，单位为f m 一，第二 列是理论值，其中既有能量的期望值，又列出了哈密顿量的平方h 2 的期望值。 其余的是由蒙特卡罗方法计算得到的数值解，并且每种算法分别运行了 l ，0 0 0 、1 0 ，0 0 0 、i 0 0 ，0 0 0 次，即分别选取了l ，0 0 0 、1 0 ，0 0 0 、1 0 0 ，0 0 0 个样本 点。所有能量的单位均为m e v 。m e t r o p o l i s 算法每次随机行走的最大距离由 变量s c a l e 控制，取值为3 f m ，该算法的效率约为6 8 7 5 4 2 6 5 。 从计算结果来看，采用直接抽样方法得到的能量值与理论值比较接近， 而且抽样1 0 ，0 0 0 次又比抽样l ，0 0 0 次得到的结果要好，抽样1 0 0 ，0 0 0 次的结 果已经非常接近理论值了。m e t r o p o l i s 算法体现得较为明显。因此，抽样1 ，0 0 0 次还远远不够，必须增大样本容量。在这一问题中，标准偏差并没有多大意 义，仅能说明精度要提高一位，计算量要增加到原来的1 0 0 倍。在实际问题 中并不能预先知道精确值，所以为检验该算法的有效性，必须做大量的计算。 1 3 周期性边界条件 对个具有无限多粒子的系统来说，要想直接描述它是一件很难的事情。 处理这类问题通常是借用经典统计力学中的方法，将整个系统划分成多个边 长为上的格子( s i m u l a t i o nb o x e so rb l o c k s ) ，作为研究的基本单元。这就是 著名的周期性边界条件( p e r i o d i cb o u n d a r yc o n d i t i o n ) 每个格子的边长及每个格子中所包含的粒子数要根据密度p 选取： 每个格子中的粒子数 ，、 p2 r 一13 j l 粒子数一般考虑晶格结构选取，当密度较高时，可以考虑成是面心立方结构 ( f a c ec e n t e r e dc u b i c ) 或者六角密排结构( h e x a g o n a lc l o s ep a c k e d ) 。 以面心立方晶胞为例，8 个顶角上的基元由8 个晶胞共享，6 个面心上的基元 4 h e 液体基态能量的蒙特卡罗方法研究 8 三+ 6 土= 4 由2 个晶胞共享，这样每个晶胞中包含8 2 个基元。如果每个格子 包含r 1 个晶胞，那么基元总数为4 玎。于是就得到一系列“幻数( m a g i c n u m b e r ) ”，3 2 、1 0 8 、2 5 6 、5 0 0 、8 6 4 大多数情况下，用蒙特卡罗方法解 决玻色子系统中的问题都以这些数字中的一个作为粒子数。 费米子系统中的幻数又与玻色子系统中的大不相同，与单个粒子的动量 的量子化以及给定费米动量之下粒子的占据状态等都有关系。动量的可能取 值用三个整数【q ，”z ，他j 来描述，这三个整数可正可负，动量的大小取决于到 原点的距离，p , h , 1 1 2 + , 2 2 + , 3 2 。原点【o ，o ，o ) 对应一个状态，到原点距离为1 的有6 个点( 1 ，o ，叭、( o ，1 ，0 ) 、( o ，0 ，l l ，对应六个新的状态。距离再远一 点，j ，有1 2 个点( 1 ，1 ，叭、( 1 ，o ，、( o ，1 ，土，以此类推。累计数分 别为：1 ，7 ，1 9 ，2 7 ，3 3 ，5 7 再乘以系统的简并度( 2 或4 ) 。 采用周期性边界条件要遵循一定的规则。粒子即可能跑出所研究的格子， 也j 能有格子外面的粒子进入格子( 例如，m e t r o p o l is 算法被接受的走法总 是有使粒子走出格子的趋势) 。为了检测粒子所走的步子是否能被接受，应该 考虑所有格子中由于浚粒子运动引起的波函数的改变。为使计算简便可行， 通常假设当距离大于兰 2 时，粒子间的相互作用因 子都相同。可以在所考虑 的粒子周围画一个半径 ， 为；的球，只计及球内 z 粒子对波函数的贡献，不 管这些 粒子是否在所考虑的格 子内。如果某一个粒子被 包含在球内，那么其它周 期同样位置的粒子就不 可能再出现在同一个球 内，也就是说每个粒子最 多被考虑一次。因此为了 计算由于粒子运动引起 q q 0 t - _ _ ； ； 、。 | j _ _ _ ，、 | o - - 。 | 。鸷 ：、 。| ，； ：。爹遥 j 一一”“。、菇 。”，。，。” ， 、一 ； 。af 图1 粒子f 的最近邻粒子 第一章变分蒙特卡罗方法 的波函数的改变，必须通过平移在所考虑的格子里找到该粒子最近邻的所有 粒子。如图l 所示。 一旦采用m e t r o p o l i s 算法将所有粒子都检测过之后，就可以计算出这种 新的分布下的局域能。这个局域能包含两项，一项是相互作用势( i n t e r a c t i o n p o t 8 i 3 t i a l ) 的贡献，另一项是由于相互作用因子( c o r r e l a t i o nf a c t o r ) 所产 生的。在生处对相互作用因子所作的假设导致了只有球内的粒子才对动能项 2 由贡献。 注意：必须保证相互作用因子及其导数在妻附近光滑，这样得到的动能 z 才正确。为了能解决这个问题，最简单可行的办法就是将三取得足够大，使得 相互作用因子在大于当时可以忽略。这样做的代价就是所要计算的粒子数增 z 多，尤其是当密度很大时。 当计算势能时就不会出现这样的问题。我们只需考虑格子中每对粒子之 间的相互作用。如果相互作用是短程力，可以不考虑某个截断距离之外的所 有粒子。但如果相互作用是长程力时，必须将所有的粒子都考虑进去。 4 h e 液体基态能量的蒙特卡罗方法研究 第二章变分蒙特卡罗方法在4 h e 系统中的应用 z 1 物理背景 当描述简单液体的性质时，一般假设系统中的相互作用主要是两体相互 作用( t w o b o d yf o r c e s ) 。从实验中得到的数据进一步证实了两体之间的相 互作用，因此理论物理学家正在尝试寻找一些新的势函数，以便能更好地解 释实验结果及应用在理论计算中。通常在高密度系统的计算中，如果还包括 三体相互作用，那么计算结果就会与实验结果吻合得很好。以液氦系统为例， 高阶项“7 1 的重要性远不如两体相互作用。 在4 m 系统的理论研究中，最常用的势函数是l e n n a r d j o n e s ( l j ) 6 - 1 2 p o t e n t i a l v ( r ) 。”( 州 占：1 0 2 2 k 仃= 2 5 5 6 a ( 3 1 ) 其中的参数是d eb o e r 和m i c h e l s 。”由从s e c o n dv i r i a lc o e f f i c i e n t s 得到 的。 除了l e n n a r d - j o n e s ( l j ) 6 - 1 2p o t e n t i a l 外，还有很多不同的势函数， 本文还将讨论另外两种势函数一一m s l 2 g 6 和h f d h e 2 详细的讨论可参阅 a z i z 2 4 等人的文章。m s l 2 g 6 是由m a i t l a n d 和s m i t h 提出的，将更一般的 l e n n a r d j o n e s ( l j ) p o t e n t i a l 加以改进为： 参数为 叶摩6 盼( 爿( 讣，埘 吲吣州：拍( 刊, r = = 2 9 6 7 盖。 t t f d h e 2 源于a h l r i c h s “5 1 等人提出的势函数，其形式比较复杂： 嘶叫4 唧一心 6 + 制蝙( 训砷) ， 第二章变分蒙特卡罗方法在4 h e 系统中的应用 f ( r ) 唧愕一- 2 ， 二d r 胛 1rd 参数值也是由a z iz 给出的。 a = 0 5 4 4 8 5 0 4 1 0 6 = 1 0 8 k 口= 1 3 3 5 3 3 8 8 4 c 6 = 1 3 7 7 3 2 4 1 2 d = 1 2 4 1 3 1 4 c 8 = 0 4 2 5 3 7 8 5 o r = 2 9 6 7 3 a c 1 0 = o 1 7 8 1 0 0 a z i z 等人的结果表明：h f d h e 2 给出的结果与实验符合的最好：m s l 2 g 6 与h f d h e 2 差不多，给出的结果也能接受；但没有一种势函数能够与实验数据 完全吻合。 r c a j 图2 ：三种势函数曲线图 2 2 变分蒙特卡罗方法 早在1 9 6 5 年，m c m i l l a n 啪1 就已经研究过4 胁系统的性质，并且这是第一 个用变分蒙特卡罗方法来研究的量子系统。他当时所用的势函数就是 l e n n a r d j o n e s ( l j ) 6 - 1 2p o t e n t i a l 。本文将同样采用l e n n a r d j o n e s ( l j ) 6 1 2p o t e n t i a l ，并同时采用m s l 2 g 6 和h f d h e 2 ，对比这三种势的结果。 两体相互作用因子的形式为： 4 h e 液体基态能量的蒙特卡罗方法研究 m ) - e 水圳 8 。， a = 1 6 e m 1 a t t s 盯 为了方便数值计算，根据式( 2 3 ) 试探波函数的平方的形式可写为： 妒；= e x r 一若( 专 5 c ，s ， 。：鱼矗：生 其中 = 5 ”，盯，”盯，取盯为长度基本单位( 下同) ，相应地局域能 毋= 若降( 而) e ， 势函数取三种不同的形式，分别为： 咻) _ 4 s b 虿1 ， 晰一蚴孔熹) 蚓 s ， 删z + s 睁卜半 矿c 而，= g _ o x p 、- 口名x , , ，一 c ：( 朝6 + c ；( 詈 8 + c l 。( 詈 ” f c b ， 。， 第二章变分蒙特卡罗方法在4 h e 系统中的应用 f k ) =唧悔1 2 ， 1 重 d a = 0 5 4 4 8 5 0 4 1 0 0 a = 1 3 3 5 3 3 8 8 4 d = 1 2 4 1 3 1 4 2 9 6 7 3 乇2 变分参数为c ，由于程序 运行一次消耗的时间很长，所 以根据m c m i l f a n 等人的结果， 当变分参数取为c = 1 1 6 ，密度 取为实验值的0 9 倍时，计算结 果就达到最小值，与实验值符 合得很好。 2 3 程序设计 本文所用到的程序均采用 m a t l a b 程序语言设计。 主程序流程图如图2 所示 所有详细程序见附录( t 5 j m 、 c l o s e m 、e l _ m 、e 2 m 和 e 3 m ) 。 主程序曲j m ：首先做计 算前的准备工作，将内存清空。 然后设置主要参数，并计算一 些以后要用到的常数备用，这 样可以节省机器时间，提高程 序运行的效率。h 为粒子数参 数，用于设定粒子数印，当 h = 3 时，n p 2 1 0 8 ，n p 必须为 占= 1 0 8 c 6 = 1 3 7 7 3 2 4 1 2 c 8 = o 4 2 5 3 7 8 5 c l o = o 1 7 8 1 0 0 闰3 ：主程序设计流程图 4 月0 液体基态能量的蒙特卡罗方法研究 前文所提到的幻数。c 为变分参数，根据m c m i l l a n 的结果取为1 1 6 。柏为格 子数，三维情况下为2 6 。”m 为随机行走的步数，至少为1 0 0 0 步。d e n s 为密 度，根据m c m i l l a n 的结果取为实验值的0 9 倍。，为格子的边长。s e a l 控制随 机行走的范围。粒子的位置坐标储存在矩阵。( 3 , 1 0 8 ) 中，3 代表三维坐标，1 0 8 为粒子标号。粒子的初始位置随机给定，随机数种子设为1 ( 种子可任意指定) 。 准备工作完成后，即可开始计算。 第一步首先计算中心格子内所有粒子局域能的总和。然后再计算下一个 格子，最后对所有格子求平均值。其中第零个格子的计算用来产生预期 的粒子分布。 在每个格子的计算中。从第个粒子开始直至最后一个粒子，即第1 0 8 个粒子依次进行检测。每个粒子在原有位置的基础上随机迈出试探性的一步， 然后检测这一步是否走出中心的格子。如果走出了格子，根据周期性边界条 件平移回中心格子。接着计算由于该粒子位置的改变对波函数的影响，根据 m e t r o p o l i s 算法决定该粒子的位置是否变化。当所有粒子都检测过遍之后， 计算这时的局域能并保存。这时所有粒子均随机行走一步，然后按照同样的 方法，再随机行走所设定的步数。 当所有的格子都算过之后，取平均值作为计算的结果，并计算方差。 子程序c l o s e m ：计算两个粒子之间的距离，用于计算相互作用因子。如 果粒子间距离大于l 2 ，均取为l 2 。所返回的值为在中心格子中两粒子之间 的“距离”，并且恒为正值。 子程序e 1 m 、e 2 m 和e 3 m ：分别对应于三种不同的势函数，用于计算 相应的局域能。动能在l 2 处被截断，而势能项的截断点在处，因为势能项 的贡献大。由于子程序c l o s e m 所返回的值恒为正值，因此只考虑该粒子周 围8 个格子中的粒子对局域能的贡献。 2 4 结果与讨论 程序运行的结果如下表所示： 表2 中第一列为实验值，第二列为格林函数方法结果。“。后六列为程序 运行给出的结果，其中l e n n a r d j o n e s6 - 1 2p o t e n t i a l 给出了几种不同步数 下的结果，可以看出随机行走1 0 ，0 0 0 步与随机行走5 ，0 0 0 步的方差相差不大， 而随机行走l ，0 0 0 步的方差比随机行走5 ，0 0 0 步的方差大很多。因此选取 5 ，0 0 0 步较为理想。同时对变分参数稍作改变后给出的能量值就变大，说明变 第二章变分蒙特卡罗方法在4 h e 系统中的应用 明变分参数选取1 1 6 也是比较理想的。从三种势函数给出的结果看，m s l 2 g 6 给出的结果最好，l e n n a r d j o n e s6 - 1 2p o t e n t i a l 次之，h f d h e 2 最差。 图4 和图5 是随机行走前后的粒子分布图。 图4 ：随机行走前粒子分布图图5 ：随机行走后粒子分布图 表2 ：三种势函数下计算结果比较 e x p e r i m e n tg f m c l e n n a r d - j o n e s6 - 1 2p o t e n t i a lm s l 2 g 6h f d h e 2 c 1 1 61 1 61 21 1 61 1 61 1 6 n m l ，0 0 0 5 ，0 0 05 。0 0 01 0 , 0 0 01 0 。0 0 01 0 ，0 0 0 - 6 8 45 8 9 135 8 8 5 55 6 3 9 95 8 9 7 86 3 9 4 85 8 7 1 8 e ( k ) 7 1 4土 士士 0 0 2 0 ，0 2 7 000 1 0 3 0 0 0 9 9o0 1 0 20 0 1 1 90 0 1 1 9 d e n s i t y ( ( r 。) 0 3 6 4 803 6 50 3 3 3o 3 3 30 3 3 30 3 3 303 3 30 3 3 3 4 l i e 液体基态能量的蒙特卡罗方法研究 附录 m a t l a b 程序 程序1 ：( t 32 m ) t 32 m v a r i a t i o n a lm o n t e c a r l om e t h o dw i t hg a u s s i a np o t e n t i a l t r i a lf u n c t i o ni sh 0 w i t hp a r a m e t e rap r o pt oe x p ( - 0 5 ( a + r ) “2 ) c l e a ra l l ； 一一一一- - d = 2 0 7 5 ； h b a r 2 m n f 5 ： n i = 10 0 0 0n u m b e r0 fs a m p l e s - - - 一一 s e e d = 1 ； r a n d ( s t a t e 。，s e e d ) ； 一 f o r 净1 ：5 a = 0 3 + 0i + i ： e = o ： e 2 = 0 ； e r r = 0 ； a 2 = a * a ； t = 1 a 2 ： 一 f o rk = l ：n i s e e d v a r i a t i o n a lp a r a m e t e ra p a r a m e t e ro fg a u s sd l s t r g e n e r a t ex f i 2 2 + p i + r a n d ； r = s q r t ( - t + l o g ( r a n d ) ) ； 一一一一一一一一一- i f r a n d 0 5 x o 、= r + c o s ( f i ) ； e l s e x ( 1 、= f + s i n ( f i ) ； e n d 一一一一- f i 2 2 + p i + r a n d ； r = s q r t ( - t 1 0 甙r a n d ) ) ； x ( 2 ) = r 奉c o s ( f i ) ； g e n e r a t ey x ( 3 、= 一s i n ( f i ) ； g e n e r a t ez 一一一一， v p = 一4 5 4 e x p ( x x v 4 3 6 8 1 ) + d + a 2 + ( 3 - a 2 + ( x + x 。) ) ；l o c a l e n e r g y 6 附录 e 5 e + v p ； e 2 = e 2 + v p + v p ； e n d 9 扣一一一一 e a f - - e h i ； e a v 2 = e 2 n i ； e r r = s q r t ( ( e a v 2 一o a v + e a v ) ( n i - 1 ) ) ； 一一一一一一一一 a ( i ，：) = 【ae a re a v 2e r r ； e n d 一一一一一 a 程序2 ：( t 4 m ) t 4 m v a r l a t i o n a lm o n t e c a r l om e t h o dw i t hg a u s s i a np o t e n t i a l t r i a lf u n c t i o nl sh o w i t hp a r a m e t e rap r o p t oe x p ( - o 5 + ( a + r ) “2 ) m e t r o p o l i ss a m p l t n g c l e a ra l l ； d = 2 0 7 5 ； n f ： h i = 10 0 0 0 0 s c a l e = 3 ； h b a r 2 m n u m b e ro fs a m p l e s 一一一一一一一一一 s e e d = l ：s e e d r a n d ( s t a t e , s e e d ) ； 一一一一- - f o r i - 1 ：5 一一一一一一 a = 03 + 0 m v a r i a t j o n a lp a r a m e t e ra a 2 2 a + a ： c l = 3 + d a 2 ： c 2 = d + a 2 + a 2 ： e a v = 0 ； e a r 2 = 0 ； n a = 0 ； 一一一一一一 x = r a n d ( 1 ，3 ) - 0 5 ； g e n e r a t r ex r = x + x ： 1 7 4 h e 液体基态能量的蒙特卡罗方法研究 f i o l d = e x p ( a 2 + r ) ： e = 0 ： e 2 = 0 ： e 名一一一一一一一一一 f o rk = l ：n i x p = x + s c a l e + ( r a n d ( 1 ，3 ) - 0 5 ) ； r p 2 x 矿x p ； f i n e w = e x p ( - a 2 + r p ) ； u = f i n e w f i o l d ； 一一一一一 i fu = r a n d x = x p ； r = r p ； f i o l d = f i n e w ； n a = n a + ! ； e n d 一一一 e n = - 4 5 * e x p ( r 4 3 6 8 1 ) + c 1 一c 2 + r ； e = e + e n ； e 2 = e 2 + e n e n ； e n d 一一一一一一 e = ：e n i ： e 2 = e 2 r f i ； e r r = s q r t ( ( e 2 - e + e ) n i ) ； 一一一一一一一一一 a ( i ，：) 2 【ae e i tn a n i ； e n d 一一一一一一- a 程序3 ：( t 51 m ) t 51 m m c m i l l a n sv a 砌a t i o n a lm o n t e c a r l oi nl i q u i dh e l i u m - 4 w l t ht h r e ed i f f e r e n tp o t e n t l a l l e n n a r d j o n e sp o t e n t i al m s l 2 g 6a n dh f d h e 2 g l e a ra l l ； 一一一一一一 8 m a i np a r a m e t e r n = 2 ： n p