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

本文最后更新于：2023/08/20 12:22

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

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

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

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

通过sed命令修改INCAR中的参数： 123456789101112def 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后会报错: 123456789101112131415161718 -----------------------------------------------------------------------------| || 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!!! = 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 ValueErrorexcept ValueError: print("Wrong input! Input should be positive integer!") sys.exit()if run_folder: folders = run_folderselse: folders = os.listdir()try: os.mkdir("freq")except FileExistsError: print("Attention: folder exists!")count = 0for 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:

“觉得不错的话，请我喝杯咖啡吧~ ୧(๑•̀⌄•́๑)૭”

微信支付

支付宝支付

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！