gemc法模拟流体汽液相平衡

更容易得到流体的临界点。科学研究需要不断创新，HlouchaM等人就尝试着对四维空间的体系进行模拟[33]。采用GEMC方法作出体系的相图，并得出四维流体的临界参数。另外，MC方法还可以用于研究溶液的活度系数。金文正等人阐述了无限稀释活度系数的模拟方法[34,35]，建立了分子蜕变偶合参数积分MonteCarloNPT系综法，模拟无限稀释水溶液的化学势差和活度系数。

(5)分子模拟研究相界面

分子模拟不仅能测定流体主体内分子的近程有序(即所谓局部组成)，还能测定汽液界面、液液界面以及气固、液固界面的分子分布。对纯流体垂直于汽液界面的分子密度剖面分布及界面层厚度均可从分子模拟方法得出的双曲正切函数表示

[36]

。虽然也可以采用密度泛函理论与SAFT状态方程中的微观参数预测纯流体(含

链状分子和氢键缔合分子)的表面张力和液液不互溶体系的界面张力。但上述计算仅能获得链节数沿垂直界面的分布，不能得出链状分子在界面附近的几何构型。后者必须借助于分子模拟。采用分子模拟可同时从微观和宏观研究固体表面上的吸附现象。MC法模拟分子在固液界面的吸附性质[37]，可测定吸附层厚度、分子构型分布、表面覆盖率等一般实验难以获得的信息。汪文川、周健用Gibbs系综的MC法模拟甲烷的吸附平衡[38]：在固定的温度下，对甲烷在1.91nm的活性炭孔中的吸附平衡进行模拟，最终得出了该温度下甲烷的吸附等温线。类似的研究还有CS2在石墨上的吸附[39]。在高压状态下，MC方法可以对气体的物理吸附进行模拟

[40]

，该研究的目的有两个：得出分子分布函数以及单层覆盖等用实验方法不易获

得的数据；在完全确定的体系下，得出相关的吸附数据。另外，在高压下，还可以对稀有气体的吸附进行MC模拟[41]。

1.1.3流体相平衡分子模拟研究进展

过去几十年中，计算机模拟由个人研究变成了有组织的合作研究，现在随着计算机技术的发展，人们已经把计算机分子模拟看作是与理论和实验这两种传统研究方式相平行的一种方法。计算机分子模拟所能研究的体系也越来越多，所研究的体系也变得更加复杂。从最初的单原子体系到双原子体系，再到碳氢化合物体系，高分子聚合物体系，多元混合物体系。流体相平衡性质研究的目的在于揭示相平衡时流体的温度、压力、各相体积、各相组成以及其它热力学函数间的相互关系及规律性。多年以来，精确描述流体的相平衡性质一直是热力学的重要课题。由于在计算机分子模拟的方法中，引入的近似条件很少，可以得到在较大范围内的可靠数据，因此用计算机分子模拟方法来描述流体相平衡的性质一直是人们所追求的目标，计算机分子模拟方法在流体相平衡研究中所起的作用也变得越来越重要，它从分子间的相互作用入手，利用严格的统计力学规律和先进的计算机编程技术，对流体的热力学性质和相平衡进行模拟计算。一开始人们一般采用

6

重庆大学硕士学位论文1绪论

计算流体化学势的方法来确定相平衡行为，给相平衡的模拟带来诸多不便，后来发展了吉布斯蒙特卡洛模拟方法、反应吉布斯蒙特卡洛模拟方法及巨正则系综蒙特卡洛模拟方法等。

随着对势能函数的深入了解，分子模拟技术越来越多地应用于实际流体的研究，如液态水等物质的蒙特卡洛模拟，得到了与实验一致的热力学性质，但目前对于实际溶液体系的研究尚存在一些困难。因此，对于化工中常见的物质，如水、碳氢化合物等实际溶液体系的分子模拟研究具有普遍的意义。目前，计算机的计算能力已大大增强，国内外出现了许多分子模拟小组，计算机分子模拟研究正在朝阳蓬勃地发展。现在Gibbs系综MonteCarlo方法在相对简单的气－液、液－液平衡方面的模拟也取得一定的成功。在分子模拟研究流体相平衡性质方向主要集中在：改进和建立各种用于流体相平衡模拟的分子间作用势能模型和分子模型；改进和建立分子模拟中的合理构象的抽样方法；采用各种分子模拟方法计算和预测流体的相平衡、优化分子力学参数及修正与改进其他相关的热力学理论方法。

1.2研究目的和内容

从前人的研究工作可以看到，计算机分子模拟已成为在分子水平上研究流体结构及其性质的一种强有力的工具，越来越受到包括化学工程在内的一些科学领域的重视。通过分子模拟不仅可以得到一些宏观性质，还可以得到实验无法得到的微观或介观结构图像，以便分析想象和机理之间的内在联系。计算机分子模拟在这些领域显示出日益重要的意义和广阔的前景。

流体相平衡性质揭示的是流体处于相平衡时，其温度、压力、体积、各相的组成及其他热力学函数之间相互关系的规律。它是蒸馏、吸收、结晶、萃取等单元操作的理论基础，因而对许多工业过程的设计、优化和操作是必不可少的。

本文依据分子模拟的基本原理，采用吉布斯系综蒙特卡洛模拟方法，运用VisualFortran95结构化程序语言，系统地模拟计算了水/甲烷、水/二氧化碳的二元体系以及正构烷烃混合物乙烷/庚烷和辛烷/十二烷的二元体系的气液相平衡数据。本论文的主要内容包括建立模型分子结构及其参数、确定势能模型及其参数、确定模拟方法、编写模拟程序及模拟结果分析与讨论。本文为以后利用分子模拟方法进一步研究预测更复杂的碳氢化合物的多元体系的相平衡的性质以及更为复杂的水的二元体系的热力学性质积累了宝贵的经验。

7

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

2.1系综原理

分子模拟的理解主要从系综开始。Gibbs在1901年提出了系综的概念：将大量宏观状态相同的系统集合在一起，而每一个系统各处在它所辗转经历的某一个微观状态中，这样的系统称为标本系统或者系综样本系统，由此而构成的标本系统的集合就叫做统计系综(简称为系综)。统计力学旨在根据物质的微观状态的动力学行为来解释宏观状态的性质，其基本原理就是系综原理[42,43]。热力学状态是由一系列参数确定的，热力学性质可由状态方程和热力学基本方程得出。虽然象扩散系数D这样的量与体系的微观结构和动力学性质有关，但它也是状态函数，其值完全由少数变量(如N、P、T)决定，而不是由瞬间的力学状态比如原子位置、动量等决定。

在相空间中某一点A，其某一瞬态性质M的值记为M(A)，体系的宏观观测值Mobs是各瞬时M(A)的平均值：

t

MobsMtimeM(A(t))time

1

limM(A(t))dttt0

(2.1)

大量分子体系的M(A(t))十分复杂，要进行时间平均是不可能的。因此采用系综来进行计算，因为标本系统的数目N很大(N的数量级通常为1024)，足以代表所有的微观状态，所以对于给定的热力学系统，某一性质的时间平均值就近似等于它的系综平均值。通常将系综看成是相空间所有点A的集合，各点按概率密度(A)分布。它由指定的宏观参数N、V、T确定，因此采用符号NVT来标记，一般用ens表示。每一点代表任一特定瞬间的系综，每一体系随时间变化，遵循力动方程而与其它系综无关。如果ens(A)表示某一平衡系综，即ens(A)/t=0，则体系的变化显得相当特殊：某一体系离开状态A()到另一状态A(+l)，就有另一体系从A(-l)状态进入来代替它。如果通过相空间的所有点中只有一个轨迹ens(A)为非零，那么每一体系最终在所有点都可预见，这称为“各态遍历”。然而我们关心的是平均性质，可以对系综的每个成员的平均来代替时间平均：

MMensM/ensM(A)ens(A)

A

N

E(pA,rAN)exp[]

KBT

ens(A)

Qens

(2.2a)

(2.2b)

8

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

Qens

E(piN,riN)

exp[]KBTi

(2.2c)

式2.2a～2.2c中，KB是Boltzmann常数；T体系绝对温度；pN、rN及E(pN，rN)分别表示体系的广义动量、广义坐标及势能；ρens(A)是体系空间构型A存在的概率；Qens是体系的配分函数。

系综等价性使我们能在各种情况下确定一个合适的相函数，体系的基本热力学量可以在任一合适的系综中计算。系综的维里就是各个粒子的坐标和作用在该粒子的力的乘积之总和的期望值。

体系的动能、势能及总能量可用相函数来计算：

EHk

(2.3)

可用维里定理计算温度和压力：

Pk

H

PkHqk

kBT

(2.4a)

qkkBT

Hqk

(2.4b)

式中：qk为普适坐标，Pk为动量。这些方程是M

kBT

M

的特例，它们qk

可由广义正则系综导出，且对任一系综都有效。对原子体系：

p

i1

N

2i

/mi2k3NkBT

(2.5)

这就是著名的能量均分原理：每一自由度平均能量为0.5kBT。原子j作用在原子i的力的公式为：

fijrijv(rij)

1fij

rij

dv(rij)(rij)rijrij

2drrijij

(2.6a)(2.6b)

压力可从维里函数计算：

1N

Wrifi

3i1

(2.7)

9

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

PkBT

WV

(2.8)

式中fi为分子i所受到的作用力。对于两两相互作用体系，维里函数可表示为与初始坐标无关的形式：

rfr

ii

i

i

ji

ij

fij

(2.9)

式中rij=ri-rj，fij为原子i对原子j作用力。

若分子模型采用点点(site-site)模型，则仍可用上式计算维里函数，只是在取加和时应遍取ia，ib对。对于分子体系：

W

Wab(rab)11

W(r)(rr)ijabij23iji3ijiabrabdU(rab)

Wab(rab)rab

drab

(2.10a)

(2.10b)

式中下标a，b指分子的不同作用点。rab=ria-rjb，为i、j分子间作用点之间的距离向量，rij是i、j分子质心之间的距离向量。

体系的压力计算式为：

PVNkBTW公式2.11中的〈

〉表示系综平均。

(2.11)

2.2系综和MC方法

分子模拟中常用的系综有微正则系综(NVE)、正则系综(NVT)、NPT系综、巨正则系综(VT)及Gibbs系综等，另外还有一些由这当中的某一类派生或拓展而来的系综，下面对几个比较常用的系综及MC方法进行介绍。

2.2.1正则系综(NVT)MC法

将大量体积为V、粒子数为N的标本体系(它们之间通过导热壁隔开，可彼此传递热量，却无粒子透过)堆积在一起，并放入一温度为T的大热浴中，使所有标本体系达到热平衡，标本体系之间有能量交换，而无物质交换。这种标本体系的集合即为正则系综[44]。

对于正则系综体系中只有粒子的位置移动，没有体积变化和物质交换。根据统计力学理论，其配分函数为：

QNVT

1

3N

drexp[u(r)]N!

(2.12)

式中u(r)为标本体系构型r时的粒子间相互作用势能。

在MC模拟方法中，结合Metropolis重要性抽样方法，其基本算法如下：在

10

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

模拟盒子中随机选择一个分子，计算此构型的能量；对此分子的位置进行随机扰动，即给此分子一随机位移，得到新的构型。粒子移动的接受概率为：

Pmovmin1,exp(U)式中，U为由于粒子位置扰动一起体系构型能的变化。

(2.13)

2.2.2等温等压系综(NPT)MC法

若标本体系之间通过柔性的导热壁格隔开，除彼此可传递热量外，还可通过体积扰动来保证标本体系内的压力不变，这种标本体系的集合即NPT系综[44]。NVP系综方法的实现与NVT系综类似，只是标本体系压力为恒定值P。它包含两种扰动：粒子的位置移动和体系的体积变化。粒子位置扰动与体积变化扰动按某一权重随机选取。根据统计力学理论其配分函数为：

1

QNPT3NdVexp(pV)drexp[u(r)]dVexp(pV)QNVT

N!

VV

Pvolmin1,expUpV(N1)ln

V

(2.14)

其中粒子位置扰动接受概率与NVT系综相同，而体积扰动的接受概率为：

(2.15)

由于实验测定通常在定压情况下进行，而且混合理论常在这一假设下得到的，所以NPT系综特别适用于模拟混合物。模拟时，通过计算第二维里系数将得到的压力与给定压力作比较，以检验计算程序的正确性。

2.2.2巨正则系综(μVT)MC法

将大量达到平衡的体积为V、温度为T和化学势为μ的标本体系从大热浴中拿出并与大子库分开，然后堆积在一起，置于一与外界绝热且不能传递物质的屏蔽层内，而标本体系之间可传递热量及物质。对于每个标本体系，其余标本体系即是它的大热浴和大子库，这样的系综称之为巨正则系综[44]。

对于巨正则系综，温度T、体积V和化学位是恒定的，但是体系的粒子数会在平衡值附近发生涨落。根据统计力学理论，该系综的配分函数为：

QVT

N

1N!

3N

exp(N)drexp[u(r)]exp(N)QNVT

N

(2.16)

巨正则系综包括三种扰动：粒子的位置扰动、粒子的增加扰动和粒子的删除扰动。

(1)粒子位置的扰动：与正则系综中的粒子位置扰动方式一样，GCMC使用一般的Metropolis方法，其接受概率为：

Pmove(mn)min1,exp[Unm/KT](2.17)

11

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

(2)粒子的增加扰动的接受概率为：

V

Pcreatemin1,3exp[(Unm)/KT]

(N1)

(3)粒子的删除扰动的接受概率为：

(2.18)

Pdelate

NV3min1,exp[(Unm)/KT]

V

(2.19)

其中Unm代表相应于各种尝试的势能变化。

从巨正则系综的定义可以看出，这类系综可应用于与环境既有能量交换又有物质交换的开放体系。为满足微观可逆性，必须使得体系内增加粒子的扰动尝试概率与删除粒子的扰动尝试概率相等。模拟中广泛采用等权重选择的方法随机选择其中一种进行尝试。

GCMC模拟的优点是通过模拟能直接得到系统的统计自由能及其相关性质，比较适合于对非均相系统进行模拟研究。缺点是当密度较高时，此方法失效，因为此时成功插入或删除粒子的概率很小。因此GCMC系综还借助特殊的抽样方法，如偏倚抽样、伞形抽样或者优先抽样等。

2.3Gibbs系综MC法

Panagiotopoulos于1987年提出了一种直接模拟计算流体汽液平衡的新方法，称为Gibbs系综MC法

[45]

，即GEMC法。GEMC法是流体相平衡模拟方面的一个

重大突破，它巧妙地回避了计算化学势的难题，汽液平衡的模拟在一次计算中就能完成。由于在吉布斯系综蒙特卡洛方法中不需要提供许多高度精确的数据，相对比较容易用计算机程序进行描述，可以减少计算机的运算时间，提高计算效率，目前Gibbs系综MC法在气-液、液-液平衡方面的模拟非常成功，已经被广泛用来获得流体混合物的相平衡性质。

2.3.1GEMC法的实施过程

GEMC模拟方法的基本示意图如图2.1所示。Gibbs系综MC法在模拟流体相平衡时引入一虚拟相界面将模拟盒分成两个模拟元胞，分别表示体系中的两相。此模拟盒相当于宏观相平衡体系中汽液相界处的一个微元，其特征是微观足够大，满足统计的要求；宏观相对较小，满足计算机计算能力的要求。

在模拟过程中，对给定体系的温度T，两相满足压力和化学位都相等的条件可通过三种扰动来实现：其一，粒子位置扰动；即随机选择的粒子的移动，达到两个区域各自内部的平衡。其二，粒子体积扰动；即两个区域的体积变化，达到两个区域的压力平衡。其三，粒子转移扰动；即随机选择的粒子从一个盒子移到另

12

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

图2.1Gibbs系综方法平面示意图，图中虚线表示周期边界Fig.2.1IchnographyforGibbsensemblemethodology,thedash

denoteperiodicboundaries

一个盒子，最终达到两个区域的每一组分的化学势都各自相等。

在吉布斯系综蒙特卡洛模拟方法中，这三个蒙特卡洛运动的接受规则最早是由Panagiotopoulos[45]通过涨落理论推导而得到的。后来，Smit[44]等人作了进一步的研究，给出了系综的统计力学的定义，并由此推导得到吉布斯系综的统计分布函数。

考虑一体系处于恒定温度T、总体积为V、总粒子数为N，体系被虚拟相界面分成的两个元胞的体积分别为VI及VII(VII=V-VI)，粒子数目分别为NI及NII(NII=N-NI)。

根据统计力学理论，体系的配分函数

Q

1

[45]

为：

VN!NINIIN1dVdVdVexp[U(N)]dIIIIIII3NN!N!NII!0

NII

exp[UII(NII)]dII

(2.20)

NVT-GEMC与NPT-GEMC模拟时，体系的某一构型存在的概率分布密度P(NI，VI，N，V，T)及P(NI，VI，N，P，T)分别正比于函数式(2.21a)与(2.21b)[4]：

13

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

p(NI,VI,N,V,T)

N!

exp[NIlnVINIIlnVIIUI(NI)]UII(NII)(2.21a)NI!Nii!

p(NI,VI,N,P,T)

N!

exp[NIlnVINIIlnVIIUI(NI)]UII(NII)

NI!Nii!(2.21b)

P(VIVII)

GEMC模拟三类运动的接受规则为：

Paccmin1,B式中2.22中，B为准Boltzmann因子。

(2.22)

UUold

Bexpnew

KBT

(2.23)

接受规则Pacc的意义是：在GEMC模拟中，老构型向新构型作尝试运动，如果B>1,那么就认为尝试运动成功。如果B<1，则产生一个[0，1]上的随机数ξ与B比较，如果B>ξ，则认为尝试运动成功。如果尝试运动成功，则新构型就作为下一次构型尝试运动的起点，否则，尝试运动前的老构型将作为下一次构型尝试运动的起点。2.3.1.1NVT-GEMC模拟

对于NVT-GEMC模拟，式2.23中的参数的表达如下：

(a)粒子位置扰动，=1；

VIV

(b)粒子体积扰动，V

I





N

VIIVVII



NII

；

(c)粒子转移扰动，

NIIVI

；

(NI1)VII

(1)粒子位置扰动的接受概率

在这一步，两模拟盒以完全相同的机理变化，现仅以模拟盒I为例进行说明。此模拟盒是一个NIVIT为常数的正则系综系统。在这个系综里，构型能为UI的状态发生的概率正比于exp(-UI)。随机选择一个分子，尝试对其作随机空间位

move

置变化移动，移动的接受概率pacc的计算公式为：

move

paccmin1,

exp(-Unew)

min1,exp(U)exp(-Uold)

(2.24)

(2)粒子体积扰动的接受概率

14

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

NVT-GEMC模拟法主要用于纯流体相平衡模拟。其两个模拟盒的体积变化协同进行。模拟盒I一次随机的体积变化量为V，由于模拟元胞盒子总体积保持不变,因此模拟盒II同时进行-V的体积变化,与模拟盒I的处理方法相同,可得模

NVT

拟盒II新旧状态发生的比率公式。模拟盒I、II新旧状态的接受概率pvexchange的计

算公式如下：

VIV

(UU)Nln()IIIIVI

min1,exp

VIIV

Nln()VII

NVT

pvexchange

(2.25)

(3)粒子转移扰动的接受概率

在两子体系之间进行分子转移或交换，使各组分在两相中的化学势相等。其基本过程是：首先，随机选择一个模拟子区域；然后，随机选择此子区域中的分子；最后，对结构简单的小分子，通过在一个模拟盒中随机插入一个某类型的分子与另一个模拟盒中随机删除一个相同类型的分子来实现分子转移。对结构较为复杂的分子，尝试随机分子插入造成大量的分子重叠，使分子转移失败，可以通过分子交换的方法进行改进。分子交换是通过把一个子体系中的小分子换成一个较大的分子，其交换的方法是：在其中一相中拆除一小分子并在拆除小分子的空间处增长出一个较大的分子来，而在另一相中拆除一个较大分子并在拆除较大分子的空间处增长出一个小分子，从总体上实现了分子间的交换。分子交换示意图

NVT

如图2.2。NVT-GEMC中粒子转移扰动的接受概率ptransfer的计算公式为：

图2.2分子交换示意图[44]

Fig2.2Sketchmapformoleculesexchange

NVT

Ptransfermin1,

NIIVI

exp(UIUII)

(NI1)VII

(2.26)

体系的压力P可用体系体积波动ΔV的函数计算：

15

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

uKTVV

P(V)BlnexpVVKBT

起的能量变化量；KB是Boltzmann常数；T是体系绝对温度；符号〈综平均。

2.3.1.2NPT-GEMC模拟

对于NPT-GEMC模拟，参数的表达如下：

(a)粒子位置扰动，=1；(b)粒子体积扰动，

VIVI

V

I





N

N

(2.27)

式中V、N和Δu分别模拟体系的体积、粒子数目和由于模拟体系体积波动ΔV引

〉表示系

VIIVIIVII



NII

exp[P(VIVII)]；

ANIINIB

(c)粒子转移扰动，；AB

(NI1)(NII1)

NPT-GEMC模拟中粒子位置扰动的接受概率和NVT-GEMC模拟中粒子位置扰动的接受概率的计算公式相同。

对于粒子体积的扰动，NPT-GEMC在模拟模拟混合物时，恒压法可以指定共存相的压力P，比恒容吉布斯系综蒙特卡罗模拟法方便。此时，总体积变化不再是常数，两个模拟盒的体积变化进行。在计算新旧状态发生的比率时，可以将两个模拟盒分别计算，也可以将两个模拟盒视作一个整体计算。模拟盒I、II新旧

NPT

状态的接受规则pvexchange的计算公式如下：

VIVIVIIVII

)Nln()1,expNPT(UIUII)Nln(V(2.28)pvexchangeminVIII

P(VV)III

NPT

NPT-GEMC中粒子转移扰动的接受概率pexchange的计算公式为：

p

NPT

exchange

A

NIINIB

min1,exp(UU)IIIAB

(NI1)(NII1)

(2.29)

式2.20-2.29中,是deBroglie波长；ξI及ξII分别表示粒子在两子区域I和II中的广意空间坐标；NI、NII、VI、VII分别为两子区域的分子数目及体积；

=1/KBT，KB是Boltzmann常数；T体系绝对温度；ΔU是由于粒子平动与转动引起的构型能变化；ΔUI、ΔUII、ΔVI、ΔVII、NIA及NIIB分别表示两子体系的位能变化、体积变化及组分A、B的粒子数目。

16

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

2.3.2抽样方法

蒙特卡洛方法是用随机变量的样本来推断其分布性质的。由于样本的总量有限，因此这种由部分推断整理的结果势必造成误差。己经证实，用蒙特卡洛方法估算的流体某一性质的标准误差与微观构型数的平方根成反比。因此要提高模拟结果的精度，抽样方法就显得非常重要。

GEMC法采用的抽样方法主要有Metropolis重要性构象抽样方法和构型偏倚MC构象抽样方法。原始的吉布斯系综MC模拟中采用的是Metropolis重要性构象抽样方法，由于其对分子交换步中合理的分子构象的搜索效率比较低，只能应用于模拟极为简单的分子体系的汽液相平衡。由于其对分子构象的搜索不在是完全随机的，构型偏倚构象抽样方法，有利于分子避免相空间的障碍，对分子的合理构象的搜索效率较高。

蒙特卡罗模拟通过体系分子空间构型的随机运动而产生大量的微观构型样本。体系的分子从某一空间构型i向空间构型i+1的运动转移通常应满足细致平衡条件。细致平衡条件可以表达为：体系从相空间的某一相点o向另一相点n转移运动的流K(on)，应等于体系从相点n向相点o转移运动的流K(no)。

K(on)K(no)

(2.29)(2.30)(2.31)

K(on)N(o)(on)Pacc(on)K(no)N(n)(no)Pacc(no)

式2.29-2.31中，N(o)和N(n)分别表示相点o和相点n存在的概率；(on)和

(no)分别表示体系从相点o向相点n转移及从相点o向相点n转移的概率；Pacc(on)和Pacc(no)分别表示从相点o转移到相点n及从相点o转移到相点n的接受概率。对于GEMC中Metropolis重要性抽样方法，(on)=(no)。而对于GEMC中的构型偏倚(Configurational-bias)抽样方法，(on)≠(no)。2.3.2.1Metropolis重要性抽样方法

Metropolis重要性抽样方法是一种有偏抽样方法。它将抽样重点放在玻尔兹曼因子贡献较大的空间区域，以exp(-βU)为概率选择构型，U为体系内能。生成的构型仅依赖于前一个构型，即由后一个构型在前一个构型的基础上作一定位移的扰动而得。由各个构型作为链节组成了构型空间中的Markov链。应用这种方法粒子可由其原始位置逐渐移动到模拟盒子中的任何位置，得到符合玻尔兹曼分布的体系构型。

Metropolis抽样的基本步骤：假定一个体积为V的矩形箱中有N个粒子，势能可表示为Upo(rN，wN)，则新的位形按以下步骤产生：

17

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

1从N个分子中随机地抽取一个分子(也可轮流)；○

2让这个分子试着随机移动一步，即将分子质心坐标和方向改变为：○

xn=x+δ1ξ1θ+δn=θ2ξ4

yn=y+δ1ξ2φn=φ+δ2ξ5

zn=z+δ1ξ3ψn=ψ+δ2ξ6

式中：x、y、z—空间位置

θ、φ、ψ—Eular角，代表分子在空间的取向δ1、δ2—最大允许平动度和旋转度ξ1…ξ6—[0，1]上的分布的随机数

3在新的位型下，计算新的位能Upn；○

4如果Upnexp(-βUpo)，于是接受新的位型作为马尔科夫

链的一个状态点；

5如果Upn>Upo，则以正比于exp[-β○(Upn-Upo)]的几率接受新位形。为此将

exp[-β(Upn-Upo)]与一个[0，1]上的随机数进行比较，如果前者大于后者，接受新的位形；反之，则将原位形再次作为马尔科夫链的新的状态点。

Metropolis抽样方法就是从某一个初始构型(比如将N个粒子按面心立方晶格放置)出发，不断重复上述步骤产生一系列构型，最后对抽取出来的构型样本的势力学性质进行系综平均。为了获得满意的模拟结果，一般需要数百万乃至数千万个构型，且参数δ试移”被接受的次数1和δ2还要选得合理。参数的选择和接受率(“进行加和，再除以总移动次数，这一比值称为接受率)相关联。2.3.2.2构型偏倚抽样方法

在Gibbs系综蒙特卡罗模拟与巨正则系综蒙特卡洛模拟(GCMC)两种模拟方法中，都要进行分子交换。在Gibbs系综蒙特卡洛模拟中，分子交换既要在模拟盒中进行也要在模拟盒之间进行。这种分子交换与任何的实际运动不相关(非自然运动)，因此不能用MD方法模拟其相平衡，MC模拟却能实现其非自然运动。然而对于分子几何结构较为复杂体系，普通的GEMC与GCMC模拟是难以成功的，原因是对于大分子来说，在模拟盒子里进行完全随机尝试插入分子的接受率很低，巨大数目的尝试插入计算量是目前计算机能力所不能满足的。构型偏倚(Configurational-bias)抽样方法是普通MC尝试移动的扩展算法，其特征是插入分子的运动不是完全随机的：分子的插入运动在一定取向的约束下进行，以至于要插入的分子有一个更大的概率来适应已存在的构型。采用Configurational-bias抽样方法的代价是使计算程序更加复杂。但是，在Configurational-bias抽样方法的帮助下，常常有利于加速计算，模拟计算分子结构更为复杂的分子体系。

对于强烈倚靠于分子相对取向的相互作用势能的分子(极性分子、氢键分子等)体系，对其进行MC模拟时，要寻找一个既不和其他分子重叠，又有可能接受的

18

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

分子取向的位置是重要的。如果完全随机巧合方法来发现一个合适的取向的概率很小，则可以使用构型偏倚尝试移动来增加构型抽样的接受概率。构型偏倚取向的基本算法如下：

以o表示分子的旧构型，n表示分子的新构型。把分子的移动分解为平动及分子的取向运动两部分。对于平动部分，采用标准的随机Metropolis移动。对于移动的取向运动部分则按照如下偏倚取向方法产生尝试取向：

1在一小的随机距离内移动分子的质心，计算粒子在新构型的能量○

upos(n)。

2产生k个尝试取向{b1，○b2，…，bk},对于每一个尝试取向，计算其能量uor(bi)。3定义Rosenbluth因子为：○

W(n)exp[uor(bj)]

j1

k

(2.32)

在k个取向以如下的概率选择一个取向,比如第n个取向：

exp[uor(bn)]

P(bn)

W(n)

(2.33)

4对于旧构型o,粒子在旧位置的取向记作b0,能量用upos(o)表示。尝试生成○

的k-1个虚假取向，记作b2，b3，…，bk，使用这k个取向可得到旧构型的Rosenbluth因子W(o)：

k

W(o)exp[u(b0)]exp[uor(bj)]

or

j2

(2.34)

5分子从旧构型o到新构型n的移动的接受概率为：○

W(n)pospos

Paccmin1,exp{[u(n)u(o)]}W(o)

(2.35)

2.3.3分子间力和势能函数

流体混合物的热力学性质，依赖于分子的本性和分子间的相互作用或者分子间力。由系综理论可知，只要知道由分子间力决定的系统的位能，原则上就可通过配分函数求得各位形的热力学性质。从力的本性来分类，目前已经被认识的有万有引力、强核力、弱核力与电磁力。万有引力太小，对能量的贡献几乎可以忽略不计。强核力和弱核力的作用距离都在10-4纳米左右，分子间的距离通常都比这要大得多；只要电磁力能够正确解释分子间的相互作用。

从现有的知识来看，分子间力主要有：

(1)长程的引力，其中包括静电力、诱导力和色散力。色散力是指由于分子中电子和原子核不停地运动，非极性分子的电子云的分布呈现有涨有落的状态，从

19

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

而使它与原子核之间出现瞬时相对位移，产生了瞬时偶极，分子也因而发生变形。分子中电子数愈多、原子数愈多、原子半径愈大，分子愈易变形。瞬时偶极可使其相邻的另一非极性分子产生瞬时诱导偶极，且两个瞬时偶极总采取异极相邻状态，随时产生的分子瞬时偶极间的作用力。诱导力是指极性分子与非极性分子相遇时，极性分子的固有偶极产生的电场作用力使非极性分子电子云变形，且诱导形成偶极子，固有偶极子与诱导偶极子进一步相互作用，使体系稳定。

(2)短程的斥力；当两个分子非常靠近以至电子云相互重叠就会产生这种很强的排斥力。

另外，一种较弱的化学作用如氢键和电荷转移，由于它们的键能和通常的分子间相互作用差不多大小，因此也常作为分子间力进行研究。

体系可分为原子体系，分子体系和晶格体系。原子体系，势能可分成单原子、双原子、三原子作用等项：

uu1(ri)u2(ri,rj)u3(ri,rj,rk)＋……

(2.36)

式中第一项表示外部场的影响，第二项是最重要的，它仅取决于rij=|ri-rj|,可记为u2(rij)。对液体，方程2.12的u3(ri,rj,rk)项无疑是重要的，约占总势能的10%。四体贡献比三体贡献小得多。因为对三分子求和很耗计算机时间，三体贡献项很少被包含在计算机模拟中。二体近似已能对液体作很好的描述，因为三体影响可被部分包括于定义为“有效”分子对势能中。于是，可将方程2.36重新写为：

uu1(rj)u2(rij)

eff

(2.37)

2.3.3.1势能函数

分子模拟的基础是准确可靠的势能函数，一对分子间相互作用的势能随分子间距的变化称为势能函数。分子的势能Up与分子间距r的函数关系可以构成势能能函数，在建立一个完整的位能函数时，分子的势能通常包括吸引力贡献和排斥力贡献。分子模拟首先要建立模型分子。分子的几何结构、分子在空间排列构型及分子体系的势能都受到分子力场的约束。分子建模就是在分子力场约束的条件下，产生分子简化了的、较为合理的分子几何结构。模拟的目的不同，所建立的分子模型简化程度也不一样。流体相平衡的模拟结果是否正确，其中一个关键问题是看所建立的分子力场是否能正确描述分子间及分子内部的各种相互作用。

(1)分子间相互作用

常用的分子间相互作用势能函数有L-J势能函数模型与Bexp-6势能函数模型：

1Lennard—Jones势能函数○

Lennard-Jones(简写为L-J)势能函数是目前最为广泛的，用于计算分子间相互

20

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

作用势能的函数。对于单原子体系，L-J势能的函数的表达式为：

126U4

rr

(2.38a)

其函数图形如图2.3。在图2.3中，r为两分子间的距离，repp表示原子间势能的排斥贡献，attrp表示原子间势能的吸引贡献。当分子距离缩小时，分子间的吸引力和排斥力都同时增大，但排斥力增加较快。当分子距离非常接近时，分子间的作用力主要为排斥力。当排斥项与吸引项相当时，U=0，r/=1。分子的碰撞直径

(也称为不可压缩直径)就是当体系的势能为零时的分子直径。势阱深度表示处于势能曲线最低点的分子与处于势能零点的分子之间的能量差。在模拟中，势阱深度通常与为玻尔兹曼常数kB组合在一起，即／kB，其单位为温度单位Kelvin。

与／kB是模拟中最基本的两个作用参数。

对于多原子分子体系，分子间的相互作用可以看成是分子的各种作用单元之间相互作用的总和，再把分子的各种作用单元简化成准原子，这样，多原子分子间的相互作用势能就可以近似地用L-J势能模型计算。其分子间的相互作用势能函数式可以表示为：

Uija

bn

m

ab

4ab

rab

abrab

12





6

abqq

40rab

(2.38b)

12

rep4()p

r

Uatr64()r图2.3Lennard—Jones势能函数图Fig2.3Lennard—Jonespotentialfunction

式中，Uij是分子i与j之间的相互作用势能;a,b分别代表分子i与j的作用基团;m,n分别是分子i、j之间相互作用点数目；qa、qb分别是分子i、j之间相互作用基团所带的电荷；0是真空介电常数；ab、ab分别是作用势的尺寸参数及能量参数。作用单元可以是原子核、原子团(CH3-、NO2-等)、孤对电子或者是化学

21

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

基团，视具体物系模型分子建立而定。比如，在处理丙烷、正丁烷等烷烃分子模型时，就把CH3-看成一个作用单元，而在处理二氧化碳分子模型时，把碳原子和氧原子分别看成不同的作用单元。

多原子分子内部及其分子间的不同作用单元之间的相互作用势能参数可采用Lorentz-Berthelot混合规则计算：

ij(iijj)2(2.39)(2.40)

ij(iijj)

2○

BuckinghamExp-6势能函数

1

2

BuckinghamExp-6(简写为Bexp-6)势能函数比L-J势能函数多一个用于表征分子的作用单元的柔性的参数，把分子间的作用单元看作柔性作用单元。因此，在排斥作用方面的描述比L-J模型具有更好的理论基础。在考虑吸引贡献时和L-J模型一样，仅考虑了六次项的影响。Bexp-6函数比L-J函数复杂，采用Bexp-6函数进行分子模拟计算时需要用到复杂的数值分析方法，难以编程，计算相当耗时，以前很少以Bexp-6势能函数计算流体的相平衡。随着计算机的发展，近年来相关方面的介绍越来越多。如，Kuijperetal[47]用Bexp-6势能模型检验了混合物之间的纯排斥作用。Dodd和Sandler[47]利用Bexp-6势能模型，采用MonteCarlo方法计算了混合物的构型能、压缩率和径向分布函数。Traveres和Sandler及Errington和Panagiotopoulos[48]使用Bexp-6势能模型检验了纯流体的气液相平衡。

Bexp-6势能函数的表达式为：

6

6rmqiqjr,exp1ru(r)16/r4rm0

forfor

rrmaxrrmax

(2.41)

式2.41中，u(r)表示分子间相互作用势能；α是表征势能排斥部分的柔性参数,

0.5

αij=(αiiαjj)；截断距离rmax是方程du(r)/dr=0的最小正根，可以通过迭代方程

du(r)/dr=0计算；rm是分子间作用最小时分子间的距离；qi、qj分别是分子i、j之间相互作用基团所带的电荷。当分子间距离右逼近rmax时，分子间的势能增长很快。当分子间的距离小于rmax时，Bexp-6函数变为减函数，不满足分子间相互作用规律。(2)分子内部作用

分子内部存在伸缩势、键弯曲势、扭曲势。当分子内部任意两个作用单元之间相隔三个或三个以上的作用单元时，可以计算其间的电荷作用及非电荷作用。分子内部作用表示为：

Uint

spring

Kl(rr0)K(0)UTorUintVDWUintZZ

bend

2

(2.42)

22

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

式中，r是键长，r0是平衡键长；是键角,0是平衡键角；UintVDW是分子内vanderWaals作用势；UintZZ是分子内部哥伦布作用；UTor是分子的扭曲势。分子扭曲势能的计算表达式有DePablo扭曲势能和OPLS扭曲势能。

DePablo扭曲势能：

UTor()c0

cc1c

[1cos()]2[1cos(2)]3[1cos(3)]222

5

(2.43)

OPLS扭曲势能：

Utor()Ckcosk()

k0

(2.44)

式中，C0、C1、C2、C3及Ck均为扭曲势能常数，φ是分子的扭曲角。2.3.3.2周期边界条件

体系的宏观性质是体系中大量分子的微观相互作用的统计平均。由于计算机硬软件的，目前的计算机还只能处理几百至几千个粒子体系的结构及热力学性质，这个数目与热力学的宏观极限相去甚远。对于这类体系不能假定可靠的边界条件，为了消除模拟分子数与实际体系分子数相差悬殊对计算产生的影响，可让模拟盒子在三维空间无限重复，即可实现从有限向无限的过渡Metropolis等人于1953年引进了一种“周期性边界条件(periodicboundaryconditions)\"。所谓周期性边界条件就是：把小体系样品看作一个中心元胞，此中心元胞被与中心元胞相同的其它元胞所包围，构成一个无限大的三维点格。当某分子离开中心元胞时，该分子将在整个点格上移动以致它从中心元胞对面重新进入这个元胞。在对分子相互作用进行计算时，两分子至少有一个分子被选在中心元胞之中。为了便于计算，所有的计算和运动均在中心元胞之中进行，并同时考虑与中心元胞相邻的影像元胞，中心元胞中的分子不仅与其中其它分子发生相互作用，也和它自身及其它分子的最小影像发生相互作用。由此可见，其效果相当于模拟盒子无边界，从而消除了边界效应。如图2.4所示，“中心盒子”为N个分子组成的标本盒子，其周围盒子是其虚拟的映像，称为“像盒子”。“相盒子”的特征是其中的分子排布构型与“中心盒子”的分子排布构型完全相同，“相盒子”中的分子与“中心盒子”中的分子存在相互作用，相盒子无限延伸。2.3.3.3相互作用的截断

在周期边界条件下，不同盒子之间存在相互作用，因此，每一个分子与其它分子的相互作用延续至无穷远，势能的计算也变得复杂起来。但是，对于文中所应用的L-J势能函数及Bexp-6势能函数而言，势能随分子间距离的增大衰减得比较快，因此，在计算势能时可以近似地截断到适当的距离。分子之间远程相互作用可采用其它方法进行补偿计算。下面就简要介绍一下常用的两种近似方法，即

23

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

“最近映像”和“势能球形截断”。

图2.4周期性边界条件的平面示意图

Fig2.4Ichnographyofperiodicboundaries

如图2.4所示：每一个分子必有而且只有N-1个近邻分子不重复地落入以该分子为中心的盒子(不是原来的盒子，但尺寸相同)中。这N-1个分子可能是中心箱中的，也可能是像盒子中的，它们是中心分子的紧邻分子，称为最近映像(分子)。在最近映像近似中，中心分子与其余分子的相互作用能只能考虑到紧邻分子，也就是说将分子间力的作用范围截断到模拟盒表面。

在采用最近映像规则后，两两相互作用势仍然包括0.5N(N-1)项。对于几百个粒子构成的体系，其计算量仍然是很大的。为了进一步减少计算量，采用球形势能截断。将截断距离rc确定为一个常数(在程序中取模拟盒边长的一半)，作用范围则是以中心分子为圆心、rc为半径的圆球，称为截断球。在该条件下，与中心分子发生相互作用的对象仅限于球内分子。此时，与中心分子发生作用的分子数就明显小于N-1个。使用截断方法会使计算的热力学性质产生一定的偏差。为此，可以引入长程校正方法来解决这一问题。2.3.3.4长程作用的校正

为了简便计算，计算机分子模拟通常在截断距离rc内计算分子间近程相互作用。为了提高模拟的准确度，对势能的计算结果进行长程校正是必要的。MPAllen对长程作用作了论述，其计算公式[45]如下：

EfullECELRCEC2Nr2(r)dr

rC

(2.45)



(PV)full(PV)C(PV)LRC

(PV)C2Nr2(r)dr

3rC

rC

(2.46)(2.47)

fullCLRCC4r2(r)dr

24

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

式中Efull，(PV)full，full，分别表示模拟所要得到的体系的总体值，Ec，(PV)c，c是体系在rc处截断计算得到的值；ELRC，(PV)LRC，LRC即为长程校正项；v(r)与w(r)分别是分子间相互作用势能函数和维里。对于吉布斯系综模拟，由于模拟盒的体积和粒子数可以改变，它们对于瞬间的能量值和压力在是有很大的影响，从而使得模拟体系的密度波动很大，所以进行长程校正是必不可少的。由于分子间非电荷长程相互作用与电荷长程相互作用的计算方法不同，下面分别加简述。

(a)分子间非电荷长程作用校正

对于L-J球形分子体系，对L-J势能函数在rc到无穷大内积分而得到其长程校正能量。其势能的长程校正公式[45]为：

Ulrc(8/9)Nrc9(8/3)Nrc6Plrc(32/9)N2rc9(16/3)N2rc6

(2.48)(2.49)(2.50)

lrc(16/9)Nrc9(16/3)Nrc6

对于L-J二元混合体系，势能的长程校正公式为：

U/N2xixjuij(r)gij(r)r2dr

i

j

rc

(2.51)

duij(r)2

P/xixjgij(r)r3dr

rc

3drij



(2.52)

i4xjuij(r)gij(r)r2dr

j

rc

(2.53)

式2.48-2.53中，uij(r)表示L-J作用点i与j之间的相互作用势能；gij(r)是L-J作用点i与j之间的径向分布函数；表示分子的数密度；=1/KBT，KB是Boltzmann常数；N是分子数目；U、P及i分别表示能量、压力及化学势的长程较正。

无论L-J流体还是Bexp-6流体，TheodorouandSuter[47]开发了计算分子间非电荷长程作用的计算程序。本文的非电荷长程作用采用计算均采用此方法进行计算，模拟时调用此子程序就可以计算体系的非电荷长程作用。

(b)分子间电荷间长程作用校正

通常，当分子间作用单元之间的距离在大于截断距离rc时，体系作用点的径向分布函数g(r)接近于1，分子的作用单元之间的非电荷相互作用已经变得很小，可用上面简述的长程校正方法计算。而电荷间的距离大于用rc时，其间的相互作对体系总势能的贡献仍然较大。分子间电荷的长程作用可以采用Ewald总和方法

25

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

进行计算。Ewald总和的目的是计算模拟元胞中电荷之间及模拟元胞电荷与其周期映象电荷之间总的相互作用。通过建立电荷及其周期映象所产生的电场E(r)的数学模型，计算模拟体系中点电荷所在的位置的电场，最终计算出电荷间作用对分子间总的相互作用势能的贡献。

电荷间对势能的贡献可由哥伦布作用uCol计算：

uCol

1N

Zi(ri)2i1

(2.54)

(ri)为位于粒子i处的静电势

(ri)

j

zjrij

(2.55)

rij是电荷j与i的距离。模拟盒中的点电荷及周期映象的电荷密度ρ(r)可表示为：

(r)

l,imagejinl

vectors

qj(rrj)

j

qjr(rjl)

n

(2.56)

为了计算方便，可认为点电荷密度ρ(r)符合高斯(Gaussian)分布。

2

(r)qj(/)3/2expr(rjl)

lj

(2.57)

点电荷间总的相互作用Utot是电荷之间倒空间用作能部分Uq、自相互作用能部分Uself及真实空间作用能部分ΔU的总和。

UtotUqUselfU

(2.58)

式(2.59)中

Uq



1212



k

4Vk24Vk2

e

k2/4V2

i,j

qiqj

e

ik(rirj)

(2.59)



k

ek

2

/42

ˆ(k)26

重庆大学硕士学位论文2吉布斯系综蒙特卡洛(GEMC)模拟的基本理论

Uself



12qj(0)

j

12

j

q2j

(2.60)

U



1212qijrijllij

r

lij

qiqj

ij

l

erfc

rijl

(2.61)

式2.54-2.61中，加粗的参数表示矢量参数；qi、qj分别是分子i、j之间相互作用单元所带的电荷量；l表示周期边界条件；是Ewald总和参数；erfc表示错误修正函数；(ri)为位于粒子i处的静电势；k是倒倒空间矢量数目；V表示积分空间。

27

重庆大学硕士学位论文3模拟过程及程序

3模拟过程及程序

3.1模拟体系的选择

采用分子模拟技术在取得化工产品和过程开发与设计所需要的热力学性质及体系的微观组成性质等领域正起着越来越关键的作用。但是，由于分子模拟技术正处于基础研究阶段并受到多方面的约束，以及目前计算机的计算能力还远不能满足所要求的计算能力，模拟的结果不是特别精确而且效率比较低，分子模拟目前还只能够模拟较为简单或者简化的体系。水和碳氢化合物的相平衡数据在化学工业中有着十分重要的地位，因此对于水以及碳氢化合物的分子模拟研究也就相当的重要。因此，本文选择了碳氢化合物中体系相对简单且较为常见的正构烷烃的二元体系及其水的二元混合物体系进行模拟，计算了乙烷/庚烷、辛烷/十二烷，水/甲烷和水/二氧化碳的二元体系汽液相平衡性质，为以后进一步研究预测更长链的正构烷烃及更加复杂的水的二元混合物的性质打下良好的基础。

3.2流体汽液相平衡的GEMC法模拟过程

Gibbs系综MontoCarlo模拟流体相平衡的基本过程如图3.1所示。首先确定分子的势能模型和参数并建立起初始构型，然后进行GEMC法的三类扰动以使模拟体系达到相平衡的热力学条件，即对随机选择的子区域中的粒子进行随机平动与转动，使两个子区域各自内部的达到力学平衡；使两个子区域的体积随机波动，使两个子区域间达到压力平衡；在随机选择的子区域中随机地选择粒子，并使其在两子区域间进行转移或交换，使各组分在两相中的化学势渐趋于相等，最后抽取合理的样本，并进行统计计算得到体系的宏观热力学性质。

对于二元体系汽液相平衡的模拟，本文采用恒压的NPTGibbsEnsembleMonteCarlo(NPT-GEMC)结合构型偏倚蒙特卡洛法进行模拟计算。模拟体系由400个分子组成，模拟盒I、II的初始分子数目分别为300个与100个，采用随机的方法产生其初始构型。MonteCarlo运动包括44%的平动、44%的转动、2%的体积波动和10%的粒子交换。尝试选择每种组分进行平动和转动的概率相同。尝试选择两种分子进行交换的概率为1:1。使用构型偏倚抽样法提高粒子交换的成功率，较快地实现体系化学势达到平衡，提高计算时间效率。两百万MonteCarlo步用于达到相平衡，四百万MonteCarlo步用于计算平均相共存性质。采用L-J势能模型或Exp-6势能模型计算其分子间的近程相互作用势能。分子间的非电荷长程作用均采用TheodorouandSuter[47]论述的方法计算。电荷间的长程作用采用Ewald总和方法计算，计算电子间长程作用所用的周期边界条件由276个K矢量构成。模拟盒I、

28

重庆大学硕士学位论文3模拟过程及程序

分子势能模型及参数的确定分子建模与初始构型安排一定机率的GEMC三类运动

两相随机平动与转动，达到内部平衡两相体积波动，达到压力平衡两相间分子交换，使各组分在两相中的化学势相等抽取合理的样本计算合理样本的宏观热力学性质统计计算体系的宏观热力学性质模拟结果分析图3.1吉布斯系综蒙特卡罗模拟基本过程示意图

Fig3.1TheprocessofGEMCSimulation

II的初始体积由计算时所给的初始分子密度计算而得。模拟体系的模型分子结构采用构型偏倚蒙特卡洛增长方法产生。分子间相互作用单元的势能模型参数采用Lorentz-Berthelot混合规则进行计算。

3.3分子模型及其参数

(1)二氧化碳

二氧化碳分子采用三作用单元的EMP2模型[48]，如图3.2所示。

29

重庆大学硕士学位论文3模拟过程及程序

OCO图3.2二氧化碳分子的模型结构示意图Fig3.2SketchofCarbondioxidemolecule

模型参数如表3.1所示。

表3.1CO2分子的势能模型参数Table3.1PotentialparametersofCO2molecule

EMP2模型

rc-o1.149

σc(Å)2.757

εc/kB(K)28.129

σo(Å)3.033

εo/kB(K)80.507

qc(e)0.6512

qo(e)-0.3256

(2)水

水分子采用固定点电荷模型[49]，如图3.3所示。

HO图3.3水分子的模型结构示意图Fig3.3Sketchofwatermolecule

H模型参数如表3.2所示。

表3.2H2O分子的势能模型参数Table3.2PotentialparametersofH2OMolecule

EXP-6模型

ROH(Å)1.0668

σ(Å)3.1947

ε/kB(K)159.8

ΦH-O-H109.47°

q(e)0.3687

12

(3)正构烷烃

正构烷烃分子视为非线性、非极性质分子；其基团CH3-与CH2-分子别被简化为分子内部的作用单元，如图3.4所示。

CH2CH3CH2CH3图3.4正构烷烃分子的模型结构示意图Fig3.4Sketchofn-Alkanesmolecule

30

重庆大学硕士学位论文3模拟过程及程序

模型参数如表3.3所示。

表3.3正构烷烃分子的势能模型参数Table3.3Potentialparametersofn-Alkanesmolecule

模型L-JEXP-6

甲基

亚甲基

ε/kB(K)

98129.6

σ(Å)

3.753.679

--16

ε/kB(K)

4673.5

σ(Å)

3.954.00

--22

在表3.1-表3.3中，ε/kB、σ和q分别代表能量参数、几何参数和电荷量。

3.4模拟的主要程序及框架

本文的模拟程序是对A.Z.Panagiotopoulos等人编写的GEMC模拟基本程序进行修改而成的。整个模拟程序中包含一个主程序及多个子程序，各个子程序分别执行不同的功能，因此程序的具有较好的可读性和移植性。整个程序主要包括以下一些程序单元：Gemc主程序、Ranf子程序、LRC子程序、Displace子程序、Transfer子程序、Volume_NPT子程序、Volume_NVT子程序和Creat子程序等。Gemc主程序用来控制对参数的输入及调用各类子程序，计算所研究体系的热力学性质；Ranf子程序是用来产生随机数；LRC子程序用来进行长程校正；Displace子程序用来进行所模拟体系的粒子的平动；Transfer用来进行所模拟体系的粒子在两个模拟盒之间的交换；Volume_NPT子程序和Volume_NVT子程序分别进行等压和等容条件下模拟体积变化；Creat子程序用来生成所模拟体系的初始构型。

3.4.1Gemc主程序

主程序是整个程序的核心部分，对子程序的调用进行控制。首先，它对各个参数进行初始化，调用输入文件以此输入势能模型参数、分子结构参数及分子力学参数，产生分子结构模型，创建模拟体系的初始构型，并计算初始构型的总能量。接着，调用三类移动步的子程序进行三类吉布斯系综的扰动。最后，计算新构型的能量，维里等相关参数，通过系综平均最终计算出体系在平衡后的各项热力学性质，从而得到我们所研究的流体的性质，达到模拟的目的。下面的程序是主程序的简化，它代表了主程序的主体计算部分。PROGRAMCallCallCall

Gibbs_main

#Gibbs主程序#创建体系的初始构型#计算初始构型能量#MC循环

input(E_patrameters,M_parameters,F_patameters)#输入模型参数Creat(Grow(molecules),system(rx,ry,rz))interaction(initial_Energy)do

icycle=1,ncycle

31

重庆大学硕士学位论文3模拟过程及程序

ran=ranf()*(nmove+nvol+nswap)if(ran.le.nmove)thencallcallelsecallendifcall

enddoend

programGibbs_main

sample()

#抽样系综平均

transfer()

#调用分子交换子程序

displace()volume()

#调用分子移动子程序#调用体积变化子程序

elseif(ran.le.(nmove+nvol))

3.4.2GEMC法的三类扰动子程序

(1)Displace子程序

Displace子程序是对随机选择的子区域中的粒子进行随机平动与转动，使两个子区域各自内部的达到力学平衡，程序简化如下：Subroutinedisplace

calldo

toterg(returni=1,nrxyzold=rxyz(i)call

single_energy(returnEold

,Wold)

#移动分子I

rxyznew=rxyzold+(2.0*ranf()-1.0)*drmaxrxyznew=rxyznew-anint(rxyznew/boxl)*boxlcall

single_energy(returnEnew,

Wnew)

deltv=Enew-Eolddeltw=Wnew-Wolddeltv=deltv/temp

if(ranf().lt.exp(-deltv).and.deltv.lt.125.0)then

En=En+deltvWn=Wn+deltwrxyz(i)=rxyznewendifend

do

#分子势能平均

En=En/real(numb1+numb2)

press=real(numb1+numb2)*temp/vol+Wn/vol.or.

32

En,Wn)

#接受的移动的条件

重庆大学硕士学位论文3模拟过程及程序

press=f(deltv)end

subroutinedisplace

#计算体系压力

(2)Volume子程序

Volume子程序的功能是对随机选择的子区域中的粒子进行随机平动与转动，使两个子区域各自内部的达到力学平衡，程序简化如下：SubroutineVolumeCall

toterg(box1,box2,en1o,en2o)VI，II(new)=VI，II(old)+(Raf()-0.5)*VmaxFact=VI，II(new)/VI，II(old)rxyzI,II(new)=rxyzI,II(old)*FactCalltoterg(box1,box2,en1n,en2n)If(acceptancerule=.true.)thenreturntoprogrammainElsereturntoendif

endSubroutineVolume(3)Transfe子程序

Transfe子程序的功能是在随机选择的子区域中随机地选择粒子，并使其在两子区域间进行转移或交换，使各组分在两相中的化学势渐趋于相等，程序简化如下：

SubroutineTransfer(Rn(x,y,z))If(ranf().lt.0.5)thenIn=1out=2ElseIn=2out=1endif

Rn(x,y,z)=ranf*box(in)&

.or.Rn(x,y,z)generatedfromGBMCCallenergy(M(Rn(x,y,z)))

NPTNPT

If(pexchange=.true.orptransfer=.true.)than

#计算box-I与box-II的初始能量#改变box-I与box-II的体积

#计算box-I与box-II的新能量

SubroutineVolume#对体积移动的合理性进行判断

#决定哪一盒子插入或拆出分子#插入分子的空间构型

#计算分子作用能

Mexchange=.true.

33

重庆大学硕士学位论文3模拟过程及程序

end

Else

ReturntotrytoinsertmoleculeEndif

SubroutineTransfer

#分子交换接受规则

3.4.3马尔科夫循环框图

GEMC模拟中的马尔科夫循环框图如图3.5所示。

34

重庆大学硕士学位论文3模拟过程及程序

STEP=1TONSTEPI=1TON计算I分子和其它分子的势能及维里分子I移动，产生新的位形计算新位形中I分子和其它分子相互作用的势能及维里TOVRLAP=.TRUE..F计算ΔU,ΔW,ΔβUFΔβU<75.0exp(-ΔβU)>RANF(DUMMY)FFTΔU≤0.0TT计算U,P平均值接受I的新位形,计算U,P的瞬时值计算接受比,调整移动距离及角度改变量T计算化学势STEP>104F图3.5GEMC模拟中的马尔科夫循环框图Fig3.5MarkovcycleinGEMCsimulation35

重庆大学硕士学位论文4模拟结果与分析

4模拟结果与分析

4.1水/甲烷和水/二氧化碳的二元体系的模拟

水、甲烷和二氧化碳在化工生产中都是比较常见的，而且水是化工生产过程中最常用的溶剂以及许多化学反应的介质，因此，研究水/甲烷和水/二氧化碳的二元体系的热力学特性具有比较重要的意义。本文采用NPT-GibbsEnsembleMonteCarlo(NPT-GEMC)方法，采用固定点电荷水分子模型、BuckinghamExp-6甲烷分子模型和EMP2二氧化碳分子模型，计算汽液相平衡组成，对水/甲烷和水/二氧化碳的二元体系在一定温度和较宽的压力范围下的汽液相平衡进行了系统的模拟计算。水/甲烷体系的模拟结果与文献数据的对照情况如图4.1，4.2和表4.1所示,水/二氧化碳体系模拟结果与文献数据的对照情况图4.3，4.4和表4.2所示，表与图中的“simu”与“exp”分别表示模拟值与文献值。从对照结果中可以看出：对于水/甲烷体系，在200bar至500bar压力下，模拟结果基本与实验数据吻合，但是，在500bar以上的高压状态下由于水和甲烷发生缔合，模拟结果与实验数据有些偏差；对于水/二氧化碳体系，模拟结果与文献数据吻合较好。因此，所采用的模拟方法和势能函数模型、LB混合规则及分子(基团)参数对这两个体系的进行汽液相平衡模拟是合适的，能够在较宽的压力范围内比较准确地预测水/甲烷和水/二氧化碳的二元体系的汽液平衡性质，为以后进一步研究预测更为复杂的水的二元混合物的性质打下良好的基础。

表4.1

水/甲烷体系的模拟结果

Table4.1ThesimulationresultsforWater/Methane

P(bar)

x(CH4)

simu

T=423K

100

25050010001002505001000

0.0008120.001450.003280.004480.001550.007540.01160.0141

0.9460.9740.9830.991

T=523K0.5060.7730.8670.934

98.1196.1588.4980.7

0.002500.01170.01410.0151

0.5300.7430.8720.91

98.1392.3588.4980.7

0.00110.00310.00460.0056

0.940.9630.9830.985

y(CH4)

P(bar)

x(CH4)

exp

y(CH4)

36

重庆大学硕士学位论文4模拟结果与分析

12001000800simu(423K)exp(423K)simu(523K)exp(523K)P(bar)600400200000.0020.0040.0060.008x(CH4)0.010.0120.0140.016Fig4.1水/甲烷的平衡组成(低组成段)

Fig4.4Thepressure-compositiondiagramofWater/Methaneatlowerconcentrations

12001000800P(bar)simu(423K)exp(423K)simu(523K)exp(523K)60040020000.50.60.7y(CH4)0.80.91Fig4.2水/甲烷的平衡组成(高组成段)

Fig4.2Thepressure-compositiondiagramofWater/Methaneathigherconcentrations

37

重庆大学硕士学位论文4模拟结果与分析

表4.2水/二氧化碳体系模拟结果

Table4.2ThesimulationresultsforWater/Carbondioxide

P(MPa)2.557.512.52.557.512.5

x(CO2)

simu

T=348.15K

0.003400.007510.1410.1610.003310.006600.009120.0132

0.9690.9850.9870.994T=393.15K0.80.9470.9590.975

98.1196.1588.4980.7

0.002800.007510.01010.0122

0.8600.9390.9610.973

25812

0.003810.007800.01530.0172

0.9750.9900.9930.994

y(CO2)

P(MPa)

x(CO2)

exp

y(CO2)

14simu(348.15K)12simu(393.15K)exp(348.15K)exp(393.15K)10p(MPa)82000.0020.0040.0060.0080.01x(CO2)0.0120.0140.0160.018Fig4.3水/二氧化碳的平衡组成(低组成段)Fig4.3Thepressure-compositiondiagramofWater/Carbondioxideatlowerconcentrations

38

重庆大学硕士学位论文4模拟结果与分析

141210exp(348.15K)p(MPa)simu(348.15K)simu(393.15K)8200.84exp(393.15K)0.860.880.90.92y(CO2)0.940.960.981Fig4.4水/二氧化碳的平衡组成(高组成段)

Fig4.4Thepressure-compositiondiagramofWater/Carbondioxideatlowerconcentrations

4.2正构烷烃二元混合物汽液相平衡的模拟

碳氢化合物在化学工业中占有着相当大的比重，正构烷烃是构造最简单的碳氢化合物，因此，利用分子模拟方法对碳氢化合物进行研究就由正构烷烃开始。本文采用NPT-GibbsEnsembleMonteCarlo(NPT-GEMC)方法，分别采用BuckinghamExp-6分子模型和L-J分子模型进行模拟计算，模拟计算了乙烷/庚烷和辛烷/十二烷的二元体系的汽液相平衡数据，并将所得到的结果与文献值进行对比。乙烷/庚烷二元体系的相平衡数据如表4.3和图4.5所示；辛烷/十二烷二元体系的相平衡数据如表4.4和图4.6所示。表与图中的“simu”与“exp”分别表示模拟值与文献值。从结果中可以看出：采用两种模型模拟计算出来的数据与文献值都能吻合，但是采用BuckinghamExp-6分子模型比采用LJ分子模型进行模拟所得到得结果更好地与文献值吻合。因此，所采用的模拟方法对这两个体系的进行汽液相平衡模拟是合适的，能够比较准确地预测的体系的汽液相平衡性质；模拟中采用BuckinghamExp-6分子模型更能准确的对简单正构烷烃的体系的二元相平衡进行预测；所得到的结果有助于以后利用分子模拟方法进一步研究预测更复杂的碳氢化合物的多元体系的相平衡的性质。

39

重庆大学硕士学位论文4模拟结果与分析

表4.3乙烷/庚烷体系的模拟结果

Table4.3ThesimulationresultsforEthane/Heptane

P(bar)

x(C2H6)

simu(EXP-6)

T=366K

1020355560P(bar)8.218.532.248.661.1

0.1380.2400.4120.6150.2x(C2H6)

exp

0.09210.2410.3980.5630.675

0.7830.9370.9450.9500.9530.8300.9250.9420.9470.950y(C2H6)

0.0990.2380.3960.6150.4

0.8530.9300.9370.9420.945

y(C2H6)

x(C2H6)

y(C2H6)

simu(L-J)

706050P(bar)EXPsimu(EXP-6)simu(L-J)40302010000.20.40.60.81x(C2H6)Fig4.4乙烷/庚烷的平衡组成Fig4.4Thepressure-compositiondiagramofEthane/Heptane

40

重庆大学硕士学位论文4模拟结果与分析

表4.4辛烷/十二烷体系模拟结果

Table4.4ThesimulationresultsforOctane/Dodecane

P(bar)

x(C8H18)

simu(EXP-6)

T=366K

360380400410420

0.6210.4110.1690.1120.072x(C8H18)

exp

0.5870.3420.1980.1100.048

0.9700.9110.7980.6010.3810.9780.9420.7790.6130.413y(C8H18)

0.4780.2430.1100.0470.028

0.9710.8630.6250.4200.169

y(C8H18)

x(C8H18)

y(C8H18)

simu(L-J)

P(bar)3623813910418

EXPSIMU(EXP-6)SIMU(L-J)430420410400T(K)39038037036035000.20.4x(C8H18)0.60.81s

Fig4.5辛烷/十二烷的平衡组成Fig4.5Thetemperature-compositiondiagramofOctane/Dodecane

41