对于较大分子,计算势的时间会随位点数目的平方而递增,从而使得计算效率大为降低.

描述非球形分子的Gay-Berne势和电多极展开势(EMP)
任何一个好的模型必须具备两个主要特征,一方面数学上要简单并易于处理,另一方面能够对模拟体系做较好的物理描述.基于这两方面的考虑,Berne和Pechukas15近似地将分子视为绕其主轴旋转的椭球体,得到了由高斯函数描述的经验势.在对非球形分子的分子模拟中,短距离吸引和排斥相互作用可通过如下方法来完成,即在分子内定义多个位点,每两个位点间的范德华相互作用势由Lennard-Jones势描述.需指出的是对于较大分子,计算势的时间会随位点数目的平方而递增,从而使得计算效率大为降低.为了解决这一问题,Gay和Berne16对描述非球形分子间的短距离排斥和相互吸引作用的高斯重叠势进行了修正,把分子处理成软的、单轴的椭球体,把GB位点放置在椭球体的质心上,用单位点势替代多位点势,这样就得到了GB势.在惯性系中,椭球体在位形空间中由其质心坐标和三个欧拉角描述,对应于由一套GB参数来描述.目前,GB势在描述液晶体系中得到了广泛应用.17-20为了更好地控制GB势的软硬度,Brene和Pechukas建议增加参数来控制经验势的柔性,由此Kabadi等21,22引入了dw参数.


当两个分子相同或者其中一个分子是球形分子时,势函数用Berne和Pechukas的原始形式表示.如果两个分子都是球形的,势函数将进一步简化为Lennard-Jones势.
长距离静电相互作用可以由电多极展开势来表示.在大于分子尺寸的距离处,分子的电势可准确地通过对其质心的电多极展开势来表示23
M=()
q,μx,μy,μz,Θxx,Θxy,…,Θzz(12)
其中q、μ、Θ分别表示电荷、偶极矩和四极矩.电多极势可表示为电荷-电荷相互作用势、电荷-偶极相互作用势、电荷-四极相互作用势等项之和.


基于粗粒化模型对有机溶剂的分子动力学模拟
许佩军1
唐媛媛1,2张静1张知博1王昆1
邵颖3
沈虎峻2,*
毛英臣1,*
(1辽宁师范大学物理与电子技术学院,辽宁大连116029;2
中国科学院大连化学物理研究所,辽宁大连116023;
3
大连海事大学物理系,辽宁大连116026)
摘要:在基于Boltzmann分布对四种基本构象进行MonteCarlo取样后,通过与全原子模型的范德华势比较
得到了Gay-Berne(GB)参数.又在对用量化计算得到的分子体系的电势进行电荷、偶极矩和四极矩的拟合后,得到了电多极展开势(EMP)参数.利用得到的粗粒化参数,基于粗粒化模型,对CHCl3及四氢呋喃(THF)两种有机溶剂进行了分子动力学模拟(MDS),并将结果同全原子模拟进行了比较.计算结果表明用粗粒化模型从整体上能重复全原子模型的模拟结果,但在某些细节的计算与全原子模型有偏差,其原因可能是目前工作仅考虑了单位点情况,为此今后在对具有复杂结构的分子进行粗粒化模拟时还应考虑合理放置及增加相互作用位点.关键词:
粗粒化模型;Gay-Berne势;电多极展开势;径向分布函数;分子动力学模拟
中图分类号:
O645;O641
MolecularDynamicsSimulationofOrganicSolventsBasedonthe
Coarse-GrainedModel
XUPei-Jun1
TANGYuan-Yuan1,2ZHANGJing1ZHANGZhi-Bo1WANGKun1
SHAOYing3SHENHu-Jun2,*MAOYing-Chen1,*
(1SchoolofPhysicsandElectronicTechnology,LiaoningNormalUniversity,Dalian116029,LiaoningProvince,P.R.China;
2
DalianInstituteofChemicalPhysics,ChineseAcademyofSciences,Dalian116023,LiaoningProvince,P.R.China;
3
DepartmentofPhysics,DalianMaritimeUniversity,Dalian116026,LiaoningProvince,P.R.China)Abstract:ToobtainGay-Berne(GB)parameters,wecarriedoutMonteCarlosamplingoffourreferenceconfigurationsbasedontheBoltzmanndistribution.AftercomparingwiththevanderWaalspotentialwithintheall-atommodelweobtainedtheGBparameters.Alsobyfittingthecharge,dipole,andquadrupolewiththeelectricpotentialobtainedfromquantumchemicalcomputationswithGaussian03weobtainedtheelectricmultipolepotential(EMP)parameters.WiththeGB-EMPparameterswethencarriedoutmoleculardynamicssimulations(MDS)forCHCl3andtetrahydrofuran(THF)basedonthecoarse-grained(CG)model.Comparedwiththeall-atommodel,theCGmodelcanreproducethesimulationresultsonthewhole,buttherearesomedeviationsinthesimulationsinsomedetails.Thereasonisthatweonlytakeoneinteractionsiteintoaccountinthiswork.Therefore,formorecomplicatedmoleculesitisnecessarytotaketheplacementoftheinteractionsitesintoaccount.Additionally,themulti-sitessituationisalsoconsideredintheMDSwithintheframeofthecoarse-grainedmodel.KeyWords:
Coarse-grainedmodel;Gay-Bernepotential;Electricmultipolepotential;Radial
distributionfunction;Moleculardynamicssimulation
[Article]
www.whxb.pku.edu.cn
物理化学学报(WuliHuaxueXuebao)
ActaPhys.-Chim.Sin.2011,27(8),1839-1846
AugustReceived:February28,2011;Revised:May9,2011;PublishedonWeb:June14,2011.∗
Correspondingauthors.MAOYing-Chen,Email:myc@lnnu.edu.cn;Tel:+86-411-82158367.SHENHu-Jun,Email:hshen@dicp.ac.cn;Tel:+86-411-84379875.
TheprojectwassupportedbytheFundamentalResearchFundsfortheCentralUniversities,China(2009QN069).中央高校专项资金优秀青年教师基金(2009QN069)资助项目
ⒸEditorialofficeofActaPhysico-ChimicaSinica
1839

   2014全国一级建造师资格考试备考资料真题集锦  建筑工程经济  建筑工程项目管理  建筑工程法规  专业工程管理与实务   
 

Vol.27
ActaPhys.-Chim.Sin.2011
1引言
分子动力学模拟(MDS)已经成为了研究复杂分子体系的重要工具,尽管一些计算方法和技术手段已经取得了突破,1-3但当前基于全原子模型对某些较大体系的动力学模拟时间尺度仍然非常短,并不能满足与实验比较的迫切需要.4-6已知在利用经验力场进行的全原子分子动力学模拟主要包括两部分:一是对分子的构象进行统计分析;二是结合相应的构象分析体系的能量及受力情况.7-9而对某一分子体系,由于非键结的范德华相互作用项和静电相互作用项的数目远远大于键长、键角和二面角相互作用项的数目,因此对范德华相互作用和静电相互作用的计算将要耗费大量的计算时间.为了简化对这两种相互作用的计算,研究人员已提出了多种粗粒化模型(CG模型).10-14CG模型的主要特点在于它将研究体系中的分子或分子的某一部分近似看成一个原子团簇,该原子团簇内部信息(键长、键角、二面角等)被屏蔽,即该团簇被处理为没有内部结构的“粒子”,该粒子对应于全原子模型中的单个原子.而在利用CG模型进行分子动力学模拟之前还必须获得描述这些粒子间范德华相互作用和静电相互作用的力场参数,应当指出这一部分工作也是构建CG模型的关键部分.CG模型结合了日益发展的计算机水平,使得对更大更复杂体系的长时模拟成为可能.
该方法的主要缺陷就在于它的计算精度要低于基于全原子模型的分子动力学模拟.10此外还需指出的是大多数的CG模型仅对某些特殊体系适用,而不能应用于其他体系,即模型的可移植性有待提高.产生这一问题的主要原因在于多数CG模型将粗粒化粒子看成是球形,这样导致了这些模型没能更好地描述分子间的范德华相互作用,以及未能对由此而导致的电势空间分布改变的静电相互作用作出更好的描述.为了克服这些缺陷,改进CG模型,Golubkov和Ren14提出了在计算粗粒化粒子间相互作用时,必须考虑由粗粒化粒子的形状而引起的Gay-Berne(GB)效应及电多极矩效应,这些改进在他们对水、苯及甲醇三种构象相对简单的分子体系的动力学模拟中得到了较理想体现.此外也应认识到的是他们工作的基石就在于较好地拟合得到了描述上述三种分子的粗粒化力场参数.为了检验如何利用CG模型对复杂分子体系进行高效计算,本文采用类似方法对CHCl3及四氢呋喃(THF)分子
的两种构象在相对复杂的有机溶剂中进行了分子动力学模拟.
2
计算与模拟方法
2.1
描述非球形分子的Gay-Berne势和电多极展开势(EMP)
任何一个好的模型必须具备两个主要特征,一方面数学上要简单并易于处理,另一方面能够对模拟体系做较好的物理描述.基于这两方面的考虑,Berne和Pechukas15近似地将分子视为绕其主轴旋转的椭球体,得到了由高斯函数描述的经验势.在对非球形分子的分子模拟中,短距离吸引和排斥相互作用可通过如下方法来完成,即在分子内定义多个位点,每两个位点间的范德华相互作用势由Lennard-Jones势描述.需指出的是对于较大分子,计算势的时间会随位点数目的平方而递增,从而使得计算效率大为降低.为了解决这一问题,Gay和Berne16对描述非球形分子间的短距离排斥和相互吸引作用的高斯重叠势进行了修正,把分子处理成软的、单轴的椭球体,把GB位点放置在椭球体的质心上,用单位点势替代多位点势,这样就得到了GB势.在惯性系中,椭球体在位形空间中由其质心坐标和三个欧拉角描述,对应于由一套GB参数来描述.目前,GB势在描述液晶体系中得到了广泛应用.17-20为了更好地控制GB势的软硬度,Brene和Pechukas建议增加参数来控制经验势的柔性,由此Kabadi等21,22引入了dw参数.
两个粗粒化粒子间的GB势可表示为
UGB()u
̂i,ûj,r̂ij=4ε()ûi,ûj,r̂ij×éë
êêêêæèççöø÷÷dwσ0rij-σ()ûi,ûj,r̂ij+dwσ012
-ùû
úúúúæèççöø÷÷dwσ0rij-σ()ûi,ûj,r̂ij+dwσ06
(1)
上式中的距离参数σ()u
̂i,ûj,r̂ij可取如下形式σ()ûi,ûj,r̂ij=σ0éëêêêê1-ìíî
ïïχα2()ûi·r̂ij2+χα-2
()ûj·r̂ij2
1-χ2()ûi·ûj2-ùû
ú
úúúüýþïï2χ2
()ûi·r̂ij()ûj·r̂ij()ûi·ûj1-χ2()ûi·ûj2
-12
(2)
1840

 

No.8
许佩军等:基于粗粒化模型对有机溶剂的分子动力学模拟
其中
σ0=di2+dj2(3)χα2
=li2-di2
li2
+dj2
(4)χα-2
=lj2-dj2lj2
+di2(5)χ2
=()li2-di2()lj2-dj2()l
j
2
+di
2
()
l
i
2
+dj
2
(6)
其中l和d分别描述分子的长和宽,这样分子的形状就可以用任意的棒状、盘状或者球状来表示.
(1)式中的阱深参数取如下形式
ε()u
̂i,ûj,r̂ij=ε0ε1ν()ûi,ûjε2μ
()ûi,ûj,r̂ij(7)
其中μ和ν均为可调指数,借鉴Golubkov和Ren14的工作,本文分别将其取值为2.0和1.0.而ε1和ε2可分别表示为
ε1()ûi,ûj=éëùû
1-χ2()u
̂i·ûj2
-12
(8)
ε2()ûi,ûj,r̂ij=1-ìíî
ï
ïχ′α′2()ûi·r̂ij2+χ′α′-2()ûj·r̂ij2
1-χ′2()ûi·r̂ij2
-üýþ
ï
ï2χ′2()ûi·r̂ij()ûj·r̂ij()ûi·ûj1-χ2()ûi·ûj2
(9)
(9)式中的χʹ和αʹ分别取如下形式χ′=
1-()
εEεS1μ1+()
εEεS1μ
(10)
α′2
=éëùû
1+()εEεS1μ
-1
(11)
其中,ε0描述了交叉结构的阱深,εE和εS分别表示
尾对尾和边对边结构的阱深.当两个分子相同或者其中一个分子是球形分子时,势函数用Berne和Pechukas的原始形式表示.如果两个分子都是球形的,势函数将进一步简化为Lennard-Jones势.
长距离静电相互作用可以由电多极展开势来表示.在大于分子尺寸的距离处,分子的电势可准确地通过对其质心的电多极展开势来表示23
M=()
q,μx,μy,μz,Θxx,Θxy,…,Θzz(12)
其中q、μ、Θ分别表示电荷、偶极矩和四极矩.电多极势可表示为电荷-电荷相互作用势、电荷-偶极相互作用势、电荷-四极相互作用势等项之和.本文在分子的质心上放置了一个电多极展开位点,同时选取
了和GB位点相同的惯性系.
两个多极位点间的相互作用可通过电多极展开势表示为UEMP=MtiTijMj
(13)
其矩阵形式可参见文献.14,24本文对静电相互作用截断值的处理方法,采用了TINKER软件的标准方法.25-27
此外还需指出的是,考虑到由于粗粒化粒子具有各向异性的形状而使得在进入排斥的范德华相互作用区域前电多极展开方法将给出不正确的相互作用势,类似于文献,14,28本文同样加入了一个阻尼函数以确保对短距离多极相互作用势的精确描述.
2.2粗粒化模型力场参数的确定2.2.1Gay-Berne(GB)参数的拟合
本文利用粗粒化模型进行分子模拟的前提是要获得GB-EMP模型的参数,为此我们首先对研究体系进行了合理取样选取,然后利用描述有机分子的AMBER力场(GAFF),29,30通过与全原子的范德华势进行比较后拟合得到了GB参数.选取GAFF力场是由于它的函数形式简单,原子类型有限,还可利用实验和理论模型计算获得力常数和部分原子电荷.在研究中,我们发现构象的选取会直接影响GB参数的拟合,进而决定GB势是否能准确地模拟全
表1氯仿和THF分子的四种基本构象
Table1FourreferenceconfigurationsofCHCl3and
THFmolecules
face-to-face
side-by-side
T-shape
C-shape
CHCl3
THF
1841

 

Vol.27
ActaPhys.-Chim.Sin.2011
原子的范德华势,因此选择可以反映真实分子液体环境的构象是十分重要的.由于在实际的溶液环境中,面对面(或者尾对尾)、边对边、T形和交叉结构这四种构象出现的几率相对较大,17,19另外考虑到在真实的液体环境中分子的构象趋于取能量最低的构象,所以本文在筛选用于拟合GB参数的构象时,对上述四种构象分别基于Boltzmann分布采用MonteCarlo方法获得到分子的构象集.这样才可确保在一定的能量范围内,构象数目随能量的升高而减少.表1给出了CHCl3和THF这两种分子的四种基本构象.
图1及图2给出了氯仿和THF分子的四种基本构象在各能量区间内的构象数目.从图中我们可以清楚地看到所选的分子构象在较低的能量范围内
较多,且构象数随着能量的升高而逐渐减小.
在参数拟合的具体过程中,5个自由参数(l,d,ε0,εE/εS,dw)由遗传算法拟合得到.表2给出用上述方案拟合得到的氯仿和THF分子的GB参数.2.2.2EMP的参数拟合
在拟合粗粒化模型的EMP参数时,本文首先应用Gaussian03软件包31的MP2方法采用6-311G**(2d,2p)机组计算得到了两种分子的电势,然后利用GDMA程序32对电势进行电荷、偶极和四极的拟合,
图2THF分子的构象数目随能量的变化
Fig.2ConfigurationnumbersofTHFmolecularchangewiththeenergy
Thedashlinesareonlyguidestotheeyes.
图1氯仿分子的构象数目随能量的变化
Fig.1
ConfigurationnumbersofCHCl3molecularchangewiththeenergy
Thedashlinesareonlyguidestotheeyes.
表2氯仿和THF分子的GB参数
Table2
GBparametersofCHCl3andTHFmolecules
CHCl3THF
l/nm0.2470.313
d/nm0.4020.403
ε0/(kJ·mol-1)3.1823.977
εE/εS2.291.37
dw0.660.61
1842

 

No.8
许佩军等:基于粗粒化模型对有机溶剂的分子动力学模拟
最后得到了两种分子粗粒化模型的EMP参数.表3中给出了氯仿和THF的EMP参数.2.3模拟体系与动力学模拟方法
为了检验本文对氯仿和THF分子粗粒化模型GB-EMP参数的拟合,我们利用粗粒化模型对两种有机分子溶剂进行了粗粒化的分子动力学模拟.为了在计算的过程中可任意地在粗粒化模型和全原子模型间进行转换,方便容易地比较CG模型和全原子模型的结果,故本文在计算中采用了TINKER软件包.此外还基于如下考虑,即TINKER也包含了对粗粒化粒子GB-EMP能量、力和力矩部分的计算,并可在惯性系中非常容易地整合GB和EMP两部分的计算结果.
氯仿和THF有机溶剂的立方体盒子由Materi-alsStudio(MS)5.0软件生成,其边长均为3nm.对GB和EMP两部分的截断半径都取为1.2nm(这考虑到在通常情况下,要求截断半径值不大于盒子边长的一半).CG模型和全原子模型的分子动力学模拟都是采用NVT系综在室温(298K)下进行的.基于全原子模型的分子动力学模拟的计算步长取为1fs.文献14对考虑了GB-EMP效应的CG模型的时间步长作了实验,推荐CG模型的计算步长可取5-20fs.考虑到由于自由度的减小,以及分子内高频运动的“屏蔽”,在本文我们取时间步长为5fs.在对运动方程进行数值化处理时,两种模型的计算都采用了TINKER推荐的BERMAN积分算法.25
3
结果与讨论
3.1
对范德华相互作用的分析
在图3中,我们比较了分别利用CG模型和全原子模型计算得到的氯仿分子的范德华相互作用势,可以看到CG模型结果与全原子模型结果总体符合得较好.对面对面及交叉结构,两种方法得到的作用曲线几乎完全相同;而对于边对边及T型构象,可以看出相比于全原子模型,粗粒化模型得到了较小
的阱深和距离参数,这一现象对T型构象尤为明显.此外,基于粗粒化模型对四种构象的计算结果整体趋势较好地符合了Gay及Berne16的结论.
通过将分子构象的取样图(图1)与范德华相互作用势能曲线图(图3)进行比较,可以发现分子构象数最多的能量区域正好对应于范德华相互作用势能曲线图的最低点区域,这也印证了分子构象都趋于位于能量最低的构象状态的结论.19
图4给出了基于两种模型对THF分子范德华相互作用的比较.从图中可以看出,相比于结构简单的氯仿分子,对THF分子进行粗粒化模型的计算结果与全原子结果有较大的偏差.这表现在对应于四种基本构象的阱深参数和距离参数都一致性地小于全原子模型的计算结果.从图3和图4中,我们也观察到两种分子的交叉结构和边对边结构的范德华作用曲线非常接近,几乎重合,这一现象很好地支持了文献16的结论.在对粗粒化参数做多次拟合后,计算结果仍然如此,由此我们认为出现这一现象的主要原因是我们仅考虑了单位点情况.这也建
表3氯仿和THF分子的EMP参数
Table3
EMPparametersofCHCl3andTHFmolecules
CHCl3
THF
Charge/C
000000
1021Dipole/(C·nm)
2.5230.5161.3731.8040.0010.000
1022Quadrupole/(C·nm2)-0.517-0.144-1.718-4.660-0.0030.000-0.1440.155-0.352-0.0033.5150.163-1.718-0.3520.3620.0000.1631.145
Inertia/(g·nm2
·mol-1)
15560.915561.530030.07031.87116.712475.5
图3
基于粗粒化模型和全原子模型对氯仿分子的范德华
相互作用曲线的比较
Fig.3ComparisonofvanderWaalspotentialofCHCl3
basedontheCGmodelandtheall-atommodel
Thecurvesoffourconfigurations,namelyfacetoface(lines1and2),sidebyside(lines5and6),T-shape(lines3and4),andcross
shape(lines7and8),arealsodisplayed,respectively.
1843