（理论物理专业论文）约束条件下的硬球流体.pdf

上海交通大学硕士学位论文 约束条件下的硬球流体 摘要 f 对于物态，一般可分为固态、气态和液态，对于前两者，我们 、 的研究方法相对比较完善了，但是对于液态，研究的手段还不是很 多，主要有以蒙特卡罗方法为代表的数值方法和积分方程方法为代 表的解析方法。除此之外，二十世纪六十年代，h o h e n b c r g 和k o h n 两人提出了密度泛函理论，称为d f t ，m c r m i n 其后把这一理论推 广到有限温度情形。密度泛函理论认为，系统的自由能是密度的泛 函，只要知道了系统密度分布的形式，也就知道了自由能的形式， 其它的物理量也就知道了。此外，d f t 认为系统的密度是由外场唯 一决定的，即一个外场唯一的决定了一个系统的密度分布，一个系 统的密度也唯一的对应了一个外场。这一方法的出发点是，先找到 系统的自由能，然后再求其它物理量。此后，d f t 理论不断得到完 善，对于硬球流体r o s c n f e l d 给出了一个泛函形式，得到了广泛的 应用，被称之为r o s c n f e l d 密度泛函理论。密度泛函理论在国外正 ； 越来越多的用于解决流体问题，而且日趋成熟寸国内目前应用密度 泛函来解流体问题的并不多，本文就是利用局域密度泛函方法对约 束条件下的硬球流体进行研究。 上海交通大学硕士学位论文 我们用r o s e n f e l d 密度泛函理论，同时用分子动力学方法，对 处于两平行硬墙之问的硬球流体的密度分布进行了计算。通过比 较两种方法的结果，发现在墙之间距离较大时，r o s e n f e l d 密度泛函 理论的结果与分子动力学模拟的结果符合很好；当墙的距离很小 时，这两个结果之间存在明显的不一致。我们还研究了约束条件下 密度分布的结构。 在本文中首先介绍了密度泛函理论，然后介绍了r o s e n f e l d 密 度泛函理论和用这一泛函形式所研究问题得到的结果，最后介绍一 下在计算过程中用到的共轭梯度方法。 关键字 一 l 硬球流体：密度泛函理论( d f t 沁分子动力学f m d 并r o s e n f e l d 密度泛函理论 上燕交通大学碗士学位论文 h a r d s p h e r e l i q u i d b e t w e e nt w ow a l l s a b s t r a c t i n19 6 4 ，h o h e n b e r ga n dk o h ng a v et w os t r i c tt h e o r e m sf o r q u a n t u m m e c h a n i c a lg r a n ds t a t e sw h i c ha r et h ef u n d a m e n to fd e n s i t y f u n c t i o n a lo ft h e o r y ( d f t ) 1 9 6 6 ，m e r m i ng e n e r a l i z et h e s et w o t h e o r e m st of i n i t et e m p e r a t u r e ，a c c o r d i n gt od f t ，t h ef r e ee n e r g yo f t h em a n y b o d ys y s t e mi saf u n c t i o n a lo fd e n s i t ya n dt h ed e n s i t yf i x e s t h ee x t e r n a l p o t e n t i a l r o s e n f e l dg a v e af u n c t i o n a lw h i c hh a sb e e n w i l d l y u s e d w eu s er o s 蛐f e l d d e n s i t yf u n c t i o n a lt h e o r yt og i v et h ed e n s i t y p r o f i l eo fh a r d s p h e r el i q u i db e t w e e nt w oh a r dw a i l s t h er e s u l tw a s c o m p a r e dw i t ht h a tg i v e nb ym o l e c u l a rd y n a m i c t h ec o m p a r es h o w u st h a tw h i l et h ed i s t a n c eb e t w e e nt w ow a l l si st o os h o r t ，t h et w o r e s u l t si s q u i t ed i f f e r e n t 0 t h e r w i s 奠t h er e s e l l t so f r o s g n f e l dd e n s i t y f u n c t i o n a lt h e o r ya n dm da r ew e l lm a t c h e d i nt h i s p a p e rw em s ti n t r o d u c et h ed e n s i t yf u n c t i o n a lt h e o r y t h e nw eg i v et h er o s e n f e l dd e n s i t yf u n c t i o n a lt h e o r ya n dt h ed e t a i l o fo u rw o r k ，a n df m a l l yi n t r o d u c et h ec a l c u l a t i o nm e t h o dw eu s e di n o u rr e s e a r c hw 0 r k 上j 每交通大学硕士学位敝 k e yw o r d s h a r d - s p h e r el i q u i d ，d e n s i t yf u n c t i o n a lt h e o r y ( d f t ) ，m o l e c u l a r d y n a m i c ( m d ) ，r o s e n f e l dd e n s i t yf u n c t i o n a lt h e o r y 上海交叠大学硕士学位论文 第一节引言 第一章密度泛函理论 1 9 6 4 年，h o h e n b e r g 和k o h n i 。j 对量子多体问题的基态证明了两个严格定理，建立了密度泛函 理论d e n s i t yf u n c t i o n a lt h e o r y ，以下简称d f t ) 。这两个定理指出量子多体系统的基态能量 是粒子密度的泛函，面密度和外势一一对应。d f t 髓后在各个方向都有推广，已经成为处理固体 电子态( 能带理论) 的理论基础【2 j 。1 9 6 6 年m e r m m p l 把它推广到有限温度，指出多体系统的自 由能是密度的泛函，且密度与外势一一对应。过去几十年里d f f 被大量地应用于罪均匀经典流 体的研究，取得了很多成粜，加滠了人们对非均匀经典流= 体酶理解h ”。同时，这一方法也用于 固液相变的研究 8 1 ，非均匀量子液体问题的研究b ”】，如磁约束玻色爱因斯坦气体的凝聚【“1 磁 约束简并赞米气体1 1 2 1 等。 密度泛函于臣式的梅造是应用d f t 理：论舶盏键只要得到了系统的自电能密度泛函舶形式， 就可以求出所有其他热力学量，如自由能【i 、化学势、直接关联函数【1 4 1 等。对于一般的多体系统， 常用的构造自由能泛函的方法有局域密度泛函近似( l d a ) 9 t 1 0 1 加权密度平均近似m 嘲等。对于 每一种具体的相互作用需要梅造相应舶密度泛函 瞄。一种简单胸相互作用系统是硬球系统。对 于这种系统yr o s e n f e l d 给出了加权密度平均形式的自由能泛函形式【l ”，通过这一泛函形式可 咀得到限制在各种不同几何体中的硬球流体的复杂的热力学性质呻j 卅，其结果与m o n t ec a r l o 或 分子动力学模拟的结果非常一致，得到广泛的应用。对于萁它形式的相互作用。也有许多很好的 泛函形式。如l e n n a r d - j o n e s 流体1 2 0 1 、带电流体睇l 】等。 本章的后面部分我们将给出密度泛函的理论基础然后在第二章中介绍我们所用到的泛函形 式及计算方法在第三章中将给出详细的i 十算结果，在附最部分中我们将介绍一下编程计算中用 到的共轭梯度方法瞄1 和分子动力学方法田】。 第二节密度泛函理论 假设n 个粒子( 举标为亏，动置为a ，质胃为珊) ，粒子之问没有相瓦作用- 处于体积y 的 连续分布外势西中。则该系统的哈密顿量为： ，：羔i 嘉+ 羔烈元) ：+ l 万( 兰子扩一i ) ) 烈尹) ”。善缶+ 善烈元) = + l 万( 善艿扩一i ) ) 烈尹) = 巧+ l 啾硼尹) ( 1 - 1 ) 上海交叠大学碟士学位论文 j v 其中z 是动能，氏( 尹) = 占( f 一亏) = 声( 芦) 是单粒子密度算符( o n e - p a r t i c l e d e n s i t yo p e * a t o r ) ， 1 = 1 在温度不变的备件下平均单粒子密度为： 以尹) = ( 声( 尹) = 行( p ( 尹) g 肌) t r ( e 一用一) = p e 一岸即l i v ( a y l v ) p 一朋3 ( 1 - 2 ) 该系统的h e l m h o l t z 自由能为： f = 一k 8 t i n t r ( e 一掣“1 = 一a r k 。，l n ( e p a 3 ) ( 方矿) p 俐】( 1 - 3 ) 其中，a 是热波长，p = ，v ，吒是珲想气体自由能。利用( 1 - 2 ) 式可以得到： r l 击：烈芦) 扩) = b r l 西钗芦) 【l n p 妒) a 一1 1 = 厶【厦尹) 】 ( 1 4 ) 可见h e l m h o l t z 自由能是与单粒子孵度有关的泛函秆极限情拥- 烈尹) _ + 0 f ， 凡。k 7 1 却【1 n p a 3 1 】：蛾( ，) 篇( p ) 是每个粒子在该条件下的自由能。 下面再来考虑经典相百作用下的非均匀系统。首先，哈密顿胃写成； ，= 巧+ 师( ，渺( ，) + 矾亏，不；矿) = + + l 师( f 渺( 尹) ( 1 - 5 ) 其中c ，j 、t 是相百作用能。 南度泛函理论以 珊( h o h e n b e r g - k o h n - m o r n u n ) 定理为基础。现存给出这一定理的表述和证 明。 对于( ! - 5 ) 式描写的系统，在室温条件下，其自由能由配分函数的对数给出 其中 f = - k 口t i n 甄( r ) ( 1 - 6 ) q n = t r e 书 2 上海交叠大学焉士学位论文 这罩 m 一田= 斧南赢氟帝。帝。e 唧 是对相卒间的积分。此外，自由能也可以写成对i t 女l j 分布的平均 f = 丹( 旧+ t i n ) ) ( 1 - 7 ) 其中w o i ie - q ，是正则分布函数。现在考虑一分布函数的泛函 冗明= t r ( w ( h + k t i n w ) ) ( 1 - 8 ) 其中是归一化的分布函数，t r w = 1 。当= w o 时，f = f 为系统的真实自由能，而当 w 时，有 冗阶瓢】= k b t n ( 眦势 。 ( 1 。) 当相互作用给定时，完全由外势决定，从而自由能及系统的平均密度 p o ( 芦) = 什( a ( 芦) ) ( 1 - 1 0 ) 也由外势决定。 另一方面，如果假定系统的平衡密度p 0 扩) 给定，则对应的外势也唯一确定。否则，假定 p 。扩) 对应于两个不同的外势珐( 芦) 和疵妒) ，则有 暇= e - 朋t q l= p 胛2 驮 r l 1 1 ) 这罩 h 。= + + 硝尹) ( f ) ：= 如十+ j 彩( 尹地( 尹) 可以证明 e = 开( 厅1 ) + rr r ( w , 如嘎) r r ( h 1 ) + k tn ( i n ) = e + 纸( 尹x 氟扩) 一妒：扩) ) ( 1 - 1 2 ) 3 上海交通大学礞士学位论文 此外还可以证明 b 鼻一却。扩脚。( 尹) 一：( 产) ) ( 1 - 1 3 ) 与前式矛盾。因此必然有破扩) = 以( 芦) 。所以p k 和外势烈尹) 就有一一对应的关系，从而m o 也 可以表示为p 。的泛函 p o 】。对于仟一给定的台珲的烈芦) ，总存在唯一的外势( 尹) 与其对应 由此可以得到一个【p 】。于是，我们可以构造密度的泛函 f p 】= t r ( w o 【p 】( + u j 、r + k t i n 硝) ( 1 1 4 ) 对于给定的外势烈尹) ，定义自由能泛函 h 用= l i p 】+ i 却( 芦( f ) ( 1 1 5 ) 则当p ( 尸) 取平衡密度时，- r 纠达到极小，即平衡密度由如下条件极值决定 研砌p l u = o ， j p ( 芦) 布= n 。 以上讨论和结论就构成了所谓唧定理的内容。 配分函数可以写成 q ：t r e 一田：t r ( e - p ( h “w ) 8 一町- 柙仰州7 ) ( i - 1 6 ) 自f f a $ i f = - k e r l n q 是外势的泛函。当外势发生变化捌( f ) ，( 尹) + 扩) + 聊( f ) ，则得到 h e l m h o 1 t z 自由能的变化： 铲= l 方 耐( 芦) 可以马上得到： 黑：反尹) 却铲) 一 对上式取二阶变分， 印( 产) 扩f d 烈，) 职芦) 6 烈尹r ) = 一卢 = i 尹) p ( f ) 一一( 尹一尹) 从f ) ( i - 1 9 ) 对关联函数g 2 ( f ，一) 和总关联函数厅2 ( 尹，) 由下式定义 = p ( 尹) p ( 尹) g 2 ( 尹，尸) = p ( 芦) p ( 芦) ( 1 + 矗往( 产，芦) ) ( 1 - 2 0 ) 因此可以得到 毒黯一肌胛渺2 协，扩。晰) ) ( 1 - 2 1 ) 由上面的式子可以得到， 布p 2 亿一) = 户伊以痧以尸) 矗啪( 尹，一) + p ( 尹) = = ( 一d p ( 产) ( 1 - 2 2 ) 从上面的式子中我们也可以得到一个对于总关联函数矗2 ( 尹，尸) 的重要约束条件： l p ( 芦) 2 ( f ，尹) = 一1 第三节直接关联函数和o r n s t e i n - z e r n i k e 方程 多体的密度一密度关联函数，或者说总关联函数是对f = h 妒伊) 】泛函求导得到。 给定的系统，若能找到 p 1 的泛函形式，由觚r 定理可以得到平衡条件 = 烈芦) - i - 旁如( 芦) ( 1 - 2 3 ) 对于一个 ( 1 2 4 ) 上式中以l a g r a n g e 园子的形式出现，侗它其实是化学势。 在哈密顿囊中，动曹b 和相互作并豫是相加的，从而在经典统计下二者葡贡献可以分别计算， 即 = k + 。 其中，厶是存没有相百作用情况下的自由能，即 双芦) = ( 1 ，) e x p 觑一妒( 芦) ) 厶= 丘却( 产) ( h 烈芦) a 3 一1 ) 上海交叠大学焉士学位论文 与密度p ( f ) 有关。厶是剩余自由能( e x c e s sf r e ee n e r g y ) 。把，代入( 1 2 4 ) 中，可以得到 或者 = 妒( 尹) + 七日，l n p ( 芦) ? + t o 印( 尹) 烈产) = pe x p ( 一j 日烈f ) + c 1 ( 砷一矗”) ( 1 - 2 5 ) 其中一1 ( 尹) = 一6 净：d 以芦) 。如果我们定义 则可把c 1 ( 尹) 改写为： 当西- - - h 0 时，有： 享( 尹) = ( 1 a 3 ) e x - p f l ( t 一( 芦) c 1 ( 尹) = i n p ( f ) ( 产) ( 1 - 2 6 ) 哿= i n ( p z ) ；= = ( 1 ) e x p ( 芦如) c 1 ( f ) 是由关于密度的泛函厶泛函求导而来。它是一阶直接关联函数。由于( f ) 是外势，因 此在单粒子情况下，用c 1 ( f ) 来表示系统内产生的相互作用。 必须注意到，c 1 ( 产) 本身也是一个关于p 扩) 的泛函。因此， p ( 尹) = pe x p 一卢铲( 尹) + c ( 芦；【p ( 尹) 】) 一c5 i ( 1 - 2 7 ) 由上式可以看到，当烈芦) - - 0 可以得到两个完全不同的解。若取c 1 ( 尹) - + c g ，可以得到均匀 流体的结果；但还可吼存存对称件被破坏的解 p ( 芦) = pe x p c 1 ( t 【p ( f ) 】) 一c o o ( 1 - 2 8 ) 上式具有平均场的形式。这个方程的稳定性是由c 1 ( 【p ( 芦) 】) 来决定。 下面我们再来讨论一下o n l s t e i n z e r n i k e 方程。对c ( f ) 进行泛函求导 c 2 ( 芦，尸) g 岔1 1 ( f ) ，t ；d ( 芦) = ( 占t ；d ( 芦”) ) l n p ( 芦) ，三( 尹) = ( r t p ( f ) ) 万( 产一i ”) 一弛芎( 对本，伊1 ) ( 1 七g ) 6 c 上海交叠大学焉士学位论文 利用 占( 产一r ) = 占i n 芋( 芦) ，d l n 手( 芦) ：f t 8 i n f ( f ) 墅贮2 i v 印伊) 8 1 n f ( f ) 根据( 1 - 2 1 ) 式可以得到： 毒堑錾羔：以尸) 烈一) ( 句( ，t ，) + 占( ，t 一尸) p ( 一- ) 占i n 三f 一1 、 、7、 则由( 1 2 9 ) 、( 1 - 3 0 ) 、( i - 3 i ) 可以得到： ( l - 3 0 ) ( 1 - 3 1 ) 2 1 扩，) = c t 封( e ，) + i l ：矿k 2 ( f ，冲2 扩”，f ) ( 1 3 2 ) 上式给出了两体总关联函数i 2 伊，一) 和两体直接关联函数c 2 ( 芦，一) 之问的关系。由于均匀流 体具有转动不变件和平移不变忡我们可以把卜式写成等价的形式+ 聪2 ( 尹) = c 2 ( 尹) + 以布。c 乳尹一尹磁2 ( 一) ( i - 3 3 ) 这就是0 r n s t e i n z e r n i k e 方程( 简称0 z 方程) 。二体直接关联函数c 2 ( 尹，一) 是关于均匀流体的 o z 方程。上式是一个卷积形式，可以对它作傅立叶变换得到， 毋( 七) = 搋正7 毋( 芦) = 0 ：x l 一1 s 2 ( 量) ) 其中； s 仕) 一l = l 托正钟( 尹) = 执正p ( 9 5 2 ) ( 尹) 一1 ) s 2 ) ( 简写作s 他) ) 是可测量的量，它是流体的静态结构因子，其定义如下： s ( 孟) = ( 1 ) 一以。 ( 1 3 4 ) 静态结构因子满足求和规则。1 ： 嘞s ( 七) = 砖e t i c r 其中f ，是绝热压缩系数。 7 上海交叠大学焉士学位论文 多体的关联函数可以通过泛函求导来表示成关于c 2 ( 芦，一；陋2 】) 的函数。比如 c 。( f ，尹。，尹。) = 8 c 2 ( 芦，产) 8 p ( 尹。) 显然它是一个 2 ( f ，一) 的泛函。因此 = 一占3 舟乞。西“尹) 每， 3 盯时，r o s e n f e l d 的密度泛函理论已经可以给出很好的结果从袭( 3 1 ) 可队看到此时寐 1 3 上海交蠢大学焉士学位论文 度的偏差小于1 。系统的很多整体性质来源于对密度分布的秘分，因此，其误差应该更小。在 得到密度分布的基础上我们可以进一步来研究系统的相变多分置体系的熵力等一系列性质。 已 o 6 4 2 0 口m d d ，t l 、v 、二 0234 z 6 图( 3 - 1 ) 墙的距离为= 8 0 ( 盯为小球的直径) ，密度为n ，= 0 6 8 8 时系统的 密度分布曲线。实线为密度泛函的结果，点是分子动力学模拟的结果。曲线和点 几乎重合，表明这两种方法的结果符合很好。 f i g ( 3 - 1 ) t h ed e n s i t yp r o f i l e so f h a r d - s p h d ef l u i db e t w e e nt w oh a r dw a l l sw i t ht h e w a l ls e p a r a t i o nl = 8 a ( 盯i s t h e d i a m e t e r o f t h e h a r d s p h e r e ) t h e l l t l m b t t d e m i t y i sp 020 6 8 8s o h dl i n ea n dt h ed o t sc 0 卿o dt ot h er e s u l t so b t a i n e db ym o t u l a r d y n a m i c sa n db yd e n s i t yf u n c t i o n a lt h e o r y ( d f t ) r e s p e c t i v e l y , t h es i m u l a t i o nd a t aa r e a l m o s tl o c a t c d0 1 1t h el i n eo f d f tr e s u l t s 8 6 蔷一 2 z a 同围( 3 - 2 ) 同图1 3 - 1 ) ，l = 5 盯 f i g ( 3 - 2 ) s a m ea sf i g 1 3 1 ) ，l = 5 仃 上海交迁大学鼋士学位论文 。 导4 & 。 z ，g 同图( 3 - 3 ) 同图( 3 - 1 ) ，l = 4 0 和l = 4 5 0 - f i g ( 3 3 ) s a m ea sf i g ( 3 - 1 ) ，l = 4 0 a n dl = 4 5 0 6 4 z ，a 同圈( 3 - 4 ) 同图( 3 - 1 ) ，l = 3 0 和l = 3 5 0 f i g ( 3 - 4 ) s a n ea sf i g ( 3 一1 ) ，l = 3 0 a n dl = 3 ，5 0 o d ，包 上海交叠大学磺士学位论文 i 吞 z a 同图( 3 - 5 ) 同图( 3 - 1 ) ，l = 2 9 0 f i g ( 3 - 5 ) s a m ea sf i g ( 3 1 ) ， 上= 2 9 0 0 8 舍叫 04 263 6465 6 l ，6 同围( 3 - 5 ) 密度风= 0 6 8 8 1 耐，墙的距离不同的系统中，所得到的第一个密度 极小处的密度与培距离的关系 f i g ( 3 - 5 ) t h e d 锄s i i y a t t l m f i r s t d e n s i t y m i n i m a l f l l n 击o n o f f i c e w a l ls e p a r a t i o n l t h en u m b e r d e n s i t y 岛。o 6 8 8 0 8 6 4 2 l ，g 同围( 3 7 ) 密度风= o 6 8 8 1 1 目 ，培的距离不同的系统中，所得到的第一个密度 极小的位置与墙距离l 的关系曲线 f i 垂( 3 - 7 ) t h ep o s i t i o na tt h ef 缸d e 啦畸埘n i m a i 鸥如毗i o no f t l ”w a l ls c p a x a t i o n t h e n u m b e r d e n s i t y 岛= 0 6 8 8 已 乏 24 56 l ，a 同图( 3 _ 8 ) 密度风= o 6 8 8 时，墙的距离不同的系统中，两堵墙中心处( ：= l 2 ) 密度大小与墙距离l 的关系 f i g ( 3 - 8 ) t h ed e n s i t y a tm i d d l eo f t h et w o w a l l s ( 二= l 2 ) a sf i m c d o no f t h ew a l l s e p a r a t i o nlt h en u m b e rd e n s i t y 风。o6 8 8 1 7 上海交蕊大学碗士毕业论文 参考文献 【1 】h o h e n b e r gp ，k o h n 、，p h v s r e vb ，1 9 6 4 1 3 6 ：8 6 4 2 】谢希德，陆栋，固体能带理论，复旦大学出版社，上海1 9 9 8 f 3 1m o r m i n n d ，i h y s 胁a ，1 9 6 5 ，1 3 7 ：l1 4 t 【4 ll o w e , n h l ，p h y s r 印，1 9 9 4 , 2 3 7 ：2 4 9 【5 l e v a n s r ，f u n d a m e n t a l s o f l n h o m o g e n e o u s f l u i d s , d e k k e r , n e w y o r k , 1 9 9 2 【6 】e v a n sr ，a d v 脚s ，1 9 7 9 ，2 8 ：1 4 3 【7 】g r o s sek u d r e i z l e rr m ，d e n s i t yf u n c t i o n a lt h e 0 1 , p l e n u mp r o s s , n e wy o r k , 1 9 9 5 【8 】r a m a k r i s h n a ntv ，y u s s o u f f m ，1 , h y s r e v b ，1 9 7 9 ，1 9 ：2 7 7 5 ( 9 】o l i v aj ，1 , h y s r e v b ，1 9 8 9 a , 3 9 ：4 1 9 7 1 0 o l i v aj p h y s r o yb ，19 8 9 b , 3 9 ：4 2 0 4 【l l l c h o u t t ，y h n g cn ，y h lh ，p t 邺r e va ，1 9 9 7 ，5 5 ：1 1 7 9 1 2 1c h e r tyj ，m a h r ，c h i n e s el h y s l e t t s ，2 0 0 1 ，1 8 ：3 3 【1 3 】t a nz ，m s a e o n iu mb ，s w o l e v a n , g u b b i n ske ，j c h e m 1 , h y s ，1 9 8 9 ，9 0 ：3 7 0 4 【1 4 ls o k o l o w s k is ，f i s c h e r j ，m 0 1 p l a y s ，1 9 9 0 ，7 0 ：1 0 9 7 r 1 5 1 d e n t o n a 1 l ，a s h e r o t t nw ，e h y sr “a ，1 9 9 1 ，4 4 ：8 2 4 2 【1 6 1k h e m a a s h e r o t tn w ，l h y s 叭e 19 9 9 ，6 0 ：2 8 7 5 【1 7 】n o w o r y t a jp ，h e n d e r s o n d ，s o k o l o w s k is ，c h a r t k y ，m 0 1 p l a y s ，1 9 9 8 ，9 5 ：4 1 5 f 1 8 】e v a n sr ，j p l a y s c o n d e n s e dm 瓶1 9 9 0 , 2 ：8 9 8 9 【1 9 】r o t hi l ，d i e t r i c hs ，p h y s r e v e ，2 0 0 0 , 6 2 ：6 9 2 6 【2 0 lr o s e n f e l dy ，l h y s r e v l e t t ，1 9 8 9 ，6 3 ：9 8 0 2 1 1 r o s e n f e l d y ，j c h e m 1 , h y s ，1 9 9 3 ，9 8 ：8 1 2 6 【2 2 】马红孺，计算物理讲义，上海交通大学，2 0 0 0 【2 3i 马红孺，研究生计算物理课程资料，上海交通大学物理系 秦克诚译理论物理学中的计算机模拟方法，北京大学出版社 【2 4 】h s m c n j ，- p ，m c d o n a l d ir ，t h e o r y o f s l m p e l i q u i d s a c a d e m a e p r e s s ，l o n d o n , 1 9 8 6 【2 5 lk i e r l i ke ，r o s i n b e r gml ，p h y sr e v a ，1 9 9 0 , 4 2 ：3 3 8 2 f 2 6 1a l l e nm ，同d e s l e yd ，c o m p u t e r s i m u l a t i o no y l i q u j d s ，c l a r e n d d ap r e s s , o x f o r d ，1 9 8 7 【2 7 】k i l n s c ，j o n e s g l ，l l a y , s r c v a ，1 9 9 0 , 4 1 ：2 2 2 2 h e n d e r s o n & l ，a s h e r o f tn w ，l , h y sr e v a ，1 9 7 6 ，l3 ：8 5 9 f 2 8 】p c r e t t sj k ，l h y s r e v l e t t ，1 9 6 2 , 8 ：4 6 2 p t c e u jk t h e e q u i l i b r i u mt h e o n ，o c l a s s i c a l f l u i d s 1 9 6 4 【2 9 lv a n l e e u v e nj mj ，c , r o e n e v e l dj ，j d eb o e r p h y s i c s ，1 9 5 9 ，2 5 ：7 9 2 3 0 】k nd m ，l a i r db b ，p r a y s m a ，1 9 0 0 ，4 2 ：4 8 0 6 【3 i 】m e i s t o r t f ，k r o l l km ，l h y s r e v a1 9 8 5 ，3 1 ：4 0 5 5 1 8 【3 2 】n o r d h o l ms ，j o n s o nm ，r r e a s i e bc ，n u s tj c h e m ，1 9 8 0 , 3 3 ：2 13 9 j o n s o n m ，n o r d h o l m s jc h e mp h y s ，1 9 8 3 。7 9 ：4 4 3 i 3 3 1t a r a z o n ap ，m 0 1 p h y s ，t 9 8 4 ，5 2 ：8 4 7 【3 4 】c u r t i nw ，a s h c r o t nw ，p h y s r e v a ，1 9 8 5 3 2 ：2 9 0 9 【3 5 】d e n t o na k ，a s h c r o l rnw ，啦s r e va ，1 9 8 9 ，3 9 ：4 7 0 9 【3 6 l i k o sc ，a s h c r o i t n w ，p h y s r e vl e t t ，1 9 9 2 ，6 9 ：3 1 6 1 9 上海交蕾大学磺士孥位论文 附录一共轭梯度方法介绍 对于一个含有雎个自变量的目标函教f ( 墨，也，屯) ，我们的目的是要求出一组 ( i ，蔓，z ) 使同标函数达到极小。为了讨论问题的方便，下而我们用x 代表，l 维空间的一个 点x ，使得目标函数f ( x ) 昂小。 对十高维晟优化问题有报多求解方法。但其基本总路是檀同的，就是从一点】c 0 出发，按j ! 【某 种j 妻i 定的方向p o ，用一维最优化的方法找到极小点x 1 ，继续这一过程直到达到最优点。若第f 步 到达码，设选定方向只，则新豹k + i 可如下求得。实际上，我们只要找一参数，i 一使得 f ( x , ：+ t l p j ) = r a i n f ( x i + t p 。) j ( 4 1 ) 则x i + i = x i + t p 。便为所求。因此，高维擐优化的不同方法实际上是选择不同的一维寻找方向。 一个育观的推测就是存每个点选择该点目标函教下降晟快的方向 a = 一订( x j ) ( 4 - 2 ) 作为寻找方向，遮一方法称为最陡下降法。直观的想象可能会认为最陡下降方向是一种坤想的寻 找方向，但事实并非如此，最陡下降方向只是表示在k 附近的下降性质，而对整个最优化过程来 说，我们寻找的是全局的最优方向。 塌陡下降方向并不是理想的下降方向，为了改进最优化的收敛速度，需要寻求其他方向。为 了判断一个寻找方向的好坏，我们应该青一个标准，注意到一般函数在最小点附近近似于二次函 数，因此如果一个方法对二次函数比较有效则可望其对般函数也比较有效，如果一个方法对 二次函数的效果都不好的话根难期望它对一般两数会有很好的效果。计算实践也证明了这一论 断。下面，我们就从二次函数寻求一种比较有效的下降方向。为此，我们先给出几个数学定理： 定理一：若维欧几里德空间中的向量口与个线性独立向量 p l ，p 2 ，p 都正交，则向量4 必为0 。 这一定理的结论是显然的。 定义：a 共斩向昔。假设一是个的对称讵定矩阵，g 是两个维向晕若， p7却=0(4-3) 则称向量p 和向量g 百为a 共辄或者可= 为a 一交。 上海交通大学磺士学位论文 定理二：若a 为n x n 对称i f 定耋巨阵，p l ，p 2 ，p 为a 共辆的n 维非零向着，则此向量 组必为线件独立。 证明：设向量组p l ，p 2 ，p 之间存在线性关系 口l p l + 口2 岛+ 。+ 口p = 0 对f = 1 ，2 ，n ，用p f a 左乘上式得到 因为p 。0 ，a 诈定，所以 从而必有 q a p , = 0 试点p 。 0 口，= 0 ，f = 1 , 2 ，一，n 因此向量组p i ，p 2 ，p 线性独立。 有了以上的数学准备，我们来考虑二次函数的最优化问题，我们的目的是找一种下降方向 能够在有限步达到二次函数的极小点。考虑二次元函数 f ( x ) = 口+ 6 7 x + 1 2 x t 且x 其中爿为一正定矩阵。上式的梯度为 g ( x ) = b + a x ( 4 - 4 ) ( 舢5 ) 记掏点的梯度为g 。= g ( ) 【i ) ，从x c , 点出发，对于第一个寻找方向- 除梯度外投有别的信息，我 们就取p o = - g o ，沿此方向找到极小点筠并可求得g l ，因为x l 是沿p o 方向的极小点，所以 g l 上p o ，从而g l 上g o 。选新的寻找方向 并要求p l 与几共耨则应有 p i = - g l c b g o ( - g i 一口。昙0 ) 7 4 = 0 ( 4 - 6 ) ( 4 - 7 ) 2 1 上海交叠大学焉士学位论文 由于 1 1 。x 0 + f o p o g t g o = a ( x l 一) 【0 ) = ，0 4 一般若第j 步的寻找方向为只，则n n - j i l t 因此，方捍( 4 7 ) 成为 由此解得 g i “一g f = a ( x i + l 一玛) = 如 ( 一g l a o g o ) 7 ( g l 9 0 ) = 0 口。：墨粤 蕊g o ( 4 8 ) ( 4 9 ) ( 4 1 0 ) 这罩假定g o 0 t 否则x 0 就是所求的极值点。把带入( 4 - 6 ) ，便得到a 沿a 求出极小点x 2 并算出9 2 。由于风和a 为彳共轭，因此由( 4 - 8 ) 式 g j ( 9 2 一9 1 ) = 0 考虑到g o 与g l 丁f 交，从而g 。与9 2 也f 交。y - g = 上p l ，则有 g r ( _ g l 一口o g o ) = 0 l i p 9 2 上g l g o ，g l 和9 2 构成一个萨交矢量组，我们可咀存这个i f 交矢量编构成的三维中间寻 求与p o a 均为a 共轭的寻找方向p 2 。令 p 2 = - 9 2 一q 9 1 一属g o 要求，2 分别与风和p 1 为a 菇轭，故必有 ( - 9 2 一口i g i 一风g o ) 7 ( g i - g o ) = 0 ( - 9 2 一口l g l 一屁g o ) 7 ( 9 2 一g i ) = 0 ( 4 1 1 ) ( 4 1 2 ) ( 舢1 3 ) 上海交叠大学磺士学位论文 由此可得 于是 ：窜生 g ig l 风= a o a l p 2 = - g ：4 - q p l 现存我们假定已经求得了p l ，p 2 ，p j ，其中 口o = 一g o p i = 一g j + 口i l n 一- ， 口“= g r g 。一t i g i - l ， f = l 2 ，j ( 舡1 4 ) ( 4 1 5 ) 这罩假定g o ，9 1 ，t 毋都不为零日构成一个萨交系。沿乃求出x j + i 并算出g j + l = g ( x j 。j ) ， 设g + i 0 ，则仿前可证g o ，g i ，g + l 也构成手交系令 p j + i = 一g j + 1 一f l o g o 一一岛一i g ，_ 一a j g j ( 4 1 6 ) 并目要求乃“与n ，p 2 ，乃为爿共轭，即 得到 ( - g j + l j o g 。一- p , 一l g ，一i n o g j ) r ( g j “一g ) = o ( 4 - 1 7 ) 口j = 盛“g m g r j g j ， p j - 。a j a 卜t 届= 吩口j l 口l ， 参n = a j a j 4 a p o 上海交_ 薯大拳焉士学位沧文 且 一p o g 。一b , g l 一一局一i g 一l 一吩g j = 川，口i ( a o g o + 9 1 ) 一屈9 2 一一e ；j g j = 口，( a l p l 一9 2 ) 一口j g j =q|p(4-18) 印 p j “= 一g j 札+ a j p i ，a i = 或n g j n g j g j t 4 2 0 ) 根据归纳法，在( 4 - 2 0 ) 中令= 0 , 1 ，2 ，所得方向为a 共轭方向。上述方法最多只要步就 可找到函数( 4 4 ) 的极小点这是因为所有的g 互相正交，即g ；g = 0 对所有 j = o ，l ，2 ，一1 都成立，二筏们的问题存维卒闻，故必有g = 0 ，亦即x n 为函敬( 4 - 4 ) 的极小点。 这就是所谓的共轭梯度方法，理论上，对于二次函数，只要次迭代就可以达到极小点，但 实际计算时，由于舍入误差的影响，一般次迭代并不能达到极小点。对于般的目标函数，有 限次( n 次) 的迭代也一般不能达到极小点，因此实际计算时，从一初始点出发迭代次后， 以所得结聚为出发点，重新开始计算，直到收敛为止。 j 上海交遗大学礤士学桓论文 附录二分子动力学方法介绍 存这罩我们所讨论的系统都明确考虑其牟部自由度，不允许存存随机1 峰因素，讨论的出发点 是经典力学范围内的表述。我们是要计算这些系统的物理量。 1 分子动力学 分子动力学( m d ) 方法的出发点是物理系统的确定的微蹋描述。它是用运动方稗来计算系 统的性质，结果得到