python脚本批量进行VASP吸附结构频率计算

前言（Introduction）

在用VASP做催化反应机理分析时，自由能变化是最为重要的分析内容之一。通常需要进行频率计算来进行自由能校正^[1]。我们常计算的体系一般是某个(系列)结构吸附一些小的分子或者中间产物，这样的情况下频率计算需要以下几个步骤：

• 复制CONTCAR为POSCAR
• 修改POSCAR，fix吸附基底，relax表面吸附的分子
• 修改INCAR参数
  • IBRION = 5
  • NFREE = 2
  • POTIM = 0.015
  • EDIFF = 1E-7
• run

虽然步骤简单，但是当要进行一系列结构吸附某个中间产物的频率计算时，手动操作比较麻烦(特别是修改POSCAR进行固定的时候)，所以尝试用python脚本批量实现上述步骤，减少工作量。

方案 (Solution)

1. 前置准备

vaspkit 1.2.5 ~ 1.3.1 (update: vaspkit 1.3.5版402功能不支持fix no atom）
python 3.8.0

2. 代码思路

step 1

脚本的最关键在于实现自动修改POSCAR使得吸附基底fix，而吸附的分子或者中间体是relax。考虑到大部分情况下，我们需要relax的原子往往都是在最顶部的(也就是POSCAR中z值最大的那几个原子)，所以我们可以先fix POSCAR中所有原子(vaspkit 402功能-> 1 -> 2 -> 0 0 -> 2)，然后根据吸附分子/中间产物的个数n，relax POSCAR中最顶部的n个原子。比如吸附的小分子是CO₂，那我们输入3，就可以relax最顶部的3个原子。具体代码如下：

def relax_top_atoms(num):
    os.popen(
        "(echo 402; echo 1; echo 2; echo 0 0; echo 2;)|vaspkit").readlines()  # generate POSCAR_FIX with Cartesian coord
    atom_num_str = os.popen('sed -n 7p POSCAR_FIX | tail -n1').readline().split()
    atom_num = sum(list(map(int, atom_num_str)))
    if num > atom_num:
        print("Relax atom number greater than total atom number! Wrong num input or wrong POSCAR/CONTCAR!")
        sys.exit()
    os.system("sed -i 's/T  T  T/F  F  F/g' POSCAR_FIX")
    lines = os.popen("sed -n '10,$p' POSCAR_FIX").readlines()
    z_coord = []
    for line in lines:
        z_coord.append(line.split()[2])
    z = list(map(float, z_coord))
    index_list = sorted(range(len(z)), key=lambda k: z[k], reverse=True)
    index_to_relax = index_list[:num]
    for i in index_to_relax:
        z_relax = z_coord[i]
        os.system("sed -i 's/{}    F  F  F/{}    T  T  T/g' POSCAR_FIX".format(z_relax, z_relax))
    os.system("cp POSCAR_FIX POSCAR")
    print("Selected top atoms set relaxed.")

这里先用vaspkit 402功能对所有原子固定，然后把POSCAR中所有原子的z轴坐标读出，排序出最顶部原子所在的z轴的坐标，最后把F F F替换为T T T，实现最顶部n个原子的relax。 最关键的代码 index_list = sorted(range(len(z)), key=lambda k: z[k], reverse=True) 很有意思:rofl:

step 2

通过sed命令修改INCAR中的参数：

def con2pos(name):
    os.chdir(name)
    os.system("sed -i 's/EDIFF[^G].*/EDIFF = 1E-7/' INCAR")
    #os.system("sed -i 's/^\s*NPAR/#NPAR/' INCAR")
    os.system("sed -i 's/IBRION.*/IBRION = 5/' INCAR")
    os.system("sed -i 's/POTIM.*/POTIM = 0.015/' INCAR")
    os.system("sed -i '/POTIM/a\   NFREE = 2' INCAR")
    os.system("cp CONTCAR POSCAR")
    print("CONTCAR in {} copied to POSCAR".format(name))
    os.system("rm slurm*")
    relax_top_atoms(relax_atom_num)
    os.chdir(cwd)

这里有一点需要注意的是： os.system("sed -i 's/^\s*NPAR/#NPAR/' INCAR")这行代码把NPAR给注释掉了，因为我一般用VASP进行结构优化时ALGO参数一般是Fast，用该算法进行频率计算算了一步之后就会报错，只能注释掉。如果想频率计算的时候也可以设置NPAR提高并行计算效率，可以在结构优化的时候设置ALGO = Normal，这样的话频率计算时就不会报错。

Update: 最近发现以前做的一些计算，一般都不会设置ISYM，因此使用的是默认设置（ISYM = 2，也就是使用对称性），但是IBRION = 5时会打破对称性，因此设置NPAR后会报错:

 -----------------------------------------------------------------------------
|                                                                             |
|     EEEEEEE  RRRRRR   RRRRRR   OOOOOOO  RRRRRR      ###     ###     ###     |
|     E        R     R  R     R  O     O  R     R     ###     ###     ###     |
|     E        R     R  R     R  O     O  R     R     ###     ###     ###     |
|     EEEEE    RRRRRR   RRRRRR   O     O  RRRRRR       #       #       #      |
|     E        R   R    R   R    O     O  R   R                               |
|     E        R    R   R    R   O     O  R    R      ###     ###     ###     |
|     EEEEEEE  R     R  R     R  OOOOOOO  R     R     ###     ###     ###     |
|                                                                             |
|      VASP internal routines  have requested a change of the k-point set.    |
|      Unfortunately this is only possible if NPAR=number of nodes.           |
|      Please remove the tag NPAR from the INCAR file and restart the         |
|      calculations.                                                          |
|                                                                             |
|      ---->  I REFUSE TO CONTINUE WITH THIS SICK JOB ..., BYE!!! <----       |
|                                                                             |
 -----------------------------------------------------------------------------

出于这个目的，代码中会把NPAR一行注释掉, 但是这样的话，计算效率会低很多。现在只要设置ISYM = 0之后，就可以正常设置NPAR了。 参考链接：https://mattermodeling.stackexchange.com/questions/9004/ncore-bigger-than-1-when-performing-phonon-vibrational-calculationibrion-5-in-VASP

step 3

最后，增加一些让用户可以选择对哪些文件夹进行频率计算的代码和一些错误提示就完事了，完整代码如下：

import os
import sys

def make_dir(name):
    try:
        os.mkdir(name)
        return True
    except FileExistsError:
        print("{} exists, skipping...".format(name))
        return False


def vasp_folder(name):
    os.chdir(name)
    count = 0
    for i in os.listdir():
        if i in ["INCAR", "POSCAR", "POTCAR", "KPOINTS", "CONTCAR", "OUTCAR"]:
            count += 1
    os.chdir(cwd)
    if count >= 5:
        return True
    else:
        print(f"{name} does not seem to be a finished VASP job folder")


def con2pos(name):
    os.chdir(name)
    os.system("sed -i 's/EDIFF[^G].*/EDIFF = 1E-7/' INCAR")
    # os.system("sed -i 's/^\s*NPAR/#NPAR/' INCAR")
    os.system("sed -i 's/IBRION.*/IBRION = 5/' INCAR")
    os.system("sed -i 's/POTIM.*/POTIM = 0.015/' INCAR")
    os.system("sed -i '/POTIM/a\   NFREE = 2' INCAR")
    os.system("cp CONTCAR POSCAR")
    print("CONTCAR in {} copied to POSCAR".format(name))
    os.system("rm slurm*")
    relax_top_atoms(relax_atom_num)
    os.chdir(cwd)
    global count
    count += 1


def relax_top_atoms(num):
    os.popen(
        "(echo 402; echo 1; echo 2; echo 0 0; echo 2;)|vaspkit").readlines()  # generate POSCAR_FIX with Cartesian coord
    atom_num_str = os.popen('sed -n 7p POSCAR_FIX | tail -n1').readline().split()
    atom_num = sum(list(map(int, atom_num_str)))
    if num > atom_num:
        print("Relax atom number greater than total atom number! Wrong num input or wrong POSCAR/CONTCAR!"))
        sys.exit()
    os.system("sed -i 's/T  T  T/F  F  F/g' POSCAR_FIX")
    lines = os.popen("sed -n '10,$p' POSCAR_FIX").readlines()
    z_coord = []
    for line in lines:
        z_coord.append(line.split()[2])
    z = list(map(float, z_coord))
    index_list = sorted(range(len(z)), key=lambda k: z[k], reverse=True)
    index_to_relax = index_list[:num]
    for i in index_to_relax:
        z_relax = z_coord[i]
        os.system("sed -i 's/{}    F  F  F/{}    T  T  T/g' POSCAR_FIX".format(z_relax, z_relax))
    os.system("cp POSCAR_FIX POSCAR")
    print("Selected top atoms set relaxed.")


cwd = os.getcwd()
exclude_folder = input("Specif exclude folders & files (using space to seperate):")
exclude_folders = exclude_folder.split()
run_folder = input("Specif folders to run (default all, using space to seperate):")
run_folders = run_folder.split()
top_atoms = input("How many atoms to be set to relaxed on top?\n")
try:
    relax_atom_num = int(top_atoms)
    if relax_atom_num < 0:
        raise ValueError
except ValueError:
    print("Wrong input! Input should be positive integer!")
    sys.exit()
if run_folder:
    folders = run_folders
else:
    folders = os.listdir()
try:
    os.mkdir("freq")
except FileExistsError:
    print("Attention: folder <freq> exists!")
count = 0
for folder in folders:
    if folder in exclude_folders:
        continue
    if os.path.isfile(folder):
        continue
    if not vasp_folder(folder):
        continue
    if folder.split("_")[-1] == "freq":
        continue
    freq_name = folder + "_freq"
    freq_name = "freq/" + freq_name
    flag = make_dir(freq_name)
    if not flag:
        continue
    os.system("cp -r {}/* {}/".format(folder, freq_name))
    #if os.popen("grep '^\s*NPAR' {}/INCAR".format(freq_name)).read():
        #print("NPAR tag should be removed for most freq calculations! Will comment out NPAR...")
    con2pos(freq_name)

if count != 0:
    print("Remember to check if atoms are fixed properly in POSCAR before freq calculation!")

自己用了一段时间，基本上没有什么问题，欢迎各位使用测试:satisfied:

参考资料 (References)

1. 校正分子和吸附分子自由能. Jincheng Liu. http://blog.wangruixing.cn/2019/04/21/freenergy/ ↩︎