计算化学公社

标题: 使用ORCA做从头算动力学(AIMD)的简单例子 [打印本页]

作者
Author: sobereva 时间: 2020-12-13 06:04
标题: 使用ORCA做从头算动力学(AIMD)的简单例子

使用ORCA做从头算动力学(AIMD)的简单例子

A simple example of performing ab initio dynamics (AIMD) using ORCA

Sobereva@北京科音

First release: 2020-Dec-13 Last update: 2021-Jul-7

1 前言

基于量子化学方法的动力学一般称为从头算动力学（Ab initio molecular dynamics, AIMD），相比于基于一般的经典力场的动力学，其关键优势在于精度高、普适性强、能够描述化学反应，代价是耗时相差N个数量级。ORCA量子化学程序有不错的做BOMD形式的从头算动力学的功能，使用很方便，而且本身ORCA做DFT的效率又高，是做孤立体系AIMD的首选程序之一。虽然有些特性不支持，比如没法像Gaussian的BOMD那样直接做准经典动力学，不能根据原子距离等标准判断什么时候自动结束任务等等，但都不是大问题。对于跑跑普通的AIMD来说，笔者感觉ORCA明显比Gaussian的ADMP或BOMD更好用（Gaussian的AIMD输入文件较为抽象，手册相应部分写得很烂，而且连个像样的热浴都没有），而且速度也明显更快。

笔者发表的18碳环（cyclo[18]carbon）的研究论文中，其单体和二聚体做的分子动力学就是用ORCA跑的，分别见《揭示各种新奇的碳环体系的振动特征》（http://sobereva.com/578）对中Chem. Asian J., 16, 56 (2021)和《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（http://sobereva.com/572）中对Carbon, 171, 514 (2021)一文的介绍。在《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（http://sobereva.com/696）中介绍的笔者的Inorg. Chem., 62, 19986 (2023)文章中还用ORCA跑了B6C6N6的高温动力学以证明其稳定性。在《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（http://sobereva.com/696）中介绍的笔者的ChemPhysChem, 25, e202400377 (2024)文章中还用ORCA跑了18氮环的动力学以考察其热分解行为。这些文章可以作为ORCA跑AIMD的典型范例，很推荐阅读和引用。

笔者在北京科音高级量子化学培训班（http://www.keinsci.com/workshop/KAQC_content.html）中会用多达两百多页的ppt专门深入详细讲AIMD的模拟，其中也包括ORCA做AIMD的各种细节、大量技巧以及诸多实例，并且此培训中还会系统讲授ORCA程序的使用，因此是使用ORCA专门做AIMD的读者一定不可错过的培训。而本文只是举一个简单的例子，帮助读者快速了解ORCA如何做最简单的AIMD。注意ORCA只适合跑孤立体系的从头算动力学，如果是做周期性体系的从头算动力学，CP2K是最佳的选择，在笔者讲授的北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）中有极为全面、系统、详细的讲解并给了大量例子。

本文内容适用于Multiwfn最新版本、VMD 1.9.3、ORCA 4.2.x及以后版本的情况。

2021-Jul-7 针对ORCA 5.0的补充说明：对于2021-Jul-7及以后更新的Multiwfn，进入本文所述的Multiwfn产生输入文件的界面后，可以通过选项-11选择适合的ORCA版本，默认为ORCA 5.0及以后版本，而下文内容对应的是4.2.x版的情况。对于>=5.0版，Multiwfn自动设的热浴是CSVR，比之前版本唯一支持的Berendsen热浴更好，而且同样普适。而且在run后面多加了CenterCOM，这是从5.0版本开始支持的消除整体质心运动的选项，而下文里提到的constraint add center这一行就没有了。

2 实例：[Al(H2O)6]3+与NH3之间的质子转移

ORCA的MD模块的开发者网站上有个真空中[Al(H2O)6]3+与NH3之间的质子转移的动画，如下所示

可见NH3向带正电的[Al(H2O)6]3+逐渐靠近，水上的一个质子转移到了氨气分子上，然后由于静电互斥，NH4+就逐渐远离[Al(H2O)5(OH)]2+了。这里我们试图重现这个过程。下面提到的文件都可以在http://sobereva.com/attach/576/file.zip中获得。

我们首先获得[Al(H2O)6]3+的基本合理的结构。当然用ORCA优化也可以，这里笔者习惯性地用Gaussian来优化。在GaussView里搭建Al(H2O)6，保存为Al(H2O)6_optfreq.gjf，将关键词改为# B3LYP/6-31G* opt freq，将电荷改为3，然后用Gaussian运行之，就得到了优化后的[Al(H2O)6]3+结构。再用GaussView打开输出文件Al(H2O)6_optfreq.out，把一个氨气分子画在一个水的旁边，如下所示，然后保存为complex.gjf。

(, 下载次数 Times of downloads: 131)

Multiwfn程序（http://sobereva.com/multiwfn）有很便利的生成ORCA常见任务的输入文件的功能，见《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》（http://sobereva.com/490），这里我们用Multiwfn生成ORCA的AIMD任务的标准输入文件。

启动Multiwfn，载入complex.gjf，然后输入
oi  // 生成ORCA输入文件
0  // 选择任务类型
6  // 分子动力学
1  // 计算级别用B97-3c
在当前目录下就得到了ORCA的AIMD任务的标准输入文件complex.inp。

在complex.inp里面将相应几行改成下面这样：
%maxcore  10000
%pal nprocs 10 end
timestep 1.0_fs
initvel 298.15_K
* xyz 3 1

现在的complex.inp的完整内容如下。以#为开头的行代表后面的设置被注释掉了，不会生效，想启用可以去掉#。此文件里也有大量Multiwfn自动添加的注释，只要一看注释马上就明白相应的行是什么含义，巨贴心，都省得查手册了。

! B97-3c noautostart miniprint nopop
%maxcore  10000
%pal nprocs 10 end
%md
#restart ifexists  # Continue MD by reading [basename].mdrestart if it exists. In this case "initvel" should be commented
#minimize  # Do minimization prior to MD simulation
timestep 1.0_fs  # This stepsize is safe at several hundreds of Kelvin
initvel 298.15_K no_overwrite # Assign velocity according to temperature for atoms whose velocities are not available
thermostat berendsen 298.15_K timecon 30.0_fs  # Target temperature and coupling time constant
dump position stride 1 format xyz filename "pos.xyz"  # Dump position every "stride" steps
#dump force stride 1 format xyz filename "force.xyz"  # Dump force every "stride" steps
#dump velocity stride 1 format xyz filename "vel.xyz"  # Dump velocity every "stride" steps
#dump gbw stride 20 filename "wfn"  # Dump wavefunction to "wfn[step].gbw" files every "stride" steps
constraint add center 0..22  # Fix center of mass at initial position
run 2000  # Number of MD steps
end
* xyz 3 1
[坐标部分]
*

下面简单说一下complex.inp里这些设置的含义、为什么这么设。
• B97-3c是一个又便宜又不错的计算级别，在ORCA里还自动会启用RIJ加速，速度很快，因此很适合跑AIMD，描述当前体系没有问题。B97-3c的介绍见《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》（http://sobereva.com/490）的相应部分。但B97-3c也不是什么时候都能用，比如18碳环用纯泛函描述皆失败，见笔者在Carbon, 165, 468 (2020) 里的讨论，显然就不能用这个了，笔者跑涉及18碳环的AIMD的时候都用的是ωB97X-D3。
• %pal nprocs 10 end代表用10核。笔者当前任务实际上是在双路E5-2696v3共36个物理核心的机子上跑的，但却故意用了10核。这是因为根据笔者以前的测试，发现ORCA做DFT的AIMD的并行效率不理想，尤其是对于小体系，用核数太多甚至反倒速度更慢（大家可以对实际情况短暂跑比如5步实测一下设多少核速度最快）。一般来说就设10核就行了，当机子里有明显更多核的时候，可以跑多个AIMD任务来充分利用计算资源，但应当对CPU内核进行绑定，否则AIMD计算速度可能显著降低，见《通过设置CPU内核绑定降低ORCA同时做多任务的耗时》（http://sobereva.com/553）。
• %maxcore  10000代表每个ORCA进程最多用10000MB。其实完全没必要这么大，普通泛函下的AIMD不怎么耗内存，给1000都绝对够了。
• restart ifexists这句被注释掉了。如果你的AIMD想续跑，且当前目录下有之前跑出来的与当前任务同名的后缀为.mdrestart的文件，可以取消注释，任务就会延续之前AIMD最后的状态续跑，新轨迹会在原有轨迹文件后面续写。
• minimize这句被注释掉了。如果取消注释的话，动力学开始前会自动在当前级别下做几何优化。
• timestep 1.0_fs代表动力学步长用1 fs。对于此例常温下的模拟，1 fs步长是合适的，改大的话可能造成动力学不稳定，而改小的话跑同样的时间长度需要更多步数导致需要更多耗时。如果追求绝对稳妥，或者是在明显更高温度下模拟，可以用0.5 fs。
• initvel 298.15_K代表根据298.15 K温度通过Maxwell分布随机初始化原子速度。no_overwrite代表如果之前已经有初速度信息了（比如可以是续跑时从之前的mdrestart文件里继承来的），就不再产生新的初速度。
• thermostat berendsen 298.15_K timecon 30.0_fs代表用Berendsen热浴控温在298.15 K（此例笔者用的ORCA 4.2.1只能用这个热浴，据悉从ORCA 5.0开始可以用更好的velocity rescale热浴），时间常数为30 fs，通常这个时间常数是合适的。
• dump position stride 1 format xyz filename "pos.xyz"代表每一步都往当前目录下的pos.xyz文件里写入当前的坐标，因此pos.xyz是多帧xyz格式的轨迹文件。此格式的介绍见《谈谈记录化学体系结构的xyz文件》（http://sobereva.com/477）。
• dump force和dump velocity开头的行都被注释掉了，这两行用于把模拟过程中的原子受力和原子速度分别保存到相应xyz文件里。dump gbw开头的行也被注释掉了，它可以在MD过程中每隔指定步数保存一次gbw文件，之后可以通过写脚本调用Multiwfn对它们进行批量分析，从而考察动力学过程中电子结构变化，得到丰富的有化学意义的信息（如动力学中的电荷转移情况、成键变化情况、电子定域性变化等等，在北京科音高级量子化学培训班里笔者会给出很多这种例子）
• constraint add center 0..22代表将当前整个体系（0号到22号原子）的质心约束在初始位置。之所以这么做，是因为尽管ORCA在产生初速度时已经去掉了整体平移的速度分量，但实际模拟过程中由于数值问题，仍然往往会看到有明显的整体漂移的现象，因而有碍观察（在VMD里观看轨迹时还得做一下align才能消掉）。因此直接把质心位置约束住就没有这个问题了。
• run 2000代表总共跑2000步，当前步长是1 fs，即最多跑2 ps。当前这个模拟的目的是观察到质子转移，跑多长时间合适并不好预估，所以可以一次先跑2000步看看，若不够到时候再续跑。值得一提的是，至少对于笔者现在用的ORCA 4.2.1，我发现如果一次跑的步数很多，达到2000步左右之后，之后每一步的耗时会有逐渐的上升趋势，原因不明。因此我建议每次跑最多不宜超过3000步，如果需要跑更长的话，最好分多段来跑。

现在用ORCA运行complex.inp，模拟过程中可以看到每一步的实时情况：

(, 下载次数 Times of downloads: 137)

Step是当前的步数，Time是当前的时间，t_SCF和t_Grad分别是算这一步的SCF过程和计算受力的耗时，二者加和就是这一步的总耗时。可见每一步耗时约6~7秒，乘以要跑的步数就可以估计跑完整个任务的耗时。后面还可以看到每一步的温度、动能、势能等信息。由于当前体系原子数很少，温度相对于热浴的参考温度波动大是很正常的事。而且由于用了热浴，所以总能量E_tot也有明显波动。

VMD是观看动力学轨迹最好的程序，可以在http://www.ks.uiuc.edu/Research/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)我的推荐设置》（http://sobereva.com/545）里说的修改VMD的初始化文件，从而添加自定义命令bt，这样在播放轨迹的时候对每一帧会自动重新判断成键关系。

启动VMD，载入pos.xyz（在本文文件包里已提供，共290帧）。然后在文本窗口输入bt，按回车，使得每帧都更新连接关系。然后再输入bw，按回车使得背景变为白色。选Graphics - Representation，将Drawing Method改为CPK。然后点击VMD界面右下角的三角播放动画，会看到如下结果。如果想在VMD图形窗口中显示出帧号或时间，看《使VMD播放轨迹时同步显示帧号》（http://sobereva.com/13）。

可见，我们跑出来的动力学现象和本文一开始的那个官方动画几乎完全一致。都是两个反应物先接近，然后形成复合状态时质子在二者之间震荡，最后NH4+跑掉。

如果想把质子转移情况通过距离随时间的变化曲线方式呈献给读者，可以在点击VMD的三维图形窗口后，按键盘上的2键（之后按r可以恢复为默认的旋转模式），然后点击两个原子正中心，二者之间就会增加Bond label（默认是以白色字显示距离，在黑背景下才看得清楚），这里笔者把N和转移过去的质子之间增加了Bond label。然后进入Graphics - Labels，切换到Bonds，选Graph，点击Show preview复选框，然后点击1:H 1:N这项，就会看到距离变化已经显示在预览窗口了，如下所示，其中红色和蓝色竖条标记的分别是最小距离和最大距离位置和数值。

(, 下载次数 Times of downloads: 160)

如果点击Graph按钮，就会把曲线显示在大窗口中，如下所示。可见，在大约60帧，也即60 fs左右，N-H键就算是基本形成了，之后N-H键不断振动。

(, 下载次数 Times of downloads: 149)

以类似方式，我们可以标记Al与N的距离，随时间变化如下所示。可见Al与N先接近，质子转移完毕后，二者就逐渐远离了。

(, 下载次数 Times of downloads: 134)

在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)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）。因此，拿ORCA跑DFT的动力学之前，先拿xtb初步跑跑，找找感觉，大体摸索出自己期望的现象能出现的条件（如温度、初始结构、反应物相对位置和碰撞方式等），然后再用DFT跑通常是比较好的做法，免得做昂贵的DFT的MD试来试去把时间都耽误了。

本文只涉及了VMD一丁点皮毛，VMD对于做动力学的人是必须玩得非常溜的。笔者在北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/workshop/KGMX_content.html）里对VMD有非常深入全面的介绍，包括tcl脚本的编写。利用VMD的tcl脚本可以对轨迹做千变万化的分析，有些分析诸如质心距离变化、平面间夹角变化、某几何变量分布统计、不同结构出现比率等，是必须靠写脚本才能实现的。

光是分析分析动力学过程的能量、结构变化是很肤浅的。利用Multiwfn，可以对ORCA跑的动力学的全过程的电子结构做深入透彻的分析，从而考察化学键、弱相互作用、电荷分布等等在动力学过程中的变化，由此能够从提供深入的视角，使研究文章信息更丰富、明显更上档次。非常建议详细阅读《详谈Multiwfn的命令行方式运行和批量运行的方法》（http://sobereva.com/612），里面第4节专门讲了怎么做这样的分析，你会发现特别容易也特别有价值。

如果Multiwfn创建ORCA做动力学输入文件的功能对你的实际研究产生了帮助，希望在写文章的时候提及诸如The input file of ab-initio molecular dynamics was prepared with the help of Multiwfn program并引用Multiwfn原文，这是对Multiwfn此功能开发最好的支持。

作者
Author: zsu007 时间: 2020-12-13 09:32
谢谢社长的分享！

作者
Author: funok 时间: 2020-12-13 12:25
谢谢社长的分享

作者
Author: fineren 时间: 2020-12-13 18:58
感谢分享，可以模拟燃烧反应吗？另外高级班预计啥时候开？

作者
Author: Lacrimosa 时间: 2020-12-13 20:23
请教一下卢老师，可以理解为AIMD和GFN-xTB方法的动力学基本是相似的，只是在精度上有差别么？另外xtb跑动力学的时候能否也输出模拟过程中的电子结构信息呢？之前简单尝试过一下，发现只输出了每一步的坐标。

作者
Author: sobereva 时间: 2020-12-13 23:00

fineren 发表于 2020-12-13 18:58
感谢分享，可以模拟燃烧反应吗？另外高级班预计啥时候开？

可以

北京科音办的培训班FAQ
http://bbs.keinsci.com/thread-5098-1-1.html

作者
Author: sobereva 时间: 2020-12-13 23:02

Lacrimosa 发表于 2020-12-13 20:23
请教一下卢老师，可以理解为AIMD和GFN-xTB方法的动力学基本是相似的，只是在精度上有差别么？另外xtb跑动力 ...

ORCA和xtb程序做的动力学都是BOMD的形式，只是计算受力的理论方法不同，ORCA做DFT的比GFN-xTB精度更好。

xtb没这个功能。想实现的话只能写个脚本，对跑完了的轨迹每隔一定帧数算一次单点产生molden文件。

作者
Author: gog 时间: 2020-12-13 23:17
看到的算例，高斯方法，都是算弱相互作用吧。无机材料，MOF，COF材料的物理性质，能用orca计算么？

作者
Author: sobereva 时间: 2020-12-14 01:25

gog 发表于 2020-12-13 23:17
看到的算例，高斯方法，都是算弱相互作用吧。无机材料，MOF，COF材料的物理性质，能用orca计算么？

Gaussian是程序，不是方法
Gaussian从来不是专拿来算弱相互作用的程序，虽然算弱相互作用是Gaussian用户常干的事之一
Gaussian和ORCA的主要应用范畴没有显著的差异

Gaussian常被拿来算什么问题，看看本版块Gaussian分类里的帖子便知。

另外建议看
谈谈学量子化学如何正确地入门
http://sobereva.com/355（http://bbs.keinsci.com/thread-4447-1-1.html）

作者
Author: alonewolfyang 时间: 2020-12-27 19:15
老师，您这里的例子在ORCA 4.0.0版不能正常计算

作者
Author: sobereva 时间: 2020-12-28 07:02

alonewolfyang 发表于 2020-12-27 19:15
老师，您这里的例子在ORCA 4.0.0版不能正常计算

当然不能。为什么要用那么老的

作者
Author: yinhang 时间: 2020-12-28 15:30
请教老师，ORCA的动力学计算能跨过势垒么，特别是激发态的动力学，谢谢您~

作者
Author: sobereva 时间: 2020-12-29 08:06

yinhang 发表于 2020-12-28 15:30
请教老师，ORCA的动力学计算能跨过势垒么，特别是激发态的动力学，谢谢您~

能

作者
Author: 喝了假咖啡 时间: 2020-12-29 09:16
请问卢老师，nitvel 298.15_K 以及 thermostat berendsen 298.15_K timecon 30.0_fs，设置参数时这里面的两个温度有什么要注意吗，比如要研究某分子在500 K的热稳定性，我觉得应该写thermostat berendsen 500_K timecon 30.0_fs，而另一个温度参数该如何设置？

作者
Author: sobereva 时间: 2020-12-29 11:10

喝了假咖啡发表于 2020-12-29 09:16
请问卢老师，nitvel 298.15_K 以及 thermostat berendsen 298.15_K timecon 30.0_fs，设置参数时这里 ...

最终要算什么温度，热浴参考温度就设多少。

产生初速度对应的温度不要设得太高，对常温产生初速度还好，如果对高温直接产生初速度的话由于某些原子运动速度太快，容易造成动力学不稳定、出现异常结构、SCF难收敛等问题

作者
Author: TDHFjiang 时间: 2020-12-31 12:52
l老师，您好，保存每一步的gbw 文件，请问如何用Multiwfn进行电荷转移分析?

作者
Author: 土拔鼠 时间: 2020-12-31 16:48
卢老师请问如果在ReaxFF下模拟了分子的热解离但是感觉和实验符合不太好我是否可以尝试AIMD去试试
和ReaxFF相比 AIMD的结果会有定性的变化吗？
问题比较蠢忘社长见谅

作者
Author: sobereva 时间: 2021-1-2 00:08

TDHFjiang 发表于 2020-12-31 12:52
l老师，您好，保存每一步的gbw 文件，请问如何用Multiwfn进行电荷转移分析?

写脚本批量计算批量提取。例如
(, 下载次数 Times of downloads: 105)
(, 下载次数 Times of downloads: 83)
(, 下载次数 Times of downloads: 108)

作者
Author: sobereva 时间: 2021-1-2 00:08

土拔鼠发表于 2020-12-31 16:48
卢老师请问如果在ReaxFF下模拟了分子的热解离但是感觉和实验符合不太好我是否可以尝试AIMD去试试
和Rea ...

可
很有可能会

作者
Author: alonewolfyang 时间: 2021-1-5 11:12

sobereva 发表于 2020-12-29 11:10
最终要算什么温度，热浴参考温度就设多少。

产生初速度对应的温度不要设得太高，对常温产生初速度还好 ...

老师，我将这两种温度设得不一样，但MD也跑完了，这种情况没事吧

作者
Author: sobereva 时间: 2021-1-5 13:12

alonewolfyang 发表于 2021-1-5 11:12
老师，我将这两种温度设得不一样，但MD也跑完了，这种情况没事吧

没事

作者
Author: alonewolfyang 时间: 2021-1-5 13:51

sobereva 发表于 2021-1-5 13:12
没事

谢谢老师

作者
Author: fineren 时间: 2021-1-25 10:35

sobereva 发表于 2020-12-29 11:10
最终要算什么温度，热浴参考温度就设多少。

产生初速度对应的温度不要设得太高，对常温产生初速度还好 ...

请问如何退火？我想从室温升至3000K，需要写脚本吗？

作者
Author: 量化小菜鸡 时间: 2021-1-26 14:43
社长，我没记错的话Gaussian的AIMD从过渡态结构开始跑可以用VENUS 96根据过渡态的振动模式计算起始速度。请问同样的方式容易在ORCA中实现吗？

作者
Author: sobereva 时间: 2021-1-26 19:20

量化小菜鸡发表于 2021-1-26 14:43
社长，我没记错的话Gaussian的AIMD从过渡态结构开始跑可以用VENUS 96根据过渡态的振动模式计算起始速度。请 ...

ORCA也可以自己在mdrestart文件里设初速度

作者
Author: sobereva 时间: 2021-1-26 19:21

fineren 发表于 2021-1-25 10:35
请问如何退火？我想从室温升至3000K，需要写脚本吗？

ORCA的MD部分手册里搜ramp

作者
Author: fineren 时间: 2021-1-27 10:41

sobereva 发表于 2021-1-26 19:21
ORCA的MD部分手册里搜ramp

谢社长

作者
Author: CYQ@ 时间: 2021-2-4 12:58
sob老师，我想请教一下：我用xtb跑了您这篇文章里的例子，我用gaussview构建了这个复合物，然后用xtb跑omd,发现它在进行分子动力学模拟前的优化中NH3向带正电的[Al(H2O)6]3+逐渐靠近，水上的一个质子转移到了氨气分子上，然后其生成的md轨迹文件是NH4+逐渐远离[Al(H2O)5(OH)]2+；而后我用xtb仅仅只对复合物做优化，发现其优化轨迹跟这篇文章过程相似；但是如果我用xtb对这个复合物不做优化，仅仅只做md的话，生成的轨迹文件是NH3一直在带正电的[Al(H2O)6]3+附近跑动,xtb的md是不能模拟这种质子转移的现象吗？

作者
Author: sobereva 时间: 2021-2-4 15:55

CYQ@ 发表于 2021-2-4 12:58
sob老师，我想请教一下：我用xtb跑了您这篇文章里的例子，我用gaussview构建了这个复合物，然后用xtb跑omd, ...

xtb可以（除非你用的GFN-FF）
我怀疑你可能没关闭xtb在MD中默认用的SHAKE约束算法，这会导致化学键键长始终固定不变

作者
Author: CYQ@ 时间: 2021-2-4 22:39

sobereva 发表于 2021-2-4 15:55
xtb可以（除非你用的GFN-FF）
我怀疑你可能没关闭xtb在MD中默认用的SHAKE约束算法，这会导致化学键键长 ...

谢谢sob老师

作者
Author: 84015917 时间: 2021-4-19 17:00
请问下ORCA4.1.1可以做这个吗？按照老师的例子，简单搞了一个自己想算的。但是一提交任务就停止了。不知道是不是版本的原因？

作者
Author: 84015917 时间: 2021-4-19 21:04

84015917 发表于 2021-4-19 17:00
请问下ORCA4.1.1可以做这个吗？按照老师的例子，简单搞了一个自己想算的。但是一提交任务就停止了。不知道 ...

自己琢磨了一下发现是 constraint add center 0..14 # Fix center of mass at initial position这行的命令有问题，删除了就行了，但是不知道是什么原因

作者
Author: sobereva 时间: 2021-4-21 15:31

84015917 发表于 2021-4-19 21:04
自己琢磨了一下发现是 constraint add center 0..14 # Fix center of mass at initial position这行的命 ...

帖子内容是对于发帖时最新版本而言的

作者
Author: 84015917 时间: 2021-4-25 16:49

sobereva 发表于 2021-4-21 15:31
帖子内容是对于发帖时最新版本而言的

谢谢老师，后续我跑了一个钙钛矿分子和一个水分子的动力学，然后发现水分子跑飞了如下图。是因为没有设置cell吗？

作者
Author: sobereva 时间: 2021-4-25 17:57

84015917 发表于 2021-4-25 16:49
谢谢老师，后续我跑了一个钙钛矿分子和一个水分子的动力学，然后发现水分子跑飞了如下图。是因为没有设置 ...

钙钛矿是无机晶体，不可能用分子模型
如果研究钙钛矿表面水分子的行为，应该用CP2K。用ORCA虽然也可以自己弄一个板来跑，但速度慢不说，还没法考虑周期边界条件，还得处理边界，很麻烦。

作者
Author: 84015917 时间: 2021-4-26 08:36

sobereva 发表于 2021-4-25 17:57
钙钛矿是无机晶体，不可能用分子模型
如果研究钙钛矿表面水分子的行为，应该用CP2K。用ORCA虽然也可以自 ...

好的，感谢老师

作者
Author: yqiusheng 时间: 2021-4-26 15:27
本帖最后由 yqiusheng 于 2021-4-26 15:32 编辑

我在用AIMD计算碳酸二甲酯在733K下的水解，水分子在5000步后经常跑飞，试了 cell fixed, cell pressure XXXX, cell elastic , 以上参数设为操作手册的默认值，但是还是跑飞，请问还有什么办法，还是我设置有问题？求解

作者
Author: sobereva 时间: 2021-4-26 16:01

yqiusheng 发表于 2021-4-26 15:27
我在用AIMD计算碳酸二甲酯在733K下的水解，水分子在5000步后经常跑飞，试了 cell fixed, cell pressure X ...

用cell sphere 施加球型限制势

作者
Author: yqiusheng 时间: 2021-4-26 16:25

sobereva 发表于 2021-4-26 16:01
用cell sphere 施加球型限制势

好的，谢谢！

作者
Author: 466468213 时间: 2021-4-28 20:34
社长，请问这种体系可以用Cp2k跑吗？

作者
Author: sobereva 时间: 2021-4-29 02:11

466468213 发表于 2021-4-28 20:34
社长，请问这种体系可以用Cp2k跑吗？

孤立体系用CP2K跑非常吃亏

作者
Author: hzfish 时间: 2021-5-10 15:16
(, 下载次数 Times of downloads: 101)

这种反应的过渡态，反应过程能不能使用ORCA做从头算动力学大致预测一下？
谢谢！

作者
Author: wzkchem5 时间: 2021-5-10 17:03

hzfish 发表于 2021-5-10 08:16
这种反应的过渡态，反应过程能不能使用ORCA做从头算动力学大致预测一下？
谢谢！

看势垒多高，如果从势垒看，ps到ns级的时间内能反应完，就可以做AIMD

作者
Author: hzfish 时间: 2021-5-10 17:55

wzkchem5 发表于 2021-5-10 17:03
看势垒多高，如果从势垒看，ps到ns级的时间内能反应完，就可以做AIMD

TS 4B势垒为75.9 kJ/mol，根据经验大约都跑多少ps或ns？
谢谢！

作者
Author: wzkchem5 时间: 2021-5-10 19:41

hzfish 发表于 2021-5-10 10:55
TS 4B势垒为75.9 kJ/mol，根据经验大约都跑多少ps或ns？
谢谢！

那没戏了。。。你这个得是ms到s量级的了
需要跑多久用Eyring方程换算一下就知道了

作者
Author: sobereva 时间: 2021-5-11 05:50

hzfish 发表于 2021-5-10 17:55
TS 4B势垒为75.9 kJ/mol，根据经验大约都跑多少ps或ns？
谢谢！

这种高度的势垒，必须用颇高的初始动能，并且恰到好处的碰撞角度（需要大量尝试）来跑AIMD，才可能得到想要的轨迹
光是把物质放一起，热浴设个实际反应温度就跑AIMD，猴年马月也跑不出来

作者
Author: sobereva 时间: 2021-7-7 13:29
新增加了一段

2021-Jul-7 针对ORCA 5.0的补充说明：对于2021-Jul-7及以后更新的Multiwfn，进入本文所述的Multiwfn产生输入文件的界面后，可以通过选项-11选择适合的ORCA版本，默认为ORCA 5.0及以后版本，而下文内容对应的是4.2.x版的情况。对于>=5.0版，Multiwfn自动设的热浴是CSVR，比之前版本唯一支持的Berendsen热浴更好，而且同样普适。而且在run后面多加了CenterCOM，这是从5.0版本开始支持的消除整体质心运动的选项，而下文里提到的constraint add center这一行就没有了。

作者
Author: Zhangpc 时间: 2021-7-10 12:07
用ORCA4.2.1(openMPI3.1.4)跑您的这个例子，发现正常跑了7步报错，好像是MPI出错：
[file orca_tools/qcmsg.cpp, line 458, Process 0]:
.... aborting the run

--------------------------------------------------------------------------
mpirun has exited due to process rank 0 with PID 117241 on
node cn23 exiting improperly. There are three reasons this could occur:

1. this process did not call "init" before exiting, but others in
the job did. This can cause a job to hang indefinitely while it waits
for all processes to call "init". By rule, if one process calls "init",
then ALL processes must call "init" prior to termination.

2. this process called "init", but exited without calling "finalize".
By rule, all processes that call "init" MUST call "finalize" prior to
exiting or it will be considered an "abnormal termination"

3. this process called "MPI_Abort" or "orte_abort" and the mca parameter
orte_create_session_dirs is set to false. In this case, the run-time cannot
detect that the abort call was an abnormal termination. Hence, the only
error message you will receive is this one.

This may have caused other processes in the application to be
terminated by signals sent by mpirun (as reported here).

You can avoid this message by specifying -quiet on the mpirun command line.
--------------------------------------------------------------------------
[file orca_tools/qcmsg.cpp, line 458]:
.... aborting the run

作者
Author: wzkchem5 时间: 2021-7-10 15:27

Zhangpc 发表于 2021-7-10 05:07
用ORCA4.2.1(openMPI3.1.4)跑您的这个例子，发现正常跑了7步报错，好像是MPI出错：
[file orca_tools/qcms ...

这些错误信息没有信息量，应该再往上多贴一些。此外把.int.log、.scf.log和.grad.log的最后几十行也贴上来

作者
Author: sobereva 时间: 2021-7-11 12:58

Zhangpc 发表于 2021-7-10 12:07
用ORCA4.2.1(openMPI3.1.4)跑您的这个例子，发现正常跑了7步报错，好像是MPI出错：
[file orca_tools/qcms ...

输入输出文件都发上来

作者
Author: sun666 时间: 2021-10-18 13:16
请问社长，这么多从头算动力学的软件。作为一个刚想接触AIMD的小白不知如何去选择软件。ORCA做AIMD的优势在哪里呢。CP2K是不是基本所有的体系都可以做？

作者
Author: 喵星大佬 时间: 2021-10-18 15:04
本帖最后由喵星大佬于 2021-10-19 09:55 编辑

sun666 发表于 2021-10-18 13:16
请问社长，这么多从头算动力学的软件。作为一个刚想接触AIMD的小白不知如何去选择软件。ORCA做AIMD的优势在 ...

ORCA的最大优势就是快，由于有RI加持，在不计算Hessian的情况下比同类软件快得多，而做BOMD用一般的走步方法只需要梯度即可。但是ORCA不能用于周期性体系，对于必须在周期性下处理的问题，需要使用CP2K/VASP/QE/CPMD等软件，CP2K也可以处理孤立体系，但是功能远不如ORCA强大和好用，此时速度也没有优势

作者
Author: sun666 时间: 2021-10-18 21:25

喵星大佬发表于 2021-10-18 15:04
ORCA的最大优势就是快，由于有RI加持，在不计算Hessian的情况下比同类软件快得多，而做BOMD用一般的走步 ...

谢谢回复，溶液中的反应应该算是孤立体系吧，我是想用来做溶液中的反应，有机萃取剂分子取代水合金属离子配位层中的水这个过程。此时是不是用ORCA应该更加合适吧。

作者
Author: sobereva 时间: 2021-10-19 01:32

sun666 发表于 2021-10-18 21:25
谢谢回复，溶液中的反应应该算是孤立体系吧，我是想用来做溶液中的反应，有机萃取剂分子取代水合金属离子 ...

是

作者
Author: ene 时间: 2021-10-19 03:08

sun666 发表于 2021-10-18 21:25
谢谢回复，溶液中的反应应该算是孤立体系吧，我是想用来做溶液中的反应，有机萃取剂分子取代水合金属离子 ...

这种情况不见得是孤立体系能够描述的，因为显式溶剂效应对体系行为的影响不一定可以忽略。这个时候用cp2k跑周期性下的包含显式溶剂的动力学显然更合理，而且用cp2k的话你还可以用更好的方式做控温控压（orca里面只有个半残废的弹性边界条件控压，这东西能不能产生正确的等温等压系综分布还不一定）。

作者
Author: 霜晨月 时间: 2021-10-28 22:14
本帖最后由霜晨月于 2021-10-28 22:22 编辑

sobereva 发表于 2021-5-11 05:50
这种高度的势垒，必须用颇高的初始动能，并且恰到好处的碰撞角度（需要大量尝试）来跑AIMD，才可能得到想 ...

不能用SMD或伞形采样吗？是不是ORCA软件不支持？

作者
Author: ene 时间: 2021-10-29 00:32

霜晨月发表于 2021-10-28 22:14
不能用SMD或伞形采样吗？是不是ORCA软件不支持？

orca5可以做metadynamics，你懂点数据后处理的话还能做伞型采样。

作者
Author: tiandikuoyuan 时间: 2021-11-12 14:14
老师，请问计算激发态动力学，是直接添加tddft模块吗？
计算s1结构的动力学，输入文件这样写对吗？

!B3LYP DEF2-SVP noautostart miniprint nopop
%maxcore 2000
%pal nprocs 24 end
%md
#restart ifexists # Continue MD by reading [basename].mdrestart if it exists. In this case "initvel" should be commented
#minimize # Do minimization prior to MD simulation
timestep 1.0_fs # This stepsize is safe at several hundreds of Kelvin
initvel 300_K no_overwrite # Assign velocity according to temperature for atoms whose velocities are not available
thermostat CSVR 300_K timecon 30.0_fs # Target temperature and coupling time constant
dump position stride 1 format xyz filename "pos.xyz" # Dump position every "stride" steps
#dump force stride 1 format xyz filename "force.xyz" # Dump force every "stride" steps
#dump velocity stride 1 format xyz filename "vel.xyz" # Dump velocity every "stride" steps
#dump gbw stride 20 filename "wfn" # Dump wavefunction to "wfn[step].gbw" files every "stride" steps
run 2000 CenterCOM # Number of MD steps. Remove motion of center of mass
end
%TDDFT NROOTS 3
IROOT 1
end

复制代码

作者
Author: wzkchem5 时间: 2021-11-12 17:03

tiandikuoyuan 发表于 2021-11-12 07:14
老师，请问计算激发态动力学，是直接添加tddft模块吗？
计算s1结构的动力学，输入文件这样写对吗？

对，这样可以

作者
Author: 胡说 时间: 2022-7-3 14:00
本帖最后由胡说于 2022-7-3 15:07 编辑

老师你好，ORCA 5.0.3试了wB97M-V和wB97X-V这两普通泛函，虽然会进入算梯度模块，但直接一开始就报错（cannot find the xc-energy file）。感觉是不是ORCA的AIMD模块本身还不支持，手册上也没找到相关说明。

作者
Author: wzkchem5 时间: 2022-7-3 15:33

胡说发表于 2022-7-3 07:00
老师你好，ORCA 5.0.3试了wB97M-V和wB97X-V这两普通泛函，虽然会进入算梯度模块，但直接一开始就报错（cann ...

贴出完整输入输出文件

作者
Author: sobereva 时间: 2022-7-4 01:16

胡说发表于 2022-7-3 14:00
老师你好，ORCA 5.0.3试了wB97M-V和wB97X-V这两普通泛函，虽然会进入算梯度模块，但直接一开始就报错（cann ...

本来就没必要用这俩泛函做AIMD

作者
Author: 胡说 时间: 2022-7-7 10:48

wzkchem5 发表于 2022-7-3 15:33
贴出完整输入输出文件

老师好，输入输出见附件。我发现这两个泛函和其它泛函不一样，它们在算梯度的时候会创建一个scfgrad.inp的文件，但其它泛函是没有的。

作者
Author: 胡说 时间: 2022-7-7 10:49

sobereva 发表于 2022-7-4 01:16
本来就没必要用这俩泛函做AIMD

是的，仅测试下，发现耗时比双杂化还高，可能是体系较难收敛。

作者
Author: 胡说 时间: 2022-7-7 10:53
另外, 请教一个关于速度初始化的问题：为什么第0步log里算出来的温度，不是输入文initvel设定的温度，而是要小得多，如下图：

作者
Author: yaoyuan0711 时间: 2022-7-7 11:12
请问sob老师，H+转移的动图是什么软件做出来的？

作者
Author: wzkchem5 时间: 2022-7-7 14:51

胡说发表于 2022-7-7 03:48
老师好，输入输出见附件。我发现这两个泛函和其它泛函不一样，它们在算梯度的时候会创建一个scfgrad.inp ...

看起来似乎是orca的bug，我有空看一下具体是怎么回事。
scfgrad.inp所有泛函都会产生，但是算完梯度以后会删掉。所以当且仅当梯度计算报错的时候会看到这个文件存留下来

作者
Author: wzkchem5 时间: 2022-7-7 15:34

胡说发表于 2022-7-7 03:48
老师好，输入输出见附件。我发现这两个泛函和其它泛函不一样，它们在算梯度的时候会创建一个scfgrad.inp ...

找到原因了。
像wB97M-V这种带VV10项的泛函，有两种计算方式：只对wB97M部分做SCF，最后算一个V的单点加到wB97M结果上（关键词叫做NL），以及对wB97M-V整体做SCF（关键词叫做SCNL）。NL比较快但是没有解析梯度；SCNL更慢，精度稍微高一点，有解析梯度。
程序本来的设计是，对于单点任务默认用NL，对于用到梯度的任务默认用SCNL。但是程序没考虑到MD也需要梯度，所以对于MD没有开SCNL，导致报错。所以解决方法是在!开头的那一行加一个关键词SCNL就行了。
下一版orca会把这个地方改过来，到时候即使不手动加SCNL，程序也能正确处理。不过正如sob老师所说，带VV10项的泛函跑AIMD性价比不好，建议改用相应的DFT-D3版本，快很多，精度基本不变。

作者
Author: 胡说 时间: 2022-7-7 16:53

wzkchem5 发表于 2022-7-7 15:34
找到原因了。
像wB97M-V这种带VV10项的泛函，有两种计算方式：只对wB97M部分做SCF，最后算一个V的单点加 ...

原来如此，多谢老师。
确实它们效率比wb97x-d3慢了一倍还要多，太不划算了。

另外，我用同样的体系也测试了双杂化泛函b2plyp,dsd-pbep86,wb97x-2。发现只有b2plyp跑出来的轨迹是合理的，也和wb97x-d3一致，其它两个泛函的结果很不合理，虽然梯度是有正常算的，体系位置却没什么大变化，这很奇怪，不知道是泛函本身如此，还是程序问题。我附上了相应的输入输出，老师如果有空的话，也可以看下怎么回事。

作者
Author: wzkchem5 时间: 2022-7-7 17:37

胡说发表于 2022-7-7 09:53
原来如此，多谢老师。
确实它们效率比wb97x-d3慢了一倍还要多，太不划算了。

可能只是另外两个泛函级别下势垒稍微高一点，导致跑MD的时间尺度下观察不到变化。势垒高1kcal/mol，速率要慢5倍左右，所以在跑MD时长不变的情况下，势垒稍微高一点就会让一个本来能看到的反应变成不反应。

作者
Author: 胡说 时间: 2022-7-7 21:30

wzkchem5 发表于 2022-7-7 17:37
可能只是另外两个泛函级别下势垒稍微高一点，导致跑MD的时间尺度下观察不到变化。势垒高1kcal/mol，速率 ...

这倒有可能，我设的温度比较低，可能没跑过去。谢谢老师。

作者
Author: sobereva 时间: 2022-7-8 05:19

yaoyuan0711 发表于 2022-7-7 11:12
请问sob老师，H+转移的动图是什么软件做出来的？

VMD的movie maker导出每一帧的图像文件，然后ffmpeg合并成动画
图省事直接用oCam在VMD播放过程中录屏也可以

作者
Author: yaoyuan0711 时间: 2022-7-8 08:27

sobereva 发表于 2022-7-8 05:19
VMD的movie maker导出每一帧的图像文件，然后ffmpeg合并成动画
图省事直接用oCam在VMD播放过程中录屏也 ...

谢谢老师！

作者
Author: 胡说 时间: 2022-7-19 15:05

胡说发表于 2022-7-7 10:53
另外, 请教一个关于速度初始化的问题：为什么第0步log里算出来的温度，不是输入文initvel设定的温度，而是 ...

发现去掉CenterCOM这个关键词那么第0步的Temp就等于Initvel指定的温度。如果加CenterCOM这个关键词，那么第0步动能对应的温度就小于指定的温度。
照理说既然在initvel阶段已经remove了COM的速度，然后再scale到对应的温度，那么加了CenterCOM为什么还会对第0步的速度有影响呢。CenterCOM不应该是对跑起来后才会有作用吗？感觉难以理解。

作者
Author: halouhapily 时间: 2022-8-26 10:45
老师，orca可以用混合基组吗？应该如何表示呢？
在orca支持镧系离子用Stuttgart赝势吗

作者
Author: sobereva 时间: 2022-8-26 13:44

halouhapily 发表于 2022-8-26 10:45
老师，orca可以用混合基组吗？应该如何表示呢？
在orca支持镧系离子用Stuttgart赝势吗

可
看手册
支持

作者
Author: 开心の土豆雷 时间: 2023-8-28 13:22
NH3分子的位置是怎么确定的？只要放在附近就会自动找到质子转移的轨迹嘛？

作者
Author: sobereva 时间: 2023-8-28 15:52

开心の土豆雷发表于 2023-8-28 13:22
NH3分子的位置是怎么确定的？只要放在附近就会自动找到质子转移的轨迹嘛？

放不远不近就完了

作者
Author: 开心の土豆雷 时间: 2023-8-28 19:53

sobereva 发表于 2023-8-28 15:52
放不远不近就完了

在调研AIMD的流程时，注意到初始结构的是一个非常关键的点。怎样根据自己的需要选择有意义的初始结构呢？
比如说我就想看一下两个不同分子之间能否发生反应以及怎样反应。从不同的初始结构开始跑AIMD肯定有不同的结论，这时候我该怎么选择初始结构呢？

作者
Author: sobereva 时间: 2023-8-28 20:02

开心の土豆雷发表于 2023-8-28 19:53
在调研AIMD的流程时，注意到初始结构的是一个非常关键的点。怎样根据自己的需要选择有意义的初始结构呢？ ...

容易发生的反应，随便放得不离谱就能发生。对于难发生的反应，各种相对位置和朝向、相对初速度都考虑。如果有静态计算得到的过渡态，可以参考着去摆

作者
Author: 开心の土豆雷 时间: 2023-8-28 21:36

sobereva 发表于 2023-8-28 20:02
容易发生的反应，随便放得不离谱就能发生。对于难发生的反应，各种相对位置和朝向、相对初速度都考虑。如 ...

谢谢sob老师，我手头是正好有过渡态的结构。

作者
Author: qaqfdmmj 时间: 2024-7-10 21:08
老师您好，我是根据您的帖子进行了类似工作的计算（仅将B97-3c换成了r2SCAN-3c），但是orca在输出每一步的时候都会有If you are NOT using ECPs, check your basis set inputs!和NOTIFICATION: Different basis set in ORCA and otool_gcp: ORCA: 619 gCP: 617的提示，请问这是为什么，是基组不适合我的体系吗？

作者
Author: sobereva 时间: 2024-7-10 23:22

qaqfdmmj 发表于 2024-7-10 21:08
老师您好，我是根据您的帖子进行了类似工作的计算（仅将B97-3c换成了r2SCAN-3c），但是orca在输出每一步的 ...

没事甭用r2SCAN-3c跑动力学，比B97-3c没明显额外好处还白多花时间

欢迎光临计算化学公社 (http://bbs.keinsci.com/)