使用CP2K进行过渡态计算经验[dimer篇]

djjj148

1. 前言

说到过渡态计算，很多人可能会更倾向于用VASP，QE，但是速度慢是硬伤，且对于大体系更是太吃资源。相比之下，CP2K免费，速度极快。撰写本文的初衷即是想让更多人感受到CP2K的强大。主要是由于我的计算资源极其有限(组里平均一个人只能分到100核)，对于100原子以上的体系算过渡态的话我就基本不会考虑用平面波程序了。而使用CP2K我就敢尝试到300原子体系，各种显式溶剂模型都敢往里塞，反正CP2K收敛性好，速度快。56核的单节点跑一个dimer的过渡态一般7-24个小时就能收敛(250原子体系)，就算用24核的小节点，也是时间多了一倍左右而已。
言归正传，此贴会分享用DIMER跑过渡态任务的整个过程，包括：初末结构的准备，过渡态结构的准备和验证，过渡态任务的实时监控，过渡态的验证。每个过程用到的脚本也会分享出来或留下出处。
PS： 可能有人会问为什么不介绍主流的CI-NEB方法。首先，我个人觉得CI-NEB太费时，就算是用CP2K。在计算资源相同下，插N个点的用时相当于跑N个DIMER计算，而且跑过渡态都用更稳健的CG方法，这是一个离子步需要多个SCF收敛的，相当于再一次扩大了CI-NEB和DIMER的用时差距。再者，我早期也用CI-NEB跑几步得到较好的DIMER任务的初猜结构再用DIMER精收敛。后面发现只要初猜得当(下文会详细讲)，直接跑DIMER就行了(起码对于我目前所研究的对象是这样的)。如果直接上DIMER发现计算实在无法步入正轨，也是推荐用CI-NEB+DIMER组合的，但不太建议一上来就全程CI-NEB(土豪请随意)。

2.1 初末结构的准备
我个人喜欢把inp文件和结构文件分开，这样方便单独对结构文件进行处理，本文的操作都是基于这样的。如果用Multiwfn生成一体化的inp也比较方便，届时请自行更改我的脚本。涉及到的python脚本都是基于python3的。
准备：配置好Multiwfn，VASPKIT
首先建好反应物和产物的初猜结构is.cif和fs.cif文件放在不同目录，借助Multiwfn和shell处理为VASP格式的POSCAR文件，比如对与is.cif：

1. cif2pos is.cif

复制代码

这会生成POSCAR文件和CP2K的结构输入文件coord.inc。计算的话这里附上一个输入文件的模板geo-opt.inp，就不详细说了。
计算结束后在两个计算目录分别运行x2c脚本，会调用VASPKIT和shell生成last和last.cif文件，分别是CP2K优化的最后一帧结构的POSCAR格式和cif格式文件。

2.2 过渡态结构的准备和验证
准备：
1. 安装pymatgen和pymatgen-diffusion

1. pip install pymatgen pymatgen-diffusion

复制代码

2. 去官网下载vtst脚本 http://theory.cm.utexas.edu/vtsttools/installation.html，解压后设置环境变量。注意运行有些脚本可能会出现报错，用dos2unix转一下一般能解决。
过渡态结构的准备：
新建一个目录，将is和fs目录的last文件分别移到该目录下并改名为POSCARis和POSCARfs。用vtst的dist.pl脚本判断要插的点数：

1. dist.pl POSCARis POSCARfs

复制代码

这会返回一个值，一般插点的数目可以是这个值除以0.8。
假设要插4个点。这里用许楠博士写的idpp.py非线性插点脚本，可得到比较理想的插点结构：

1. idpp.py POSCARis POSCARfs 4

复制代码

再用vtst的脚本把插点的结构合并成一个xyz文件：

1. nebmovie.pl 0

复制代码

生成movie.xyz文件，可用VMD查看插点的轨迹，并根据轨迹确定选哪一个结构作为DIMER计算的输入结构。比如想计算一个甲醇的亚甲基氢(H1)从C1上转移到一个表面模型的N1上的过渡态，就可以直接显示H1和C1的键长，以及H1和N1的键长，一般取两个键长大致相等的结构作为过渡态计算的初猜。或者对比is和fs的能量看看是吸热还是放热反应，灵活选取初猜。
接下来要准备DIMER_VECTOR文件，定义初始的dimer方向。先把插点的每个POSCAR文件分别复制成xyz文件：

1. for i in {0..5};do cp "0"$i/POSCAR.xyz $i".xyz";done

复制代码

假设选取第二个插点的结构(不算is和fs)用做计算：

1. cp 2.xyz coord.inc
   
2. sed -i '1,2d' coord.inc

复制代码

这里开始要用到赵亚帆博士写的一些python代码，年份较老，是基于python2的，我自己转成了python3并修改了一点点，加上自己写的一些shell脚本打包成了dimer.tar.gz。解压后，假设所在目录为$HOME/cbin/dimer，则需要设置环境变量和pythyon的环境变量：

1. export PATH=$PATH:$HOME/cbin/dimer
   
2. export PYTHONPATH=$PYTHONPATH:$HOME/cbin/dimer

复制代码

用反应物和产物生成DIMER_VECTOR文件：

1. xyz2mode.py 0.xyz 5.xyz

复制代码

至此，DIMER计算的所有文件都准备好了，可以用dcc脚本生成一个DIMER振动方向的xyz文件：

1. cp POSCARis POSCAR
   
2. dcc

复制代码

dcc中最后一行的sz命令是因为我自己用的是windows系统的Xshell，需要把Linux服务器上的文件传到windows用sz，不需要的可以直接注释掉，下同。
把dimmode.xyz文件用VMD打开即可观察到DIMER计算的方向。
确认振动方向是自己心中的反应方向无误后即可提交计算，这里传一个模板dimer.inp。

2.3 过渡态任务的实时监控
计算中需要监控能量和DIMER的方向，发现不对劲及时调整以免白白等待。
对于能量，直接cp2kdimerfe.py cp2k.out && transform_fe.py && dimer_opt_eV即可(cp2k.out是我指定的out文件名，不同的提交任务命令会不同)，注意这里会调用gnuplot画图，所以需要安装gnuplot。我个人很喜欢gnuplot，速度快。gnuplot的编译很简单，网上一搜即可，这里不细说了。
运行成功后会出现这样的窗口：

这里是跑了200多步的DIMER，每一步的能量变化(单位是eV，我以前写的是a.u.这里没改过来，请忽略)，理想情况下能量都是单调递减的，像上图这种突然有一步上升了就要注意了，可能优化方向已经偏离。

对于DIMER的方向，运行dm脚本即可，会生成dimmode.xyz文件，这是最后一步的DIMER振动方向的xyz文件，放到VMD里可视化即可。如果想查看其中一步的也可以，比如上图在191步能量升高了，可以运行dm2 cp2k-1_191.restart查看第191步的DIMER方向。
如果DIMER方向和能量都没问题，即可放心继续跑任务。

2.4 过渡态的验证
判断过渡态是否收敛，可以直接

1. grep "Informations at step =" cp2k.out -A20

复制代码

会输出类似的内容：

这里显示的是第214步的收敛信息，和Gaussian很像，4个收敛标准都是yes即收敛。
DIMER计算结束后，跑个频率计算即可(该固定的原子记得固定)。可以直接用Multiwfn处理cp2k-1.restart生成频率计算的inp文件。详见http://sobereva.com/587(使用Multiwfn非常便利地创建CP2K程序的输入文件)
可以grep -E "Freq" cp2k.out查看频率，目录下还出现了freq-VIBRATIONS-1.mol，里面用molden格式记录了振动矢量，可以用Molden程序观看振动模式。
如果振动模式连接反应物和产物，且只有一个较大的虚频，过渡态计算成功。

3. 总结
这里借用了很多VASP的预-后处理脚本和程序，因为它们在结构处理上都是通用的。另外此贴可能不适合萌新，过渡态计算的基础知识公社里已经有很多帖子了，可以自行查看互补。

============================================================
能力有限，请批评指正。

ljc050512

好贴~

丁越

老师的方法非常赞，使得找过渡态工作量大大简化了。但是编程小白还有个小小的建议，老师能不能把最后的绘制能量曲线的脚本简化成一个脚本呢？有点难记住，每次还得回来翻帖子

djjj148

丁越 发表于 2021-8-21 08:26
老师的方法非常赞，使得找过渡态工作量大大简化了。但是编程小白还有个小小的建议，老师能不能把最后的绘制 ...

设置alias即可，网上一搜就知道怎么操作了，我习惯调用脚本时只打2-3个字母以减轻手指负担，写全(cp2kdimerfe.py cp2k.out && transform_fe.py && dimer_opt_eV)是想让大家知道调用的具体脚本，以方便想改脚本的人。

丁越

djjj148 发表于 2021-8-21 11:20
设置alias即可，网上一搜就知道怎么操作了，我习惯调用脚本时只打2-3个字母以减轻手指负担，写全(cp2kdim ...

我滴天，我咋忘了可以设置alias呢谢谢老师

sinco

楼主用的是哪个版本的cp2k呀，安装不上去，编译不好

djjj148

sinco 发表于 2021-8-30 12:19
楼主用的是哪个版本的cp2k呀，安装不上去，编译不好

popt

安装问题公社有很多帖子详细讲了，可以搜索学习。

Aridea

djjj148 发表于 2021-8-30 12:40
popt

安装问题公社有很多帖子详细讲了，可以搜索学习。

哈哈，第一次看大佬帖子不解其中意，只给了1分，后来想给10分，追悔莫及。。

风飞

老师，您好，请问你这变最终生成DIMER_VECTOR 文件是一个所有插点XYZ结构格式的文件吗？不知是否可以看一下您生成的DIMER_VECTOR文件的范例。
对于插点结构，我用的苯办法是， 利用反应物 和产物搭建QE的NEB输入文件，然后使用XCRYSDEN将各个插点导出成XSF文件，让后 用VESTA转换成xyz 文件，接着我就不知该怎么做了，所以想参考一下您的范例文件

djjj148

风飞 发表于 2021-9-11 00:14
老师，您好，请问你这变最终生成DIMER_VECTOR 文件是一个所有插点XYZ结构格式的文件吗？不知是否可以看一下 ...

是这样的。

风飞

djjj148 发表于 2021-9-11 10:32
是这样的。

老师，您好， 我尝试使用：xyz2mode.py 0.xyz 5.xyz  生成DIMER_VECTOR  结果出现如下错误提示
请问该怎么弄呢？

djjj148

风飞 发表于 2021-9-11 11:14
老师，您好， 我尝试使用：xyz2mode.py 0.xyz 5.xyz  生成DIMER_VECTOR  结果出现如下错误提示
请问该 ...

你应该是用python2运行的脚本，帖子里有说明，需要用python3。

风飞

djjj148 发表于 2021-9-11 12:48
你应该是用python2运行的脚本，帖子里有说明，需要用python3。

嗯，你好，请问如果安装python3，会不会对其他用python2程序有影响呢？

djjj148

风飞 发表于 2021-9-11 16:17
嗯，你好，请问如果安装python3，会不会对其他用python2程序有影响呢？

使用恰当就不会，可以用Anaconda构建虚拟环境，搜索"Anaconda"即可。

风飞

djjj148 发表于 2021-9-11 17:30
使用恰当就不会，可以用Anaconda构建虚拟环境，搜索"Anaconda"即可。

老师，您好，我安装了Anaconda， 也安装了python3.7  请问  这个目录怎么修改呢？
export PATH=$PATH:$HOME/cbin/dimer
export PYTHONPATH=$PYTHONPATH:$HOME/cbin/dimer

我的Anaconda目录是home/room/anaconda3/bin