使用ORCA做从头算动力学(AIMD)的简单事例
Sobereva@上海科音
Firstrelease:2020-Dec-13Lastupdate:2021-Jul-7
1序言
基于量子物理方式的动力学通常称为从头算动力学（Abinitiomoleculardynamics,AIMD），相比于基于通常的精典力场的动力学，其关键优势在于精度高、普适性强、能够描述物理反应，代价是历时相差N个数目级。ORCA量子物理程序有不错的做BOMD方式的从头算动力学的功能，使用很便捷，但是本身ORCA做DFT的效率又高，是做孤立体系AIMD的首选程序之一。其实有些特点不支持，例如无法像Gaussian的BOMD那样直接做准精典动力学，不能按照原子距离等标准判定哪些时侯手动结束任务等等，但都不是大问题。对于跑跑普通的AIMD来说，笔者觉得ORCA显著比Gaussian的ADMP或BOMD更好用（Gaussian的AIMD输入文件较为具象，指南相应部份写得很烂，但是连个像样的热浴都没有），但是速率也显著更快。笔者近日发表的文章中对18碳环（cyclo[18]carbon）单体和二聚体做的分子动力学就是用ORCA跑的，分别见Chem.Asian.J.(2020)DOI:10.1002/asia.202001228和《全面探究18碳环奇特的分子间互相作用与pi-pi堆积特点》（）中对Carbon,171,514(2021)一文的介绍，可以作为ORCA跑AIMD的范例。
笔者在上海科音（）的中级量子物理培训班中会用多达两百多页的ppt专门深入详尽讲AIMD的模拟，其中也包括ORCA做AIMD的各类细节、大量方法以及众多实例。本文只是举一个简单的事例，帮助读者基本了解ORCA怎么做AIMD。弄明白这个反例，结合指南举一反三，跑其它体系也就没哪些难的了。假如对ORCA一无所知的话，点击左侧的“ORCA”分类，笔者之前写过不少相关文章。ORCA只适宜跑孤立体系的从头算动力学，若果是做周期性体系的从头算动力学，CP2K是最佳的选择，在笔者讲授的上海科音CP2K第一性原理估算培训班（）中有极为全面、系统、详细的讲解并给了大量反例。
本文内容适用于Multiwfn最新版本、VMD1.9.3、ORCA4.2.x及之后版本的情况。
2021-Jul-7针对ORCA5.0的补充说明：对于2021-Jul-7及之后更新的Multiwfn，步入本文所述的Multiwfn形成输入文件的界面后，可以通过选项-11选择适宜的ORCA版本，默认为ORCA5.0及之后版本，而下文内容对应的是4.2.x版的情况。对于>=5.0版，Multiwfn手动设的热浴是CSVR，比之前版本惟一支持的Berendsen热浴更好，并且同样普适。并且在run前面多加了CenterCOM，这是从5.0版本开始支持的清除整体刚体运动的选项，而下文里提及的constraintaddcenter这一行就没有了。
2实例：[Al(H2O)6]3+与NH3之间的质子转移
ORCA的MD模块的开发者网站上有个真空中[Al(H2O)6]3+与NH3之间的质子转移的动漫，如下所示
#a:7:e:6:9:4:3:5:7:3:8:3:5:7:b:e:0:f:7:c:c:6:f:f:7:4:e:b:0:0:0:6#
可见NH3向带正电的[Al(H2O)6]3+渐渐紧靠，水上的一个质子转移到了气体分子上，之后因为静电互斥，NH4+就渐渐远离[Al(H2O)5(OH)]2+了。这儿我们企图再现这个过程。下边提及的文件都可以在中获得。
我们首先获得[Al(H2O)6]3+的基本合理的结构。其实用ORCA优化也可以，这儿笔者习惯性地用Gaussian来优化。在GaussView里搭建Al(H2O)6，保存为Al(H2O)6_optfreq.gjf，将关键词改为#B3LYP/6-31G*optfreq，将电荷改为3，之后用Gaussian运行之，就得到了优化后的[Al(H2O)6]3+结构。再用GaussView打开输出文件Al(H2O)6_optfreq.out，把一个气体分子画在一个水的门口，如下所示，之后保存为complex.gjf。
#0:a:e:f:d:e:0:b:8:1:3:2:b:f:c:6:3:0:4:0:a:f:9:f:8:6:7:9:8:8:0:9#
Multiwfn程序（）有很便利的生成ORCA常见任务的输入文件的功能，见《详谈Multiwfn形成ORCA量子物理程序的输入文件的功能》（），这儿我们用Multiwfn生成ORCA的AIMD任务的标准输入文件。
启动Multiwfn，载入complex.gjf，之后输入
oi//生成ORCA输入文件
0//选择任务类型
6//分子动力学
1//估算级别用B97-3c
在当前目录下就得到了ORCA的AIMD任务的标准输入文件complex.inp。
在complex.inp上面将相应几行改成下边这样：
%maxcore10000
%palnprocs10end
timestep1.0_fs
initvel298.15_K
*xyz31
现今的complex.inp的完整内容如下。以#为开头的行代表旁边的设置被注释掉了，不会生效，想启用可以去除#。此文件里也有大量Multiwfn手动添加的注释，只要一看注释马上就明白相应的行是哪些含意，巨贴心，都省得查指南了。
!B97-3cnoautostartminiprintnopop
%maxcore10000
%palnprocs10end
%md
#restartifexists#ContinueMDbyreading[basename].mdrestartifitexists.Inthiscase"initvel"shouldbecommented
#minimize#DominimizationpriortoMDsimulation
timestep1.0_fs#ThisstepsizeissafeatseveralhundredsofKelvin
initvel298.15_Kno_overwrite#Assignvelocityaccordingtotemperatureforatomswhosevelocitiesarenotavailable
thermostatberendsen298.15_Ktimecon30.0_fs#Targettemperatureandcouplingtimeconstant
dumppositionstride1formatxyzfilename"pos.xyz"#Dumppositionevery"stride"steps
#dumpforcestride1formatxyzfilename"force.xyz"#Dumpforceevery"stride"steps
#dumpvelocitystride1formatxyzfilename"vel.xyz"#Dumpvelocityevery"stride"steps
#dumpgbwstride20filename"wfn"#Dumpwavefunctionto"wfn[step].gbw"filesevery"stride"steps
constraintaddcenter0..22#Fixcenterofmassatinitialposition
run2000#NumberofMDsteps
end
*xyz31
[座标部份]
下边简单说一下complex.inp里这种设置的涵义、为什么那么设。
•B97-3c是一个又实惠又不错的估算级别，在ORCA里还手动会启用RIJ加速，速率很快，因而很适宜跑AIMD，描述当前体系没有问题。B97-3c的介绍见《详谈Multiwfn形成ORCA量子物理程序的输入文件的功能》（）的相应部份。但B97-3c也不是哪些时侯都能用，例如18碳环用纯泛函描述皆失败，见笔者在Carbon,165,468(2020)里的讨论，其实就不能用这个了，笔者跑涉及18碳环的AIMD的时侯都用的是ωB97X-D3。
•%palnprocs10end代表用10核。笔者当前任务实际上是在双路E5-2696v3共36个数学核心的机子上跑的，但却故意用了10核。这是由于按照笔者曾经的测试，发觉ORCA做DFT的AIMD的并行效率不理想，尤其是对于小体系，用核数太多甚至反而速率更慢（你们可以对实际情况短暂跑例如5步实测一下设多少核速率最快）。通常来说就设10核就行了，当机子里有显著更多核的时侯，可以跑多个AIMD任务来充分借助估算资源，但应该对CPU内核进行绑定，否则AIMD估算速率可能明显增加，见《通过设置CPU内核绑定减少ORCA同时做空任务的历时》（）。
•%maxcore10000代表每位ORCA进程最多用10000MB。虽然完全没必要如此大，普通泛函下的AIMD不怎样耗显存，给1000都绝对够了。
•restartifexists这句被注释掉了。假如你的AIMD想续跑，且当前目录下有之前跑下来的与当前任务同名的后缀为.mdrestart的文件，可以取消注释，任务都会延续之前AIMD最后的状态续跑，新轨迹会在原有轨迹文件前面续写。
•minimize这句被注释掉了。假如取消注释的话，动力学开始前会手动在当前级别下做几何优化。
•timestep1.0_fs代表动力学步长用1fs。对于此例常温下的模拟，1fs步长是合适的，改大的话可能导致动力学不稳定，而改小的话跑同样的时间宽度须要更多步数造成须要更多历时。假如追求绝对稳当，或则是在显著更高湿度下模拟，可以用0.5fs。
•initvel298.15_K代表按照298.15K气温通过Maxwell分布随机初始化原子速率。no_overwrite代表假如之前早已有初速率信息了（例如可以是续跑时从之前的mdrestart文件里承继来的），就不再形成新的初速率。
•thermostatberendsen298.15_Ktimecon30.0_fs代表用Berendsen热浴温控在298.15K（此例笔者用的ORCA4.2.1只能用这个热浴，此外从ORCA5.0开始可以用更好的velocityrescale热浴），时间常数为30fs，一般这个时间常数是合适的。
•dumppositionstride1formatxyzfilename"pos.xyz"代表每一步都往当前目录下的pos.xyz文件里写入当前的座标，因而pos.xyz是多帧xyz格式的轨迹文件。此格式的介绍见《谈谈记录物理体系结构的xyz文件》（）。
•dumpforce和dumpvelocity开头的行都被注释掉了，这两行用于把模拟过程中的原子受力和原子速率分别保存到相应xyz文件里。dumpgbw开头的行也被注释掉了，它可以在MD过程中每隔指定步数保存一次gbw文件，然后可以通过写脚本调用Multiwfn对它们进行批量剖析，因而考察动力学过程中电子结构变化，得到丰富的有物理意义的信息（如动力学中的电荷转移情况、成键变化情况、电子定域性变化等等，在上海科音中级量子物理培训班里笔者会给出好多这些反例）
•constraintaddcenter0..22代表将当前整个体系（0号到22号原子）的刚体约束在初始位置。之所以如此做，是由于虽然ORCA在形成初速率时早已除去了整体平移的速率份量，但实际模拟过程中因为数值问题，一直常常会听到有显著的整体甩尾的现象，因此有碍观察（在VMD里观看轨迹时还得做一下align能够消去）。因而直接把刚体位置约束住就没有这个问题了。
•run2000代表总共跑2000步，当前步长是1fs，即最多跑2ps。当前这个模拟的目的是观察到质子转移，跑多长时间合适并不好预估，所以可以一次先跑2000步瞧瞧，若不够待会儿再续跑。值得一提的是，起码对于笔者现今用的ORCA4.2.1，我发觉假如一次跑的步数好多，达到2000步左右以后，然后每一步的历时会有渐渐的上升趋势，缘由不明。为此我建议每次跑最多不宜超过3000步，假如须要跑更长的话，最好分多段来跑。
如今用ORCA运行complex.inp，模拟过程中可以看见每一步的实时情况：
#9:0:6:6:9:1:a:2:0:3:a:1:5:f:3:5:c:e:3:f:0:7:5:e:4:c:3:e:8:f:4:1#
Step是当前的步数，Time是当前的时间，t_SCF和t_Grad分别是算这一步的SCF过程和估算受力的历时，两者加和就是这一步的总历时。可见每一步历时约6~7秒，除以要跑的步数就可以恐怕跑完整个任务的历时。旁边还可以看见每一步的体温、动能、势能等信息。因为当前体系原子数极少，气温相对于热浴的参考气温波动大是很正常的事。并且因为用了热浴，所以总能量E_tot也有显著波动。
VMD是观看动力学轨迹最好的程序，可以在免费下载。建议在模拟过程中，隔一阵子就用VMD把pos.xyz载入进去，瞧瞧当前的动态行为怎么、跑成哪些结构了。挪到289步的时侯，笔者看了一眼pos.xyz，发觉质子除了早已完全转移，但是NH4+都早已跑走一定距离了，所以就没必要继续跑了，故把ORCA直接杀掉了。
模拟过程中ORCA还形成其它一些文件。complex.md.log是ORCA的MD模块输出的信息，相当于整个输出文件中的子集。complex.mdrestart是拿来续跑的文件，每一步就会往上面写入当前步的时间、坐标、速度、受力等信息。complex-md-ener.csv是把每一步的时间、耗时、温度、能量等信息以csv格式保存的文件。还有其它一些零碎的文件，一般不是通常用户关心的，这儿就不说了，都可以放心删除，留着也没用。
如今我们来看模拟轨迹。建议你们按照《VMD初始化文件(vmd.rc)我的推荐设置》（）里说的更改VMD的初始化文件，因而添加自定义命令bt，这样在播放轨迹的时侯对每一帧会手动重新判定成键关系。
启动VMD，载入pos.xyz（在本文文件包里已提供，共290帧）。之后在文本窗口输入bt，按回车，致使每帧都更新联接关系。之后再输入bw，按回车促使背景变为黄色。选Graphics-Representation，将DrawingMethod改为CPK。之后点击VMD界面右下角的三角播放动漫，会看见如下结果。假如想在VMD图形窗口中显示出帧号或时间，看《使VMD播放轨迹时同步显示帧号》（）。
#a:a:1:f:d:2:7:0:6:2:a:4:1:3:6:e:8:a:6:e:8:e:c:0:c:5:e:9:c:0:8:e#
可见，我们跑下来的动力学现象和本文一开始的那种官方动漫几乎完全一致。都是两个反应物先接近，之后产生复合状态时质子在两者之间回落，最后NH4+飞走。
假如想把质子转移情况通过距离随时间的变化曲线形式呈送给读者，可以在点击VMD的三维图形窗口后，按按键上的2键（以后按r可以恢复为默认的旋转模式），之后点击两个原子正中心，两者之间都会降低Bondlabel（默认是以红色字显示距离，在黑背景下才看得清楚），这儿笔者把N和转移过去的质子之间降低了Bondlabel。之后步入Graphics-Labels，切换到Bonds，选Graph，点击Showpreview复选框，之后点击1:H1:N这项，都会听到距离变化早已显示在预览窗口了，如下所示，其中蓝色和白色条纹标记的分别是最小距离和最大距离位置和数值。
#b:e:8:8:7:c:d:d:7:a:c:0:8:5:4:4:4:2:0:3:c:f:8:e:8:e:9:3:e:d:4:d#
假如点击Graph按键，都会把曲线显示在大窗口中，如下所示。可见，在大概60帧，也即60fs左右，N-H键即使是基本产生了，然后N-H键不断震动。
#b:7:d:7:f:e:2:4:3:b:3:e:b:3:f:8:1:d:7:d:3:9:1:e:b:8:9:5:6:1:1:9#
以类似方法，我们可以标记Al与N的距离，随时间变化如下所示。可见Al与N先接近，质子转移完毕后，两者就渐渐远离了。
#4:c:b:f:1:9:3:4:3:8:b:b:6:5:5:6:4:f:a:a:3:5:d:7:c:c:0:a:c:7:a:3#
在Labels窗口里还可以点击save，把距离变化数据保存到文本文件里，然后可以导出Origin等程序里勾画曲线。
用VMD还可以检测角度、二面角的变化，分别是在图形窗口里按3之后点三个原子、按4之后点四个原子进行标记，然后在Graphics-Labels里观看。
3总结&其它
通过前面的事例，可以看见ORCA做AIMD是相当容易的，只要把Multiwfn支持的富含结构信息的文件（如pdb/gjf/xyz/mol/mol2/fch等等，见Multiwfn指南2.5节）载入Multiwfn，敲几下按键形成标准AIMD任务的输入文件，之后按照实际情况稍稍改几个设定就可以跑了。
以几十核的通常双路服务器的运算能力，ORCA里用B97-3c跑几十原子有机体系的几十ps的动力学不是非常困难的事。不过，能跑的时间尺度仍远远比不上xtb跑半经验层面DFT的GFN-xTB方式的动力学，xtb跑位热学的肤浅介绍和反例看此文的相应部份：《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构型搜索》（）。因而，拿ORCA跑DFT的动力学之前，先拿xtb初步跑跑，找找觉得，大体摸索出自己期望的现象能出现的条件（如体温、初始结构、反应物相对位置和碰撞方法等），之后再用DFT跑一般是比较好的做法，免得做高昂的DFT的MD试来试去把时间都耽搁了。
本文只涉及了VMD一丁点皮毛，VMD对于做动力学的人是必须玩得十分溜的。笔者在上海科音分子动力学与GROMACS培训班（）里对VMD有特别深入全面的介绍，包括tcl脚本的编撰。借助VMD的tcl脚本可以对轨迹做千变万化的剖析，有些剖析例如刚体距离变化、平面间倾角变化、某几何变量分布统计、不同结构出现百分比等，是必须靠写脚本能够实现的。
光是剖析剖析动力学过程的能量、结构变化是很浅显的。借助Multiwfn，可以对ORCA跑的动力学的全过程的电子结构做深入透彻的剖析，因而考察物理键、弱互相作用、电荷分布等等在动力学过程中的变化，由此才能从提供深入的视角，使研究文章信息更丰富、明显更上档次。特别建议详尽阅读《详谈Multiwfn的命令行形式运行和批量运行的方式》（），上面第4节专门讲了如何做这样的剖析，你会发觉非常容易也非常有价值。
假如Multiwfn创建ORCA做动力学输入文件的功能对你的实际研究形成了帮助，希望在写文章的时侯提到例如Theinputfileofab-initiomoleculardynamicswaspreparedwiththehelpofMultiwfnprogram并引用Multiwfn原文，这是对Multiwfn此功能开发最好的支持。
