ASE消虚频的小脚本

ASE消虚频的小脚本

做计算,虚频是令我们又爱又恨的一个东西。结构优化的时候,出现虚频,往往需要我们花费大量的时间和精力去消掉它们,最后搞得精疲力尽。算过渡态的时候,它又不出来,又得折腾个半死才出现。再不幸一些,除了过渡态对应的虚频外,还出现了一个伴生的姐弟虚频,我们又得继续消下去。本节简单介绍一个结构优化时消虚频的小脚本。主要是根据虚频的振动方向将体系中对应原子的坐标进行微调,然后用于重新优化。

简单介绍下原理,先看下VASP的OUTCAR中虚频对应的部分,

1)7220行带有f/i说明有虚频存在,

2) X, Y, Z为体系中原子的坐标,也就是POSCAR中的内容。

3) dx, dy,dz为虚频对应的原子振动方向上的位移。

4)我们将(dx, dy,dz)这个位移乘以一个介于0-1之间的校正因子,然后跟POSCAR中的坐标加起来即可。

5)以z方向为例,校正因子设为0.1,微调后的z方向坐标为: 10.709980 + (-0.676634) $\times$ 0.1

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
7219
7220 3 f/i= 9.313471 THz 58.518262 2PiTHz 310.663930 cm-1 38.517432 meV
7221 X Y Z dx dy dz
7222 6.003159 11.629641 6.067545 0 0 0
7223 7.651085 3.894943 9.628260 0 0 0
7224 2.424442 5.785382 9.658882 0 0 0
7225 5.384679 7.755200 0.000000 0 0 0
7226 0.004928 9.643165 0.000000 0 0 0
7227 8.286347 7.755200 0.877073 0 0 0
7228 2.944132 9.643165 0.971480 0 0 0
7229 6.823299 9.643165 1.459420 0 0 0
7230 1.481083 7.755200 1.553827 0 0 0
7231 4.619917 7.755200 2.336493 0 0 0
7232 9.762503 9.643165 2.430900 0 0 0
7233 6.756153 7.724844 3.627352 0 0 0
7234 1.491401 9.684490 3.647528 0 0 0
7235 8.239374 9.666853 5.136018 0 0 0
7236 3.012611 7.740104 5.091365 0 0 0
7237 6.016826 7.737975 6.049788 0 0 0
7238 0.788563 9.668849 6.080188 0 0 0
7239 8.321833 7.730633 7.184644 0 0 0
7240 3.081015 9.671865 7.187374 0 0 0
7241 0.815711 7.728916 8.180983 0 0 0
7242 6.063762 9.648841 8.115270 0 0 0
7243 9.839838 9.660462 8.962566 0 0 0
7244 4.631258 7.730218 9.011537 0 0 0
7245 7.668770 7.725888 9.623941 0 0 0
7246 2.424191 9.657803 9.623832 0 0 0
7247 2.095436 3.874640 10.694551 0 0 0
7248 0.099192 11.590822 10.709980 -0.736317 0.002031 -0.676634

脚本内容 (命名为: vib_correct.py )

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
# -*- 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)。

脚本使用:

  1. 脚本跟计算文件在同一个目录下:python3 vib_correct.py
  1. 如果赋予了脚本可执行权限,并放到了~/bin目录下,直接运行: vib_correct.py 命令即可得到校正后的结构文件:POSCAR_new
  2. 重新计算的时候,新建一个文件夹,把POSCAR_new 复制到文件夹里面,记得把名字改成POSCAR
  3. 这个脚本只适用于普通结构优化时候的虚频。如果过渡态计算的时候出现好几个虚频,不要使用这个脚本,因为脚本默认使用最大虚频对应的位移进行校正。
  4. 不用ASE,自己直接写个python脚本也能实现。不过使用ASE可以很方便地获取POSCAR中原子数目,保存原子固定的信息。参考如下。
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
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
# 读POSCAR
with open('POSCAR') as f_pos:
lines_pos = f_pos.readlines()

# 获取虚频对应的OUTCAR部分
num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()])
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行

model_positions = []
vib_dis = []
for line in vib_lines:
position = [float(i) for i in line.rstrip().split()[:3]]
vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
model_positions.append(position)
vib_dis.append(vib_infor)

# 微调结构
model_positions = np.array(model_positions)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。

###保存结构
f_out = open('POSCAR_new','w')
f_out.writelines(lines_pos[:8])
f_out.write('Cartesian\n')
for i in new_positions:
f_out.write(' '.join([str(coord) for coord in i]) + ' F F F\n')
f_out.close()