计算化学公社

标题: 使用CP2K进行过渡态计算经验[dimer篇] [打印本页]

作者
Author: djjj148 时间: 2021-6-3 19:28
标题: 使用CP2K进行过渡态计算经验[dimer篇]
本帖最后由 djjj148 于 2021-6-4 15:03 编辑

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：

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

pip install pymatgen pymatgen-diffusion

复制代码

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

dist.pl POSCARis POSCARfs

复制代码

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

idpp.py POSCARis POSCARfs 4

复制代码

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

nebmovie.pl 0

复制代码

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

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

复制代码

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

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

复制代码

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

export PATH=$PATH:$HOME/cbin/dimer
export PYTHONPATH=$PYTHONPATH:$HOME/cbin/dimer

复制代码

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

xyz2mode.py 0.xyz 5.xyz

复制代码

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

cp POSCARis POSCAR
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的编译很简单，网上一搜即可，这里不细说了。
运行成功后会出现这样的窗口：
(, 下载次数 Times of downloads: 95)
这里是跑了200多步的DIMER，每一步的能量变化(单位是eV，我以前写的是a.u.这里没改过来，请忽略)，理想情况下能量都是单调递减的，像上图这种突然有一步上升了就要注意了，可能优化方向已经偏离。

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

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

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

复制代码

会输出类似的内容：
(, 下载次数 Times of downloads: 105)
这里显示的是第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的预-后处理脚本和程序，因为它们在结构处理上都是通用的。另外此贴可能不适合萌新，过渡态计算的基础知识公社里已经有很多帖子了，可以自行查看互补。

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

作者
Author: ljc050512 时间: 2021-8-5 11:14
好贴~

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

作者
Author: djjj148 时间: 2021-8-21 11:20

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

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

作者
Author: 丁越 时间: 2021-8-23 08:25

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

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

谢谢老师

作者
Author: sinco 时间: 2021-8-30 12:19
楼主用的是哪个版本的cp2k呀，安装不上去，编译不好

作者
Author: djjj148 时间: 2021-8-30 12:40

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

popt

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

作者
Author: Aridea 时间: 2021-9-3 10:01

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

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

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

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

作者
Author: djjj148 时间: 2021-9-11 10:32

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

是这样的。

作者
Author: 风飞 时间: 2021-9-11 11:14

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

老师，您好，我尝试使用：xyz2mode.py 0.xyz 5.xyz 生成DIMER_VECTOR 结果出现如下错误提示 (, 下载次数 Times of downloads: 97)
请问该怎么弄呢？

作者
Author: djjj148 时间: 2021-9-11 12:48

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

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

作者
Author: 风飞 时间: 2021-9-11 16:17

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

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

作者
Author: djjj148 时间: 2021-9-11 17:30

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

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

作者
Author: 风飞 时间: 2021-9-11 20:58

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

作者
Author: djjj148 时间: 2021-9-12 09:14

风飞发表于 2021-9-11 20:58
老师，您好，我安装了Anaconda，也安装了python3.7 请问这个目录怎么修改呢？
export PATH=$PATH:$H ...

这取决于你把dimer.tar.gz解压到了哪里，和用那个python版本没有关系

作者
Author: Riszzz 时间: 2021-10-14 10:19
您好，请问如果想通过CI-NEB和dimer结合的方法的话，DIMER_VECTOR文件该怎么设置呀

作者
Author: 丁越 时间: 2021-11-22 20:07
请教老师一个问题：我看了sob老师对dimer算法介绍的博文（http://sobereva.com/44），该算法有两个过程，旋转和平移过程。对于下图中第一个dimer过程结束后，这里的 ENERGY| Total FORCE_EVAL ( QS ) energy [a.u.]: -1593.699241516124857 应该就是体系在第一个dimer过程后的体系总能量了吧？可是在第一步的信息中（Informations at step）中，出现了个Total Energy = 0.0056990086，这个能量是啥呀？我咋不太理解，还请老师解答一下。

作者
Author: bajia 时间: 2021-11-23 16:09
老师，我进行dcc操作后没有出现dimmode.xyz，这个可以怎么修改代码得到，我是在虚拟机上进行操作的

作者
Author: seven111 时间: 2021-12-2 21:20
你好我想问一下，我在插点的这一步出现了这种问题，请问怎么解决啊？

作者
Author: seven111 时间: 2021-12-6 21:06

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

你好，请问你这个问题现在解决了吗？我也碰到了这个问题，但是我也安装了python3的呀，想问下你最终怎么解决这个问题的？

作者
Author: liudong123 时间: 2022-3-6 11:26
运行dcc脚本为啥有这个错误啊
File "/home/wangyang/hsw/320Li/test1tran/dimer/dimer/dcc", line 6
while getopts ":h" optname
^
SyntaxError: invalid syntax

作者
Author: bblovelp 时间: 2022-5-16 12:59

用DIMER计算了NiPd金属吸附解离氧气的过渡态，用30个核跑了3天才跑了14次迭代，感觉已经绝望！！

作者
Author: planet5460 时间: 2022-11-10 03:27
有点复杂哈，其实只要化学直觉较好，给出一个比较合适的初猜，基于这个初猜计算频率（固定住大多数原子，只放开参与反应的原子），然后选择虚频中需要的振动模式，保存为dimer_VECTOR即可。

作者
Author: h840473807 时间: 2022-12-1 23:37
请问CP2K能实现始末态晶胞参数不同的过渡态搜索吗,如果可以的话怎样载入始末态的晶胞参数呢?

作者
Author: qzm 时间: 2023-1-10 16:14
你好，我是CP2K新手，看了您的博文，想问下如何
使用脚本 cif2pos
生成POSCAR文件和CP2K的结构输入文件coord.inc。

作者
Author: sobereva 时间: 2023-2-13 05:41

qzm 发表于 2023-1-10 16:14
你好，我是CP2K新手，看了您的博文，想问下如何
使用脚本 cif2pos
生成POSCAR文件和CP2K的结构输入文件c ...

Multiwfn载入cif文件，主功能100的子功能2就能导出POSCAR
按下文就能产生CP2K输入文件，如果你想把坐标部分抠出来单独弄个文件就自己抠
使用Multiwfn非常便利地创建CP2K程序的输入文件
http://sobereva.com/587（http://bbs.keinsci.com/thread-21668-1-1.html）

作者
Author: sobereva 时间: 2023-2-13 05:52

h840473807 发表于 2022-12-1 23:37
请问CP2K能实现始末态晶胞参数不同的过渡态搜索吗,如果可以的话怎样载入始末态的晶胞参数呢?

原理上不行。变胞优化是专门的CELL_OPT任务，dimer所属的GEO_OPT任务不考虑变胞
但也可以考虑对前后结构取中点结构和晶胞参数去用dimer搜索过渡态结构当个近似，然后可以再冻结反应部分的原子做变胞优化，反复弄几个来回直到收敛

作者
Author: sobereva 时间: 2023-2-13 06:00

Riszzz 发表于 2021-10-14 10:19
您好，请问如果想通过CI-NEB和dimer结合的方法的话，DIMER_VECTOR文件该怎么设置呀

可以自己用TS相邻的两个点的结构求差，这近似对应于TS位置的反应路径的方向

作者
Author: sedmond 时间: 2023-3-11 23:37

sobereva 发表于 2023-2-13 06:00
可以自己用TS相邻的两个点的结构求差，这近似对应于TS位置的反应路径的方向

社长这个思路更贴合一些，楼主用的初态和终态获取vector，感觉是下猛药。另外，装了anaconda 1.11.0后（内置python3.9.13, perl 5.26.3），采用附件里提供的dcc时（想看初猜结合向量的效果），一直出现报错（/usr/bin/env: 'perl\r': Not a directory）。尝试采用”dos2unix *“，问题依旧，不知道大家有没有遇到这个情况。

作者
Author: lao7 时间: 2023-7-11 07:16
楼主，看到你的帖子受益匪浅！请问，如果过渡态计算完了，想做一个类似于高斯的scan，这个能做吗？那是怎么扫描呢？

作者
Author: joke122 时间: 2023-9-8 10:01
楼主请问一下，vtst脚本是可以单独使用还是必须和vasp绑定使用？

作者
Author: djjj148 时间: 2023-9-8 18:53

joke122 发表于 2023-9-8 10:01
楼主请问一下，vtst脚本是可以单独使用还是必须和vasp绑定使用？

前者

作者
Author: maoxinxina 时间: 2024-6-3 16:27
在运行dcc的时候报错dimer/dcc: 37: dimmode.pl: not found

作者
Author: maoxinxina 时间: 2024-6-7 09:30
DIMER计算结束后，跑个频率计算即可(该固定的原子记得固定)。可以直接用Multiwfn处理cp2k-1.restart生成频率计算的inp文件. 想问下计算结束后，不应该是跑最后一步的如cp2k-1_182.restart么。固定的原子应该是哪些呢。

作者
Author: lao7 时间: 2024-8-13 21:31
请教楼主：我看到您制作了初始文件和过渡态两个文件，如果仅仅用初始文件，稍微摆弄其接近过渡态是否可以？我的意思输入文件只有一个初始文件，而无过渡态文件。因为这个cp2k输入文件太复杂，用multiwfn只能生成包含一个结构的dimer计算文件。

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