流体熵相关性质的Monte Carle摸拟新方法

1、流体熵相关性质的Monte Carle摸拟新方法 韩 敏 (无机化学专业 MG0124006) 摘要提出了用于计算实际体系的熵性质的Monte Carle多级分子取样摸拟方法.应用这一方法,对硬球流体的化学势及Helmholtz自由能进行了计算,得到了初步的结果.计算化学势时不存在通常实验粒子方法所遇到的高密度问题.该方法特别适合规律性的系统研究,较之普通模拟方法要有效得多.模拟得到的硬球体系无限稀释组分的超额化学势与对比直径的关系,在相变区域为一条双凹曲线;无论是在相变区还是在单相区,Canahan-Starling公式对这一关系的描述均有较大偏差.关键词 Monte Carle 分子模拟

2、硬球 化学势 Helmholtz自由能 熵 计算机分子模拟是伴随着计算机的出现而发展起来的一门技术,在科学研究中的地位正日益提高,在一些领域如相平衡的研究中,常常需要计算如自由能熵化学势等与体系的配分函数 性质有关的量.这些性质称为熵相关性质.通过分子模拟, 体系的系踪平均性质如内能压力对配分函数直接作计算可以得到, 但要计算熵相关性质,目前倘有一定的难度.许多研究工作者为此作了大量的工作,概括起来有以下几种方法.第一方法是设计一个可逆途径,使所研究的状态与以知熵相关性质的状态相连,两个状态的自由能之差可由 (A/ V) N,T=-P (1)(A/T)/( (1/T)N,V=E (2) 积分求

3、得1,该方法非常直观,但需要进行多个分子模拟以求得可逆途径上的压强或能量,并且如果以可逆途径通过相变区域时,这种方法就不适用了.第二种方法是利用“伞形取样”2或“接收比”3等方法来求算实际体系与某一具有相似势能的参考体系的自由能之差.这种方法目前已广泛应用于非理想流体混合物性质的计算,但若参考体系与实际体系差别较大,仍需通过多个模拟,才能保证计算的精度.第三种方法是直接根据模拟得到的构型采用试验粒子的方法求算化学势4,5,该方法的一个缺陷是难以应用于高密度流体,虽然已有一些旨在克服高密度问题的努力,如Mezei提出的孔基巨正则系踪方法6,但都不能从根本上解决这一问题. 针对分子模拟中熵相关性质

4、计算量大的难题,提出采用多级取样模拟的方法,使得在单一分子模拟过程中完成所设计途径上各状态点的模拟,从而求得体系的熵相关性质.1.理论 在正则系踪中,氦姆霍兹自由能与配分函数有如下关系: A=-KTlnQn,v,t-(3)可以设置一个描述体系某一特征的参数l ,使得l =1时代表实际体系, l =0代表某一参考体系 ,若l在0,1区间取n个值l i(i=1,2.n)使l i=0, l i+1=1,任意相邻的参数l i和l i+1,所描述的两个体系自由能之差可以表示为:A(l i+1)-A(l i)=-KTlnQn,v,t(l i+1)/Qn,v,t(l i)-(4)将配分函数的比值表示成系踪平

5、均的形式,可以得到 A(1)-A(0)=-KTln<exp-b(u(l i+1)u(l i)> l i-(5)上面的<> l i表示在参数为li时进行模拟的系踪平均,由上式可以看出,实际体系与参考体系自由能之差可以通过l1 . l 2 . l 3 . l n等系踪平均计算得到。显然,通过n个分子来计算自由能之差是一项艰苦的工作.为减轻工作量,采用多级取样模拟的方法,让参数l在模拟过程中逐级变化，从而得到各个参数下体系的系踪平均.应用多级取样模拟方法,可以推算体系的化学势.自由能.混合物的超额热力学性质等.对化学势的计算,可以表示成等容下N+1个粒子体系与N个粒子体系的自

6、由能之差的计算,可通过变化N+1个粒子尺寸的多级取样模拟完成.对自由能,可采用所有粒子尺寸的多级取样模拟,得到实际体系与理想气体或某一已知自由能的参考流体的自由能之差.对于混合物的超额自由能,则可通过NPT正则系踪里的多级取样模拟进行计算.2.模拟过程本文对硬球流体的化学势和亥姆霍兹自由能进行了计算.模拟的粒子数为54,采用立方周期性边界条件,初始构型为面心立方晶格,在开始多级取样模拟前经2-100000构型达到平衡.下面对化学势和自由能的计算作简要的说明.2.1化学势的模拟对于硬球体系化学势的模拟,参数可选为一个粒子与其他粒子的直径比,即l =-d/ d0,采用多级取样法计算当该粒子直径为d

7、和0时体系自由能之差.该粒子直径为0时的粒子体系与等容下N-1个粒子体系的自由能之差可严格由 AN(0)-AN-1=Ukl-KTln(1-rpd03/6)-(6) 得到.这里, AN(0), AN-1分别表示中心粒子直径为 0的N粒子体系和N-1个离子体系的亥姆霍兹体系的自由能.因此,硬球流体 的超额化学势为: mex = m-mid=AN(1)-AN(0)-KTln(1-rpd03/6)-(7)再求算化学势的多级取样模拟中,本文人为地选定标号为1的粒子为中心粒子,并 逐级变化其直径.由于计算自由能之差的被积函数只于中心粒子的作用有关,为提高该函数系踪平均精度,采用欧围齐的优先取样法7,使体系

8、中各粒子 的选择移动几率和其与中心粒子的距离有关,本文定义的几率函数形式为: W(ri)=1/(ri2+d02)-(8)则粒子的移动几率为 WI=W(ri)/W(rj)-(9)其中为粒子与中心粒子的距离.对于硬球流体,计算自由能之差的被积函数可简化为: 1 minri>li+1d0 EXP-b(U (li+1)-U(li)= -(10) 0 minri/>li+1d0若定义为在参数的这级模拟中各粒子与中心粒子的最近距离大于的构型数可得(EXP-b(U(li+1)-U(li) li=NK(i+1)/NK(i)-(11)对于硬球体系,在每一级的取样模拟中,若参数为,理论上可以计算得到所

9、有的,则和所描述体系的自由能之差,可以利用所有的各级模拟所获得的结果,按式(1 2)计算得出,以提高统计精度.AN(lI+1)-AN(lI)=-KTlnNK(i+1)/NK(i)-(12)3. 2亥姆霍兹自由能的模拟硬球体系亥姆霍兹自由能的模拟与化学势的模拟不一致.不同之处为在这里代表第i级模拟中粒子之间的最近距离大于ljd0的构型数.另外,在模拟中无须采用owicki优先取样法.在l=0时体系的自由能即为理想气体的自由能,由此即可得到各个值时的自由能.自由能是一种广延性质,与体系大小成正比.因此,为保证一定的统计精度,体系越大,在间划分的级数就越多.若总的模拟构型一定,则各级的模拟构型就越少

10、,相应地就会造成大的误差涨落.为获得实用的模拟结果,需要合适的数据处理方法.本文采用三次平滑技术,得到了较光滑的曲线.本文模拟了对比密度r*=rd3在0-1.2范围内的亥姆霍兹自由能,d0即为r*=1.2时的粒子直径,在范围内取值,使2的取值间隔为0.001d02.在每级模拟中产生1k个构型,其中前1/3的构型用于平衡而舍去,剩余构型用于计算.由于l在变化范围中存在相变区域1,为保证模拟的顺利进行,本文使由1逐级减小至0的途径进行模拟.3.结果与讨论3.1 化学势的模拟化学势的 大小如何,总能选择适当的级数n以获得足够精确的结果 .图1表示了本文模拟的结果,与近似结果8进行 了比较,偏差较大.

11、如果增大模拟的粒子数和构型数,可使精度进一步提高.在本文所模拟 的 密度范围r*=0.851.00内,通常所用的粒子实验法已难以应用,而采用本法可得到较合理的结果。图1中在相变区(r*>0.944)随密度的增加，可以看出硬球流体化学势先增大后减小的趋势。Pade近似在该区域已不适用，故与模拟结果有较大偏差。 多级取样模拟的一个优越之处是可以同时求得所设计途径上各状态点的熵相关性质，以便对其作系统的研究。在求算图1中各密度下的化学势的过程中，同时可得到对应密度下半径为的无限稀释组分的化学势。图2表示了本文模拟的无限稀释组分超额化学势与对比直径的关系曲线，同时还给出了相关密度下由Carnah

12、anStarling公式9计算的结果。在相变区，模拟的结果表明，无限稀释组分超额化学势与对比直径的关系为一条具有双拐点的双凹曲线，这一特征是CarnahanStarling公式无法描述的，由此可知，根据这类曲线的形状可以鉴定体系是否分相。3. 2亥姆霍兹自由能图3是模拟得到的硬球流体的氦姆霍兹自由能与密度的关系以及由公式计算得到的结果。可以看出，在0.90的密度范围内多级取样模拟结果与公式的计算结果相差不大。若在多级取样模拟中增加各级构型数，就会得到精确的结果。图3还给出了Frankel计算得到的硬球面心立方晶格的自由能10，与本文模拟的结果相近。图2，3所得的结果在多级取样模拟中通过一次普通长度的模拟就可得到。因此，对有规律性系统的研究比普通模拟方法要有效得多，一方面是由于在多级取样模