Improved Dimer Method计算过渡态的脚本

发表于 更新于 分类于 阅读次数： Disqus：

今天分享一个Improved Dimer Method （IDM）的小脚本。IDM和传统的Dimer方法，NEB一样，都用于优化过渡态。Dimer方法的话，VTST（https://theory.cm.utexas.edu/vtsttools/）已经包含很多相关的脚本，咱们就不重复造车了。NEB结合IDM优化过渡态也很高效，成功率非常高，有兴趣的可以尝试下。分享脚本前主要介绍以下几点：

1） IDM已经编译在VASP里面了，不需要额外编译。

2） IDM跟平时的优化计算差不多，主要修改IBRION = 44即可。

3） IDM计算主要包括5步：

3.1： 生成初始的过渡态结构：可以有以下三种做法:

3.1.1 经验多的，可以直接手动初步搭一个，保存成POSCAR；

3.1.2 胆儿不大不小小的，可以借助NEB的方法，先准备初末态的结构，然后用VTST的nebmake.pl插点，挑出来接近过渡态结构的POSCAR；

3.1.3 胆儿最小的，可以先NEB跑上30步左右，把能量最高的结构检查下，感觉过渡态结构合理的话，把CONTCAR保存成POSCAR。

3.2：用3.1步中获得的POSCAR做一个频率计算；主要有以下几个注意事项：

3.2.1： NWRITE = 3 一定要注意，必须是3。

3.2.2: 这一步频率计算用不着那么高的精度，差不多就OK。比如我主要做表面催化反应，这一步计算的时候，

i) slab的原子都直接固定住，

ii) 分子大的话只放开对应键断开、生成的原子，分子其他部分都固定住。

iii) gamma点直接算。

这样的话，可以在很快的时间内完成频率计算，得到键断开、生成时过渡态所对应的虚频振动方向。如果有好几个小虚频，不用尝试着去消虚频，因为过渡态结构本来就是一个粗糙的，这里消虚频也没必要；不过一定要有过渡态对应的虚频，如果没有的话，重新搭一个过渡态结构，继续算频率。

3.3 读取频率计算的结果，生成IDM计算的POSCAR。这也是本节中脚本的作用；

3.4 准备IDM计算的INCAR, KPOINTS, POTCAR, 提交任务，等待结束。

3.5 精度稍微高点的频率计算，验证IDM优化出来的过渡态结构。

3.6 以上是正常计算的一个流程，具体细节，大家根据自己情况随机应变，官网的介绍只有一页，很短，一定要认真读一下：https://www.vasp.at/wiki/index.php/Improved_Dimer_Method

脚本部分：

脚本get_dimer.py可以在我的Github页面下载: https://github.com/bigbrosci/q-robot/tree/main/actions。 q-robot是本人这些年来计算常用的一些小脚本，有兴趣的也可以全部下载下来。根据自己需求酌情修改使用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: qli
Created on Thu Jun  2 22:09:18 2022
Creat POSCAR for improved dimer calculations from frequency calculations.
1) run frequency calculation: 
IBRION = 5
POTIM = 0.015
NFREE = 2
NWRITE = 3 ### Must be 3
2) run this script:
get_dimer.py 
3) use the POSCAR for IDM calc.
NSW = 100           
Prec=Normal
IBRION=44           !  use the dimer method as optimization engine
EDIFFG=-0.05
POTIM = 0.05
    
"""

import numpy as np
from ase.io import read, write
import os
from sys import exit

model = read('POSCAR_relax') ### POSCAR_relax is the POSCAR before freq calculations, that means some atoms are not fixed.
#model_positions = model.get_positions()
model.write('POSCAR_dimer', vasp5=True)

# print(model_positions)
# print(len(model))

l_start = 0 # the number of line which contains  'Eigenvectors after division by SQRT(mass)' 
with open('OUTCAR') as f_in:
    lines = f_in.readlines()
    for num, line in enumerate(lines):
        if 'Eigenvectors after division by SQRT(mass)' in line:
            l_start = num

if l_start == 0:
    print('''Check Frequency results and then rerun this script.\n**Remember**: NWRITE must be 3. BYEBYE! ''' )
    exit()
    
freq_infor_block = lines[l_start:]        
l_position = 0
wave_num = 0.0
for num, line in enumerate(freq_infor_block): 
        if 'f/i' in line:           
            wave_tem = float(line.rstrip().split()[6])
            if wave_tem > wave_num:
                wave_num = wave_tem  
                l_position = num+2

pos_dimer = open('POSCAR_dimer', 'a')
pos_dimer.write('  ! Dimer Axis Block\n')

vib_lines = freq_infor_block[l_position:l_position+len(model)]
for line in vib_lines:
    infor = line.rstrip().split()[3:]
    pos_dimer.write(' '.join(infor)+'\n')

pos_dimer.close()
print('''
      DONE!
      Output file is named as: POSCAR_dimer and can be used for dimer calculations.
      Don't forget to rename POSCAR_dimer to POSCAR before you run the dimer jobs.      
      ''')

需要注意的有以下几点：

1） 前面介绍了，频率计算为了图快，我们固定了一些原子。脚本会读取一个POSCAR_relax,也就是没有固定这些原子的结构。用ASE读取，输出一下，保证POSCAR里面的结构干净，避免出错。

2）脚本会读取最大的虚频，然后把它所对应的振动方向附到POSCAR_relax的后面，然后IDM的POSCAR就准备好了。

3）输出为POSCAR_dimer， 计算的时候记得把它改成POSCAR。

4）IDM任务的INCAR在脚本里面有说明，改下IBRION = 44。 剩下的基本跟平时优化是一样的。

5）直接在频率计算的目录里面运行：get_dimer.py 即可。

IDM结合NEB的成功率很高，非常推荐。另外，计算出的结果极少出现多个虚频的情况。因此，这个方法也可以用来消虚频。差不多就这些了，剩下的就是多操作练习，熟能生巧的事情了。欢迎大家分享自己的经验，放到自己的网站，公众号等平台，帮助更多的人解决VASP计算相关的问题，有兴趣的也可以把自己的经验心得放到大师兄网站上，直接给我发邮件即可（lqcata@gmail.com）。