ASE消虚频的小脚本

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

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()