Reverse Monte Carlo计算机模拟

D0第糝第期sn1001-053x.1996.光0紧科技大学学报 Vol.18 No.1 19962 Journal of University of Science and Technology Beijing Feb.199% Reverse Monte Carlo计算机模拟 郭春泰 北京科技大学物理化学系，北京100083 摘要介绍了Reverse Monte Carlo计算机模拟原理，及其在应用时应注意的问题，并且详细 地给出了应用举例 关键词计算机模拟/溶液结构，蒙特卡罗法 中图分类号0242.1 液体和溶液的结构通常用径向分布函数(RDF)来表征.获得RDF的常规方法是X一射 线、中子和电子衍射实验测定及分子的Metropolis Monte Carlo(MMC)及分子动力学 (MD)计算机模拟计算·综合实验与计算能获取较为详细的结构信息· 将计算得到的RDF与实验测得的RDF相比较.如果二者符合较好，由模拟计算获得 的精细结构及三维信息则具有使用、参考价值·.但上述方法存在两个问题：(1)建立一个符 合实际体系的势函数是极其困难的，甚至是不可能的；(2)由实验测得的S(⑨)在变换为g(r) 时，有些结构信息已经损失，致使计算与实验的比较在某种程度上失去真实性·最近，MC一 Greevy等人l发展了一种新的模拟方法-Reverse Monte Carlo(RMC)计算机模拟.此方法 在欧洲已得到较为广泛的应用2~41.本文简介RMC的基本原理，并给出应用举例. 1RMC原理 Metropolis Monte Carlo(MMC)模拟是以粒子间作用势能函数为基础，而Reverse Monte Cao(RMC)方法则是以由实验测得的结构因子SE(Q)为基础，所以克服了上述MMC的两个 何题，除此之外，RMC和MMC的模拟过程基本相同， 在RMC方法中，认为实验测得的结构因子S()只包含统计误差，且此误差服从正态 分布，计算得到的结构因子S(Q)与S(Q)之差为： e,=S(2)-S(2,) (1) 并且e,具有概率： P(e)=(2π)2c(g,》'exp(-e/2o(2:)) (2) 式中σ(Q)为正态分布的均方差·系统的总概率为： P-ire)-(aonreap(-g2nOy) (3) 1995-09一07收稿 第一作者男40岁教授

第 18 卷 第 1 期 北 京 科 技 大 学 学 报 1臾场 年 2 月 oJ u 了n a 1 o f U in v O sr iyt 0 f S a e n c e a dn eT hc n 0 fo gy eB ij in g V d 。 18 N 匕 f 功 。 1望场 R e ve sr e M o n et C a ir o 计算机模拟 郭春泰 北京科 技大学物理化学系 , 北京 10 ) 刀 3 摘要 介绍 了 R e说n 七 M o 业 ( 泊d匕计 算机 模拟 原理 . 及 其在 应用 时 应 注 意 的问题 , 并 且 详 细 地给 出了应用 举例 关锐词 计 算机模拟 / 溶液结构 , 蒙特卡罗 法 中圈分类号 0 24 2 . 1 液 体和 溶液 的结构 通常 用径 向分布 函 数 ( R D F )来 表征 . 获得 R D F 的常规 方法 是 X 一 射 线 、 中子 和 电子 衍 射实 验 测 定 及 分 子 的 M e it D p o ilS M o n et aC r fo ( M M C ) 及分 子 动 力 学 ( M D ) 计算 机模 拟计算 . 综合实 验与计算能 获取 较为 详细 的结构 信息 . 将计算得 到 的 R D F 与实验测 得 的 R D F 相 比较 . 如 果 二 者 符 合较 好 , 由模 拟 计算 获得 的精 细结 构及 三维信息则具 有使 用 、 参考 价值 . 但上述方 法存在两 个 问题 : ( l) 建 立一 个符 合实 际体 系的势 函 数是 极其 困难 的 , 甚至是 不 可能 的 ; ( 2) 由实验 测 得 的 s (Q ) 在 变 换 为 g( r ) 时 , 有 些 结 构信 息 已 经损 失 , 致使计算 与 实验 的 比较 在某种程度上失去真实性 . 最近 , M c 一 G吠竹 等 人「’ l 发展 了一 种新 的模拟 方 法 一 R e ve n 七 M o n et aC lor ( R M C ) 计算 机 模拟 . 此方 法 在 欧洲 已得 到较 为广泛 的应 用「卜 4 } . 本 文简介 R M C 的基本 原理 , 并 给 出应 用举例 . 1 R M C 原 理 M e t or po ilS M o net aC lor (M M C )模拟是 以 粒子 间作 用 势能 函数为基 础 , 而 eR ~ M o net aC lor ( R M C )方 法则是 以 由实验 测得的结构因子 S “ (Q ) 为基础 , 所以克服 了上述 M M C 的两 个 间题 , 除此之 外 , R M C 和 M M C 的模 拟过 程基本相 同 . 在 R 五左C 方 法 中 , 认 为实 验测得 的结 构 因 子 S “ (Q ) 只包 含统计 误 差 , 且此 误差 服 从正态 分布 . 计算得 到 的结 构 因子 CS (Q ) 与 梦(Q ) 之差 为 : e ` = S C (Q ` ) 一 S E (Q : ) ( l ) 并 且 e ` 具有 概 率 : p ( e ` ) = (( 2 二 ) ’ ` ’ 。 ( Q ` ) 一 ’ e x p ( 一 。矛/ Z a (Q ` ) ’ ) ( 2 ) 式 中 叮 Q ` ) 为 正 态分 布 的 均 方 差 . 系 统 的 总 概 率 为 : 用 用 p 一 (少( “ ` ) 一 恢厄玉七刃 ) ’ e x p ( 一 冬 . e 子 2 叮 (Q . ) ’ ( 3 ) 19 5 一 0 9 一 0 7 收 稿 第 一 作者 男 4 0 岁 教 授 DOI: 10. 13374 /j . issn1001 -053x. 1996. 01. 007

Vol.18 No.I 郭春泰：Reverse Monte Carlo计算机模拟 29. 式中m是S(⑨)在2空间的实验点数，并且， -(ie, (4) 为了使用S(⑨)来模拟体系的结构，希望建立一个原子统计系综，并且该系综服从上述 概率分布.将(2)式中指数项写成， x2=2(s(2)-s(2,》1a(2 (5) =l 则总概率Po∝e2.我们可立即看出x2/2相当于MMC中的U/KT(U为势函数)， RMC模拟计算的步骤如下： (1)建立具有N个粒子且服从周期性边界条件的原始构型.N个粒子在构型中可以是 随机分布的，也可以是按某种晶体结构或其它任意形式分布. (2)计算统计这个原始构型的径向分布函数： i0=。 (6) 式中n()为以某个粒子为中心，半径r到r+△r球壳中的粒子数，然后取以每个粒子为中 心时的平均值，p是粒子数密度，构型尺寸L应该足够地大，以使g>L2)=1.径向分布函数 g(r)只计算rx名，则该运动为概率接受(xp(-(x-x)/2),否则该运动被拒 绝· (7)重复步骤(5)和(6)，直到x2达到一个平衡值. 模拟计算结果应该是一个与实验测得的结构因子相符的三维构型·上述RMC模拟方法 适于单组元体系的数据· 对于k个组元体系 x=Σx=∑∑(S(2,)-S(2,》P1o.(e)P (9) 对于中子散射 S(Q)=>>C.C;b.bn(S.(Q)-1 (10)

丫b l . 18 N b . l 郭春泰: R e w 晓℃ M o ln 匕 O 川。 计 算机模拟 式 中 m 是 S “ (Q ) 在 Q 空 间的实 验点数 , 并且 , 万一 (众 · (Q ` ) ’ ` “ ( 4 ) 为 了使用 梦 (Q J 来模 拟体 系 的结 构 , 希 望 建 立 一 个原 子 统 计 系 综 , 并 且 该 系 综服 从 上 述 概率分布 . 将 ( 2) 式 中指数 项写成 , x , = 艺( s C (Q ` ) 一 s E (Q ` ) ) , / 。 (Q ( 5 ) 则 总 概率 尸 oc ` 水 . 我们 可立 即看 出 x ’ / 2相 当于 M M c 中的 U/ K 八u 为 势 函数 .) R M C 模 拟计 算 的步 骤如下 : ( l) 建立具有 N 个粒 子且服 从周 期性 边界条件 的原 始 构 型 . N 个粒 子 在 构型 中 可 以 是 随机分 布 的 , 也可 以是 按某 种晶体结构 或其 它任意 形式 分布 . (2) 计算 统计这 个原始构 型 的径 向分 布 函数: ( r ) , 、 _ n J L r ) 气r ) 一 下二了万丁一一一, 兮兀 r 一 。 r ’ P ( 6 ) 式 中 叮 ()r 为 以某 个粒子 为 中心 , 半径 ; 到 ; + △ ; 球壳 中的粒 子数 , 然后 取 以 每 个 粒子 为 中 心 时 的平 均值 , p 是粒 子数密度 , 构 型尺寸 L 应 该足够地大 , 以 使 g (r > 乙/ 2 ) = 1 . 径 向分布函数 抓)r 只计 算 r 瑞 , 则 该 运 动 为 概 率 接 受 (e x p ( 一 (嵘一 瑞 ) /2) , 否 则 该 运 动 被 拒 绝 . ( 7) 重复步 骤 ( 5) 和 ( 6) , 直 到 x Z 达到一 个平衡值 . 模 拟计算 结果 应该是 一个 与实验 测得 的结构 因子相符 的三 维构型 . 上述 R M C 模 拟 方 法 适于单组 元体 系 的数 据 . 对于 k 个 组元体 系 x , = 艺 x 卜 艺工( s 万(Q ` ) 一 s 奢( Q ` )) , / a * (Q ` ) , ( 9 ) 对 于 中子 散射 S f (Q ` ) 一 万落 C · C , 石 · 万 , k ( S · , (Q ` ) 一 ` ( 10 )

·30 北京科技大学学报 1996年No.1 对于X射线衍射 S(2,)=∑∑C,Cgi(S,(2,)-1 (1I) 对EXAFS数据 对EXAFS数据 5(Q.)=E4xp(g.o(r)-1)(..r)dr (12) 在式(10)和(II)中，C.,C:为在试样中x和B组元的浓度，bx、f和b陆是在试样中a和B 组元的相干散射因子，S为a-B对结构因子.在式(12)中下标由a取代，因为EXAFS 谱是在组元a的吸收限附近测定的，有(Q,)是组元B在距离r处对EXAFS谱的贡献. 9(r)是x一B对的偏径向分布函数.RMC还可用不同种类实验测得的散射因子组合来求偏 径向分布函数（例如，由X射线、电子、电子衍射实验组合，求偏组元的RDF) 2 液体半导体MgBi2的RMC模拟s1 由中子散射测得的S(Q)和由RMC模拟计算的S(Q)示于图1.由图1可见，实验值与 计算值符合甚好，在低Q处的小肩峰也被模拟出来，图2和图3是模拟计算得到的偏结构因 子S,(Q,)和偏径向分布函数9，).由此可见，在二元系情况下仅由一个总RDF,便可 751012.51517.520 Q,、10nn Qj×10nm 图1熔融g:的结构因子 图2熔融g围：的偏结构因子 0.40 3.5 035 030 2.5 025 0.20 1.5 015 0.10 0.5 005 000 810 14 16 r/x10-nm- 配位数 图3熔融g,Bi,的偏径向分布函数 图4熔融Mg,Bi,的配位数分布

北 京 科 技 大 学 学 报 1望天i 年 N o . 对于 X 射线 衍射 ( 1 1 ) 对 E X A F S 数据 S f (Q , ,一 万落认 马 五飞` 工` S · 。 (。 1 ” 一 ` 一 厂 _ 对 E ~ 数据 S犷( Q ! , 一 今 4 “ 户 { r “ ( “ · “ ( r ’ 一 ” 寿( Q ! , r ’ “ r ( 12 ) 在 式 (l 0) 和 (l 1) 中 , 叹 、 C , 为 在 试 样 中 “ 和 刀组 元 的 浓 度 , b 。 、 、 人又和 练 * 是在 试样中“ 和 刀 组 元 的 相 干 散 射 因 子 , 凡, 为 : 一 刀对 结 构 因 子 . 在 式 (l 2) 中下标 由 : 取 代 , 因 为 E x A F S 谱 是在 组元 戊 的 吸 收 限 附 近 测 定 的 , 石(Q , , )r 是 组 元 方 在 距 离 ; 处 对 E x A F s 谱 的 贡 献 . 物 ()r 是 , 一 刀对 的偏径 向分 布 函 数 . R M C 还 可用 不 同 种类 实验 测得 的散 射 因子组合来求偏 径 向分 布 函 数 (例如 , 由 X 射 线 、 电子 、 电子衍射 实验 组合 , 求 偏 组 元 的 R D F , ) . 2 液体半导体 M g 3 悦 2 的 R M C 模拟l ’ l 由中子散 射测 得 的 S( Q) 和 由 R M C 模 拟 计 算 的 S( Q) 示 于 图 1 . 由图 1 可 见 , 实 验 值 与 计 算值符 合甚 好 , 在低 Q 处的小肩 峰也 被模拟 出来 . 图2 和 图 3 是模 拟计 算得 到的偏 结构因 子 S : 刀 (仓 ) 和 偏径 向分 布 函 数 乳刀 ()r . 由此 可 见 , 在 二 元 系 情 况 下 仅 由 一 个 总 R D F , 便 可 一 0 2 一 0 4 } , 」 } ’ … { 一 、 勺 / } } 一 一 } ; 一 … 习 12 5 1 5 1 7 5 2 0 一 } 一 10-0 入 ō é门气` 图 I 熔融 、 娘 , 段 : 的结构 因子 图 2 熔融 凡娘 , 代 : 的偏结构因 子 一 : { _ 。 i …i1 . { 。 门 040352530105 …日 擎瓣 丢 }仁二刃 l 5 1 O 5 0 心 、 } , ,户尧言 8 1 0 1 2 14 1 6 r / x 1 0 一 ’ n ln 一 ` 4 配 位 数 图 3 熔融 M g 3 B i : 的偏径向分布函 数 图 4 熔融 M g , B i Z 的配位数分布

Vol.18 No.1 郭春泰：Reverse Monte Carlo计算机模拟 .31· 通过RMC,求出3个偏RDF,这对无法测定偏结构因子的体系是十分有益的.图2 表明MgBi2具有典型的离子化合物（例NaCl,ZnCl)结构形式.Mg为正离子，Bi为 负离子.在低Q处，S+-负向峰和S,+,S的正向峰的排列形式是由于正负离子耦合以 达到最大的离子屏蔽· 1.400r RMC的另外一个优点是，除可正常计算 Mg-BI-Mg 1200 配位数以处，还可给出近邻配位的分布，因此 会1.000 如何定义配位壳层便显得不是十分重要了·图 言0.800 4为MgBi,液体中的配位分布.由图可知，对 于ne-,在Mg周围有3~4个Bi时配位的可能 慧dipthe 0.400 性最大， 0.200 由于RMC可给出三维立体结构，因此可 0000L 统计计算出配位对的取向或化学键的键角·图 -1.000-0.600 -0.2000.200 0.6001.000 Cos0(°) 5是键Mg-Bi-Mg的角度分布，由图可知 图5熔融Mg2监2g脱2和MgB风.键角的分布 Mg-Bi的角度为60°. 3结束语 尽管RMC方法已经受到了广泛的重视，但它的多维结构信息的唯一性问题还需进一步 检验证实.由统计力学和实验证明，对于偶势(Pairwise potential))Φ(z)占绝对主导作用 的液体，RMC是正确的.但对于多体势占主导作用的液体，在应用RMC时应多加注 意，此时最好将RMC和MMC相比较， 参考文献 I McGreevy R L,Pusztai L.Reverse Monte Carlo Simulation.Mol Sim,1988(1):359-367 2 McGreey R L,Howe M A.Condensed Matter,The Structure of Molten LiCl.J Phs,1989(1):9957~9962 3 Keen D A.McGreevy R L.Structural Modelling of Glass Using RMC.Nature,1990,344:423~424 4 Howe M A,McGreevy R L.Determination of Three Body Correlations in Liquids by RMC.Phys Chem Liq,1991,24:1~12 5 Chuntai Guo.Ph.D.Thesis],University of Bristol,1994 Reverse Monte Carlo Computer Simulation Guo Chuntai ng,19902:2773 Department of Physical Chemistry.USTB,Beijing 100083,PRC ABDTRACT Reverse Monte Carlo Computer(RMC),simulation is introduced in defail,and a brief comment is given.Application of RMC is also introduced by way of an example. KEY WORDS computerized simulation/solution structure,monte carlo metlod

Vb l . 18 N O . 1 郭春泰 : R e w 晓祀 M o n t e C妞ir o 计算机模 拟 通 过 R M C , 求 出 3 个 偏 R D F , 这 对 无 法 测 定 偏 结 构 因 子 的 体 系 是 十 分 有 益 的 . 图 2 表 明 M g 3 iB Z 具 有 典 型 的 离 子 化 合 物 (例 N a C I , Z n C 1 2 ) 结 构 形 式 . M g 为 正 离 子 , iB 为 负 离 子 . 在 低 Q处 , S + 一 负 向 峰 和 S 、 、 , S 一 达 到 最 大 的 离 子 屏 蔽 . R M C 的另外 一 个 优点 是 , 除 可正 常计 算 配位数 以处 , 还 可给 出近邻 配位 的分布 . 因此 如 何定 义 配位壳层便 显得不 是十分 重要 了 . 图 4 为 M g 尹1 2液 体 中的配位 分布 . 由图可知 , 对 于 n gM 一 。 , 在 M g 周围有 3 一 4个 iB 时配 位 的可 能 性最 大 . 由于 R M C 可给 出三维立 体结 构 , 因 此 可 统计计算 出配位 对 的取 向或化 学键 的键角 . 图 5 是 键M g 一 iB 一 M g 的角度 分 布 . 由 图 可 知 M g 一 iB 的角度 为 60 0 . 的 正 向峰 的 排 列 形 式 是 由 于 正 负离 子藕合以 { g M g 一 B I 一 州 … 尹入 沁协沙州 一 一 { \ 一从 ( ( D )屯/) 团 co · g 一 1 0 00 一 0 . 600 一 0 2 0 0 0 2 0 0 0 6 0 1 . 因 0 oC s (口/( 。 ) ) 图 5 熔融 M忿妙 . 2) 、 M g J残 和 M g 3 . I ` M 键角的分布 3 结束 语 尽管 R M C 方法 已 经受到 了 广泛 的重视 , 但 它的多 维结 构信 息的唯 一性 问题还需 进 一 步 检验 证 实 . 由统计力 学和 实验证 明 , 对 于 偶 势 (aP i“ 喊哭 p o ent 血l) 中 (r . 2 ) 占绝 对 主 导 作 用 的液 体 , R M C 是 正 确 的 . 但 对 于 多 体 势 占 主 导 作 用 的液 体 , 在 应 用 R M C 时 应 多 加 注 意 , 此 时最 好将 R M C 和 M M C 相 比较 . 参 考 文 献 1 M d 3 代笼竹 R L , P 出川旧1 L . R e记眨祀 M o n te C妞r】0 Sir n 创巨山n . M o l S而 , 19 8 ( l ) : 3 59 一 3 6 7 Z M以3 代犯 y R L , H o 忱 M A . 〔b 以北m ed M a t et r , , 11℃ S tr u Ct 切 re o f M o l te n L IC I . J P比 , 1989 ( l ) :卯5 7 一 望拓2 3 K“ 泊 D A, M困 er 理 R L . Str u tC t皿 I M以允 ha 唱 of G al s sU ign R M C . 卜抽 t uer , l 望刹〕 , 少抖 : 4 2 3 一 4 24 4 H o 讹 M .A M 以3 代先 W R L . L兄 t e r 面n a it o n of T b 代兄 佳劝y 山讹h t i o 飞 in L 峋山ds 勿 R M C . P h那 C 比m qjL , l卯l , 24 : l 一 12 5 hC un at i G ou . 【hP . D . T l l ` is 」 , U in vers ity o f B出ot l , l卯4 R e ve sr e M o n t e C a r l o C o m P u t e r S im u l a ti o n G z止0 C h u n t a i ng , l男X2) : 27 7 3 eD p a r t m e n t o f p h y s i e a l C h e m i s t 即 , U s T B , B e ij i n g 10 0 0 5 3 , P R e A B D T R A C T R e v e rs e M o n t e C a r l o C o m P u t e r ( R M C ) , s im u l a t i o n 1 5 i n t r o d u e d i n d e fa il , a n d a b r i e f co m e n t 1 5 g i v e n . A P P lica t i o n o f R M C 1 5 a l s o i n t r o d u e d b y w a y o f a n e x a m P l e . K E Y W O R D S co m P u t e r i z e d s im u l a t i o n / s o l u t i o n s t ur e t u er , m o n t e c a r l o m e t l o d