计算化学公社

 找回密码 Forget password
 注册 Register
Views: 26650|回复 Reply: 25

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

  [复制链接 Copy URL]

231

帖子

4

威望

2815

eV
积分
3126

Level 5 (御坂)

发表于 Post on 2021-6-3 19:28:16 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 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:
  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的编译很简单,网上一搜即可,这里不细说了。
运行成功后会出现这样的窗口:
202106031746339150..png
这里是跑了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
复制代码

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

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






geo-opt.inp

2.34 KB, 下载次数 Times of downloads: 170

cif2pos

1.97 KB, 下载次数 Times of downloads: 123

x2c

1.23 KB, 下载次数 Times of downloads: 126

dimer.tar.gz

68.5 KB, 下载次数 Times of downloads: 211

dimer.inp

2.66 KB, 下载次数 Times of downloads: 204

idpp.py

1.25 KB, 下载次数 Times of downloads: 166

评分 Rate

参与人数
Participants 26
威望 +1 eV +97 收起 理由
Reason
qzm + 4 谢谢分享
ggdh + 5 谢谢
FTW + 4 谢谢分享
skiters + 1 好物!
1143310726 + 3
DZW + 5
Y-Chaos + 4
bajia + 2 谢谢
昼夏の忧郁 + 5 好活!!!!
思辰 + 4 好物!
ztchem + 5
ymygca + 5 谢谢
obaica + 5
biogon + 5
邓苏微 + 5 好物!
玉米猫 + 1 谢谢
丁越 + 5 谢谢
ljc050512 + 4 赞!
April_362 + 2 谢谢分享
乐平 + 5 好物!

查看全部评分 View all ratings

86

帖子

0

威望

654

eV
积分
740

Level 4 (黑子)

发表于 Post on 2021-8-5 11:14:08 | 显示全部楼层 Show all
好贴~

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

发表于 Post on 2021-8-21 08:26:46 | 显示全部楼层 Show all
老师的方法非常赞,使得找过渡态工作量大大简化了。但是编程小白还有个小小的建议,老师能不能把最后的绘制能量曲线的脚本简化成一个脚本呢?有点难记住,每次还得回来翻帖子
自由发挥,野蛮生长

231

帖子

4

威望

2815

eV
积分
3126

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-8-21 11:20:28 | 显示全部楼层 Show all
丁越 发表于 2021-8-21 08:26
老师的方法非常赞,使得找过渡态工作量大大简化了。但是编程小白还有个小小的建议,老师能不能把最后的绘制 ...

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

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

发表于 Post on 2021-8-23 08:25:14 | 显示全部楼层 Show all
djjj148 发表于 2021-8-21 11:20
设置alias即可,网上一搜就知道怎么操作了,我习惯调用脚本时只打2-3个字母以减轻手指负担,写全(cp2kdim ...

我滴天,我咋忘了可以设置alias呢谢谢老师
自由发挥,野蛮生长

1

帖子

0

威望

63

eV
积分
64

Level 2 能力者

发表于 Post on 2021-8-30 12:19:53 | 显示全部楼层 Show all
楼主用的是哪个版本的cp2k呀,安装不上去,编译不好

231

帖子

4

威望

2815

eV
积分
3126

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-8-30 12:40:47 | 显示全部楼层 Show all
sinco 发表于 2021-8-30 12:19
楼主用的是哪个版本的cp2k呀,安装不上去,编译不好

popt

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

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-9-3 10:01:48 | 显示全部楼层 Show all
djjj148 发表于 2021-8-30 12:40
popt

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

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

534

帖子

2

威望

2385

eV
积分
2959

Level 5 (御坂)

发表于 Post on 2021-9-11 00:14:01 | 显示全部楼层 Show all
老师,您好,请问你这变最终生成DIMER_VECTOR 文件是一个所有插点XYZ结构格式的文件吗?不知是否可以看一下您生成的DIMER_VECTOR文件的范例。
对于插点结构,我用的苯办法是, 利用反应物 和产物搭建QE的NEB输入文件,然后使用XCRYSDEN将各个插点导出成XSF文件,让后 用VESTA转换成xyz 文件,接着我就不知该怎么做了,所以想参考一下您的范例文件

231

帖子

4

威望

2815

eV
积分
3126

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-9-11 10:32:32 | 显示全部楼层 Show all
风飞 发表于 2021-9-11 00:14
老师,您好,请问你这变最终生成DIMER_VECTOR 文件是一个所有插点XYZ结构格式的文件吗?不知是否可以看一下 ...

是这样的。

DIMER_VECTOR

414 Bytes, 下载次数 Times of downloads: 35

534

帖子

2

威望

2385

eV
积分
2959

Level 5 (御坂)

发表于 Post on 2021-9-11 11:14:19 | 显示全部楼层 Show all

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

231

帖子

4

威望

2815

eV
积分
3126

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-9-11 12:48:25 | 显示全部楼层 Show all
风飞 发表于 2021-9-11 11:14
老师,您好, 我尝试使用:xyz2mode.py 0.xyz 5.xyz  生成DIMER_VECTOR  结果出现如下错误提示
请问该 ...

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

534

帖子

2

威望

2385

eV
积分
2959

Level 5 (御坂)

发表于 Post on 2021-9-11 16:17:00 | 显示全部楼层 Show all
djjj148 发表于 2021-9-11 12:48
你应该是用python2运行的脚本,帖子里有说明,需要用python3。

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

231

帖子

4

威望

2815

eV
积分
3126

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-9-11 17:30:52 | 显示全部楼层 Show all
风飞 发表于 2021-9-11 16:17
嗯,你好,请问如果安装python3,会不会对其他用python2程序有影响呢?

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

534

帖子

2

威望

2385

eV
积分
2959

Level 5 (御坂)

发表于 Post on 2021-9-11 20:58:58 | 显示全部楼层 Show all
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

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-2 21:46 , Processed in 0.232983 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list