分子动力学入门第三章Word文件下载.docx
- 文档编号:22393657
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:16
- 大小:69.46KB
分子动力学入门第三章Word文件下载.docx
《分子动力学入门第三章Word文件下载.docx》由会员分享,可在线阅读,更多相关《分子动力学入门第三章Word文件下载.docx(16页珍藏版)》请在冰豆网上搜索。
在分子动力学模拟过程中,系统在某个状态具有一特定的密度,温度,压力(这些量的计算会在下面介绍):
模拟材料的相图可以被研究。
然而,我们怎样让系统到达我们想要的区域?
换句话说,我们怎样控制系统?
在标准的算法中,密度的控制是由算选的盒子的体积控制的。
程序通常在开始计算时根据需要对体积有规定。
体积的变化要适度,最多只有几个百分点的变化。
压强会在跑的过程中测定,在下面会介绍。
另一方面,更多的分子动力学模拟,用户可选择压强,而盒子的体积是可以变化的(在过程中测量)。
温度的变化可以通过设备中的代码调节速度获得。
速度的Verlet算法在2.3.1中
已讨论,下面公式
v(t+=t/2)=v(t)+1/2a(t):
t
可以替代为:
v(t+:
t/2)=Tov(t)+1/2a(t)-:
YT(t)
To是想要的温度,T(t)是模拟过程中瞬时温度。
这种改变说明我们不在遵循牛顿方程,并且总能量也不再守恒。
该过程不能收集重要数据:
这种控制温度的方法只能用于把系统从一个状态转换到另一个状态。
我们必须要等到系统在一个恒定的能量下达到平衡才能开始收集数据。
3.3平衡
每次当系统的状态发生改变时,系统都会破坏平横一会儿。
我们这里指的是热力学平衡。
由鉴于此,这就意味着系统状态的指示量一一大多数下面都会介绍一一不是稳定的(在一个固定值周围扰动),而是relaxing向一个新值(在一个随时间变化不大的值周围扰动)。
系统状态的变化可能是我们自己或自动变化引起的。
可以通过我们改变模拟的参数而引起变化一一例如改变温度,密度一一从而扰动系统,接着我们等待到达一个新的状态。
自动变化,例如当体系进行相变时,从而从一个平衡状态到达另一个平衡状态。
在所有的情况,我们都希望在我们进行测量之前系统要达到平衡。
一个物理量A一般接近它的平衡值是时间的指数:
A(t)=Ao+Cexp(-t/.)
A(t)表示通过一段很短时间摆脱瞬间扰动过程中的一个物理量的平均值,不是很
长时间的扰动.•是扰动持续的时间(再次平衡所需的时间,很短.)我们有时需要•有几百步的时长的数量级,使A(t)收敛于A。
,从而直接测量平衡值。
在另一方面.太大超过了我们模拟的数量级,例如,一秒的数量级。
在这种情况,我们不在跑的过程中能看能任何relaxation,分子动力学不能用于这种情况。
在上面两种情况之间的情况,我们也许能看到扰动,但是我们要等很长时间,A(t)才收
敛于A。
在这种情况,即使最终的点任然离A。
很远,我们任然可以利用可利用的数据按照上式估计A。
3.4观察原子(lookatatom
一种最简单的”探针“(至少在概念上)我们用来了解我们的系统正在发生什么是”
“lookingatit“。
我们可以指定一个半径给原子(例如原子间对势变得强烈排斥时的距离的一半),原子看成球,球的半径,例外有个坐标从一个观察到的特定的点构建系统的“图像”。
最终的图像可以在图形屏幕上观察到或打印出来。
该过程一般都是利用程序解决的。
期中的一个程序叫“Ballroom”,是免费应用
的。
除了别的以外,该程序可以标记单独原子的特性(例如,扩散,势能,坐标等)成不同颜色,可以形象的表示不同特性在空间的分布。
在一些情况下,例如分析复杂结构的原子布置,目测是不行的。
在例外的一些情况下,给的信息很少,有用的数据必须通过处理原子的坐标或速度来提取,例如获得相关函数,接下来会讨论。
3.5测量简单的统计量
测量分子动力学物理量通常指的是通过系统的轨迹求物理量的时间平均值。
物理
性能通常是离子的位置和速度的函数。
因此,例如,我们可以概括一个在时间t
的物理性能的值:
A(t)=f(ri(t),,,\(t),Vi(t),,,Vn(t))
则平均值为:
1Nt
<
A>
=—A(t);
Ntta
t指的是从第一步到总的Nt步所跑的时间。
在实际情况下有两种等价的方法来计算物理量:
(1)A(t)在跑的过程中每步都被MD程序计算。
A(t)每步都会更新。
当跑
完后,平均值可以通过除以总的步骤立刻获得。
当物理量计算简单或相当重要,优先考虑这种方法。
一个例子就是求系统的温度。
(2位置(也可能是速度)在“轨迹文件中”周期性更替,当程序跑动时。
一
个独立的程序,在模拟程序跑完后,会处理轨迹文件并计算想要的物理量。
这种
方法对硬盘的需求非常大:
用64位精度记录每步每个粒子的数据需要48个字节。
但是,这种方法经常用于计算复杂体系,或结合不同时间段的动态相关性,或者物理量依赖于当MD程序跑动时其他不知道的例外的参数。
在这种情况下,模拟只跑动一次,但是结果轨迹可以反复处理。
接下来将要讨论最常见的要测量的物理量。
3.5.1势能
势能V的平均值是通过求它的瞬间值的平均值获得的,通常这是直接获得的,并同时计算力。
例如,two-body相互作用
V(t)…"
(|ri(t)-rj(t)|)
iji
尽管完全不需进行时间积分(力的计算是需要的),但是必须要验证能量的守恒。
在MD模拟时这是一个重要的检查。
3.5.2动能
瞬间的动能公式:
K(t)=1/2、mi[Vi(t)]2
i
这是非常容易计算的。
3.5.3总能量
总能量E=K+V在牛顿力学中这是一个守恒的物理量。
然而,在实际中普每步中都要计算它的值,这样做是为了检查它的确是随时间变化是守恒的。
换句话说在
分子动力学过程中能量在势能和动能之间不断的转换,引起K(t)和V(t)的
扰动,但是他们的和是固定不变的。
在实际中总能量是有小的扰动的,具体的扰动的量不超过10=这种扰动通常是由时间积分的算法的误差引起的,如果每步的时间认为过多,可以通过减少每步的时间而减少误差的数量级。
Ref充分的分析了用不同的时间积分算法引起的总
能量扰动。
总能量的慢速漂移有时在长时间的模拟过程中被观察到。
这种漂移可能也是由于过多的厶t。
能量漂移比能量扰动更危险,因为系统的状态伴随着能量的变化,因而整个过程的时间平均值不能指单一的状态。
如果经过一个很长时间能量漂移
发生了,他是可以被阻止的,例如,把长时间划分成更小的步骤,并且恢复能量到正常值在一步与下一步之间。
一个通常的调节能量的方法是通过调节速度来改变动能。
最后的一个警告是:
我们尝试通过不断降低步长而获得非常好的能量守恒,但是
这会增加计算机计算时间。
实际情况中,会允许存在很小的能量扰动和慢的能量漂移最为一个可以接受的大的.-:
t的代价。
具体看:
[3]M.P.AlienandD.J.rildcslcy,Computersimulationofliquids,Oxford^1987.Thisisadassichowtobook,coutiuningalargenumberoftechnicaldetailsonavastarrayoftcdiiiiqucs.
[6]J.M.Haile,MoleculardynaTnicssiTnulation^Wiley,1992.Anexcellentandrecentgeneralintrodnetiontomoleculardynamics.Notasthoroughas[3inexjjlaiuingtechniques〕butverydeepincxplciiningthefuudaiiicuttLk,including^philosophicaPaspect他Probablythebeststartingpoint,AgoodaiuuuiitufFortran77codeisincluded.
3.5.4温度
温度T直接与动能有关,根据能量均分定律,在每个自由度上动能为■-bT/2:
K=3/2(N■■-BT)
因此温度可以直接从平均动能K获得。
为了实际的需要,通常也定义一个瞬间
温度T(t),与瞬间动能温度K(t)成正比,具有一个类似上面的公式。
3.5.5热量曲线
一旦总能量E(也叫内能)和温度T在不同的过程和不同的热力学状态下被
测定了,我们可以构建一个热量曲线E(T)。
这对监测相转变非常有用的,E(T)
曲线上的跳跃暗示一级相变,E(T)导数上的一个跳跃暗示二级或更高级相变。
然而,这种策略不适合准确估计熔点Tm,在3.6节讨论。
最普遍的一阶相变就是融化,非常容易通过模拟观察到。
当系统放弃晶体结构变成无序的液体,我们可以看到E(T)的一个跳跃,与潜在的融化热一致。
这通
常发生在温度比真实的Tm高的模型,由于滞后因子:
需要一个液相的种子。
只有但液相种子包含少量的原子形成后,液体开始以牺牲固体为代价成长。
在通常的分子动力学时间范围内,我们必须“overshoot”并且设置温度比Tm高20-30%,为了观察到融化的发生.
3.5.6均方偏移
按照定义原子在模拟过程中的均方偏移很容易计算
MSD=v|r(t)-r(0)|2>
….>
在这里表示所有原子的平均(或者所有原子在一个给定的亚型)。
必须注意不要认为粒子的“jumps”,当用边界条件时原子折回盒子里,做为扩散的贡献。
MSD包含了原子扩散的信息。
如果系统是一个固体,MSD在一个有限值附近,
如果系统是液体,MSD会随时间线性增加。
在这种情况下,可以用于表征系统的扩散曲线斜率,即扩散系数D:
12
D=limv|r(t)-r(O)|2>
96t
公式中的6,如果在二维空间必须替代为4.
3.5.7压力
在分子动力学过程中测量压力是依据Clausiusvirial函数
N
W(几,,,「N)吃ri*FiTOT
i丄
fTot表示作用在原子i上面的力的矢量和。
通过分子动力学轨迹,可以统计vw>
的平均值:
1tN
W>
=lim-d「二"
ri(■)*miri()
经常利用牛顿定律。
部分积分
!
ftN
0〉=一丿鑒7/dT工呵冋Of・
这是平均动能的两倍。
按照统计上的能量均分定律
{W}=-DNkHT
(2)
D是系统的维度2或3,N是系统中粒子的数目,kB是波尔兹曼常数
现在,我们可以认为一个粒子上面的力可以分为两个部分:
玲0丁=匕十F『XT
Fi是内力(来源于原子间的相互作用),卩匚已乂丁是外部力来源容器的墙,如果粒子
被包围在一个平行六面体中,各边为Lx,Ly,Lz,体积为VLxLyLz,并且
坐标原点在一个角上,部分vWExT>
来源于容器,可以利用下式计算
(yyEXT〉=S(-厂―)+Zy(-PZXL.)+S(-叫Lj=-DPV
-PLyLz,是yz墙沿x方向作用在x=Lx上的外部力fExt,因此上面的公式
(2)可以写成:
1/N
PV=NkBT+万〈工叶巧
(3)
\i=l
这就是重要的virial公式。
在模拟过程中除了压强P所有的物理量都很容易得到,因此上式(3)构成了测量P的一种方法。
注意上式(3)怎样推导成著名的理想气体方程(如果粒子间没有相互作用)。
在成对的相互作用是通过势能(r)的情况下,读者自己证明(3)式可以变成
该表达式比(3)的优点在于,很好的用于边界条件存在的情况下:
只要考虑下
rij的规定就行啦3.6测量熔点温度
均方偏移最为时间的函数可以区别固体和液体。
我们可能尝试通过升高温度直到晶体系统出现扩散现象来确定模拟物质的熔点Tm,热量曲线的跳跃表明吸热。
尽管这的确是固液转换的指标,但是在MD过程中的温度要比熔点要有不同度的高。
实际上,熔点是定义为固液共存时的温度。
然而在确少液态种子(液体从这点开始增加),超过熔点温度普遍发生。
在这个区域,系统是热力学亚稳定状态,不过在模拟的时间范围内是稳定的。
当它的力学不稳定点到达时,过度加热的大体积晶体破碎。
该点一般比熔点高20-30%。
中的一种就是设计一个大的样本含有50%的固体和50%的液体。
这种样本可以开始由人工构造,它存在固液分界面。
我们接下来试着建立固液共存的平衡态。
为了这个目的,系统在接近热量曲线的跳跃点附演化。
我们假设系统起始温度为T。
,这温度将身高通过一未知的有关能量E的函数,一般T。
与Tm是不同的。
我们假设T°
>
Tm.因此液体比固体稳定,从而固液界面开始移动,增加液体部分,
减少固体部分。
也就是说,一些材料融化了。
但这发生了,系统伴随着吸收热量,既然总的能量不变,吸热就意味动能的减少。
因而温度要下降。
驱动固液界面运动的力也会变弱,随着时间的推移,界面的位置和温度以指数级的变化达到平衡状态,这时的温度就接近Tm。
如果To<
Tm,反过来:
系统释放热量,动能增加,温度也再次以指数级变化达Tm:
不过是从Tm以下达到。
熔点依赖于系统的压强。
在最后的状态也需要测量压强来表征热力学性质。
经常,上面的步骤在一个恒定的压强技术上实现。
在这种情况下,盒子的体积会随着融化或结晶而自动变化,在恒定压强下,调节两相的不同密度。
3.7真实空间的相关性(realspacecorrelationS
真实的空间相关性有以下的形式:
A(r)A(0)>
可以由分子动力学直接获得:
我们必须从几个构象的位置和速度计算感兴趣的
A(r),构建每个构象的相关性,以及计算所有构象的平均值。
最简单的成对相关性函数g(r),本质上是密度与密度的相关性。
g(r)表示找到一对相距r的概率,指的是在同样的密度下一个统一随机分布的期望:
该函数包含了系统的结构信息。
对于晶体,它表示的是一系列峰值围绕给定原子所形成的图形上的位置。
峰值的位置和数量级是晶体结构的属性。
对于液体,g
(r)展示的峰值靠近原子与邻近原子分离距离的平均值,并且在最大距离处有明显的峰值振荡。
峰值随距离指数衰减,并且g(r)接近1.在所有情况,g(r)在距离低于一定值情况下突然消失,在这个距离排斥足够强,从而原子不能靠的太近。
经常通过g(r)计算的一个物理量是在r1与r2之间的平均原子数,
这公式允许确定所选的数目,尽管有时存在无序情况。
g(r)是一个0(N2)的算法,因而考虑优化分子动力学程序降低计算量。
如果在最远距离r处的行为不重要,因而很方便定义一个截断距离,借用快速力计算的技术,利用短程势能函数降低计算量。
应该还要注意边界条件引进的一个截断L/2,L是盒子间的最小距离。
远距离的效应通过SIZE效应消失。
3.8reciprocalspacecorrelation倒易空间相关性
结构信息也可以通过计算倒易空间获得。
最基本的方法就是密度函数
:
「(r)=7「(r-rj的傅里叶转换:
p(k)=£
^kr>
t=i
这个物理量,
础:
(对于给定的结构很容易很快的获得,)是构建静态结构因子的基
S(k)=i(p(k)p(-k))
=寺件严皿〉
这个物理量与散射实验的结果比较极其有用。
在液体时(各向异性的),S(k)
只依赖与波矢量的模,本质上是成对相关性函数的傅里叶转换:
S(k)=l+pg(k).
计算S(k)比g(r)快很多。
利用类似的观察量A(k)替代「(k),其他的相关函数也可以在倒易空间中得到计算。
必须要注意的是,在处理边界条件时,结果应该不依赖于所选的一个特许粒子的
图像。
因此我们必须利用,Lx是重复盒子的长度,这就
暗示了所选的允许的波矢量存在一个限制:
2%
nx是整数,对于y与z也一样。
因而,允许的波矢量由于边界条件是量化的。
3.9动力学分析
与蒙特卡洛相比,分子动力学是考虑系统的时间演化,因而,动力学系统的信息都是现成的。
实验的数据(例如,对于晶体的声子色散曲线)在波矢量-频率空
间可以利用。
通过MD在真实空间和时间中获得的信息,通过傅里叶转换与真实实验数据比较。
为了进行MD分析,我们经常需要程序周期性大量储存粒子的的速度和位置。
这些会产生一个很大的文件:
为了计算动力学物理量我们需要大量的存取空间。
这些数据会被其他程序利用来计算我们感兴趣的物理量。
在某
些情况,MD直接进行分析,并且把分析好的数据保存到磁盘空间。
一个典型的物理量是VanHove相关函数G(r,t),表示在时间t在位置r处找到一个原子的概率如果有一个原子在r=0在t=0时。
该物理量进行两次傅里叶转换,变成动态结构因子S(k,w)(实验报道):
S(虬3)=[/处泓』(kL況
对于静态的物理量,最方便的方法处理计算量就是利用密度函数的空间傅里叶转
换:
p(k』)二工严皿
这是单一的物理量计算非常快,通过这,我们可以构造一个时间依赖的密度-密
度相关函数(中间散射函数)通过处理所用的MD轨迹:
F(k,i)=1(p(k,t+io)p(-k,to))
这个平均是对所有可能选择的to。
最大的时间t(与频率分辨率相关)跑的时间的限制,然而周期性的数据收集一一经常是一些时间步一一决定了最大频率的范围。
动态结构因子最终通过时间的傅里叶转换获得:
S(虬s)=[dt
再次,允许的K是量子化的(边界条件的结果)。
S(k,w)在(k,w)平面的峰与分散的传播模式相关。
在不同的温度下进行计算,我们可以研究声子谱的温度效应,包括非简谐效应:
在这种情况,这种算法优于标准的晶格动力学技术。
MD最主要缺陷就是不能直接提取特征向量。
当我们处理振动模型和不同的极化情况,利用上面的方法就不行啦,必须利用:
j(k"
)=工v他严血
代替'
(k,t)
3.10退火和淬火:
MD是一个优化工具annealingandquenching:
MDasanoptimizationtool
分子动力学也可以作为优化的工具。
我们假设N个粒子有许多可能的平衡构象。
这些构像具有不同的能量,但其中的一个是最优的。
他们都是局部能量最小,每个之间都存在能垒。
寻找最优结构基于传统的最小化技术(最陡下降法,共轭梯度法等等)是存在欺骗的。
这些方法不能越过能垒,从而最终落在最近的局部最小周围。
因而我们不得不尝试几个不同的起点,在能量图中伴随不同的“attractionbasins“,是每一个达到最低。
因而最优结构就是其中的一个,具有最低能量,如果我们有足够的数量。
在分子动力学算法中(蒙特卡洛)温度提供了一条路径“flyoverthebarriers”:
具有能量为E的状态被访问的概率为exp(-E/kBT)•如果温度足够大,系统可以“see”模拟存在许多不同的最小值,并且花费更多时间在最小的那个。
温度从T慢慢降到0,系统能有很好的机会抓住最小的那个,并且停留在哪。
这就模拟的退火方法,系统在一定的温度下达到平衡,然后慢慢冷却到T=0.然而这种
方法不能确保真实的最小能过达到,但是它能找到最好的最小,并且由于不存在
先验的有关最优结构的假设,很难对产生的结构进行判断。
这种方法经常用来优化原子结构,但是有效性很一般:
给定一个目标函数zc-1,
色,。
an)依赖N个参数,我们可以认为每个参数是一个自由度,给他一个质量,用分子动力学或蒙特卡洛的算法让系统跑动进行模拟退火。
早期有关这方
面的应用可以在下面的文献中看看到(用来解决旅游商人问题)。
TW1.I丄,J.XUiUJ.kJkJFJ_|.1■JUJ)爭J*1_JLUli上J*UJJ.7\J.1fTWVJFJLJ.U1kJSJAJ-/8JJLJ.UXJ.J.J.UJ
21]SKirkpatrick,GD.Gelatt,andP.Vecchi,Science220,671(1983).
3.11其他的统计体系(otherstatisticalensembles)
我们已讨论了标准的分子动力学步骤,依赖于时间的牛顿方程的积分,以及总能量守恒。
用统计的说法,这些模拟是在微正则体系中进行的或者NVE体系:
粒
子的数目,体积以及能量是恒定的物理量。
物理量A的微观平均是通过对轨迹的时间平均获得的:
I
〈4〉nve=-r-力4(『(±
))
在这里我们用-(t)指示系统的相空间的坐标(3N个位置和3N个速度),
也存在一些其他重要的体系,这里我们只做简单的介绍。
一个基本的想法就是对其他方程进行积分而不是对牛顿方程进行积分,通过这种方法取样是在其他的统计体系中进行。
然而在新的体系中获得的物理量的平均与以前介绍的时间平均一样。
Andersen发展了一种在等焓等压体系(NPH)进行模拟的方法。
在这里,添加一个自由度代表引进的盒子的体积,所有的粒子相对盒子都给定坐标。
盒子的体积V在模拟过程中是可变的。
焓H=E+PV是恒定的,P是外部压力。
Parrinello和Rahman发展了一种不同的方法.盒子的形状随体积一起变化,这种方法新增了9个自由度而不是一个:
三个向量新生成MD盒子。
它们每一个都是动态可变的,根据运动方程(来自一个合适的拉格朗日因子)演化。
这种策略可以用来研究结构的相转变,(作为P的函数),例如系统放弃一个特定的晶体结构而去适合一
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 分子 动力学 入门 第三