脚本分享(一)——消除VASP过渡态和稳态计算过程中的非理想虚频

type

Post

status

Published

date

Dec 11, 2022

slug

summary

分别介绍了可以用于消除VASP计算产生非理想虚频情况的几个脚本

category

技术分享

icon

password

URL

Property

Dec 3, 2023 06:35 AM

在使用VASP进行过渡态和稳态的结构优化后，进行零点能校正和吉布斯校正前需要进行频率计算。

查看频率时，可进入完成频率计算的路径，在终端中输入：


grep cm-1 OUTCAR

但是在进行频率计算后可能会产生我们意料之外的非理想虚频：

过渡态计算产生了两个及以上大于的虚频，或者没有虚频。稳态计算产生大于的虚频

之前，大师兄在他的微信公众号和书里分享了消除稳态计算虚频的脚本:

ASE消虚频的小脚本

ASE消虚频的小脚本。

https://mp.weixin.qq.com/s/J7KaUjclZoWB_XNecfpHJQ

如果有安装ase包，可以使用代码：


 # -*- coding: utf-8 -*-
 """
 Qiang Li
 Command: python3 vib_correct.py
 This is a temporary script file.
 """
 import numpy as np
 from ase.io import read, write
 import os

 #获取虚频开始行
 l_position = 0  #虚频振动方向部分在OUTCAR中的起始行数
 with open('OUTCAR') as f_in:
     lines = f_in.readlines()
     wave_num = 0.0
     for num, line in enumerate(lines):
         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
 # ASE读POSCAR
 model = read('POSCAR')
 model_positions = model.get_positions()
 num_atoms = len(model)
 #print(model_positions)

 # 获取虚频对应的OUTCAR部分
 vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行
 #print(vib_lines)
 vib_dis = []
 for line in vib_lines:
     #model_positions = [float(i) for i in line.rstrip().split()[:3]]
     vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列
     vib_dis.append(vib_infor)
 vib_dis = np.array(vib_dis) #将振动部分写到一个array中。

 # 微调结构
 new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子，即虚频对应振动位移的0.4，具体多大自己根据经验调。
 model.positions = new_positions

 ###保存结构
 write('POSCAR_new', model, vasp5=True) # POSCAR_new是微调后的结构，用于下一步的计算（别忘了把POSCAR_new改成POSCAR）。

问题在于，大师兄的脚本只能用于稳态的情况，但其实对于过渡态的情况，我们只需要对大师兄的脚本做一点小小的修改：

过渡态频率计算没有大于的虚频:

这种情况其实很简单，大师兄的脚本是消虚频，我们放大虚频就好了~~ 把大师兄的脚本这一行：


 # 微调结构
 new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子，即虚频对应振动位移的0.4，具体多大自己根据经验调。

中的+改成-即可，但是请记得做可视化，确定虚频方向是反应的方向！！！

过渡态频率计算产生了两个及以上大于的虚频：

这种情况我们希望把最大虚频以外的虚频优化到以下，我们只需要把大师兄脚本里需要消除的虚频改成次最大虚频即可：


# -*- coding: utf-8 -*-
import numpy as np
from ase.io import read, write
import os

#获取次最大虚频开始行
l_position = 0  #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
    lines = f_in.readlines()
    wave_num = 0.0
    wave_tem = 0.0
    for num, line in enumerate(lines):
        if 'f/i' in line:
            exwave_tem = wave_tem
            wave_tem = float(line.rstrip().split()[6])
            if wave_tem > wave_num:  #最大的虚频
                wave_num = exwave_tem
                l_position = num + 2
# ASE读POSCAR
model = read('POSCAR')
model_positions = model.get_positions()
num_atoms = len(model)
#print(model_positions)

# 获取次最大虚频对应的OUTCAR部分
vib_lines = lines[l_position - num_atoms - 3:l_position - 3]  #次最大虚频的振动部分
#vib_lines = lines[l_position:l_position + num_atoms]  #最大虚频的振动部分
#print(vib_lines)
vib_dis = []
for line in vib_lines:
    #model_positions = [float(i) for i in line.rstrip().split()[:3]]
    vib_infor = [float(i)
                 for i in line.rstrip().split()[3:]]  # dx， dy， dz对应的那三列
    vib_dis.append(vib_infor)
vib_dis = np.array(vib_dis)  #将振动部分写到一个array中。

# 微调结构
new_positions = model_positions + vib_dis * 0.4
# 0.4是微调的校正因子，即虚频对应振动位移的0.4，具体多大自己根据经验调。（消除虚频+）
model.positions = new_positions

###保存结构
write('POSCAR_new', model,
      vasp5=True)  # POSCAR_new是微调后的结构，用于下一步的计算（别忘了把POSCAR_new改成POSCAR）。

代码放在了我的github中，或者可以直接下载下面这个文件：

vib_correct_ts.py

1.5KB

用得开心~