脚本分享(五)——VASP的简单工作流

type

Post

status

Published

date

Dec 3, 2023

slug

summary

VASP的工作流脚本

1. 引言

入学三个月了，研究生生活最难以适应的还是在实验组很难找到人可以一起讨论，超算之类的还要自己去找吧。最近忙着改本科的各种稿子，发现很久没有写一些自己的东西了。正好今天睡了难得的懒觉，写一篇博客，给自己一点正反馈。

写这几个脚本的原因是：节省时间！

其实VASP表面计算的主要流程：建模——准备输入文件——结构优化收敛——固定原子——计算频率做热力学校正。建模的工作量无法避免，所以想办法简化一点后面的工作。

现在在实验组，除了帮同门做计算，还得做自己的实验课题。上述每个步骤基本还是需要自己调整，时间上安排不过来，所以想整一个简单的工作流节省一下自己的时间。自己也没空去学aiida之类的，于是有了后面的一系列东西。

由于现在自己组里没有服务器，现在正在使用的是THU的超算探索1000，老板在里面充了1w，本着给老板省钱的想法，在结构优化这一步采用了粗收敛，精收敛的步骤，增加了不少工作量。

使用vaspkit的402功能固定原子可能无法准确固定基底，尤其是在缺陷位吸附之类的情况

来这边之后，除了配置服务器之外，写了下面的这几个脚本:


探索1000上提交任务的脚本：
svasp
job.vasp

准备输入文件的脚本：
new_miss_dir.py
auto_vibm.py
auto_fixatom.py

进行自动化工作流的脚本：
subvasp
check_convergence.sh
handle_INCAR_structural_relaxation_fine.py

2. 正文

首先介绍实现脚本需要的文件结构，之后介绍几个脚本的具体内容。

2.1 文件结构

对于一个系统的研究，通常需要对比多条路径，例如现在的四条：


.
├── 00moluar
├── 01GasPhasePath
├── 02GasPhasePath2
├── 03SurfacePath
└── 04H_Path

进入04H_Path文件夹，其结构为：


(VASP) [zhaoxb@ibcln01 04H_Path]$ tree -d -L 1
.
├── 01Ag_H1
├── 02Ag_H2
├── 03Ag_H3
├── 04Ag_H4
├── 05Ag_BHMF_OH
├── 06Ag_BHMF_CH
├── 07Ag_OH1
├── 08Ag_OH2
├── incar
└── vaspkit_command

由于每条路径基本可以使用相同的参数，所以将其可能会用到的INCAR以及vaspkit命令放到incar与vaspkit_command文件夹。


(VASP) [zhaoxb@ibcln01 incar]$ tree -L 1
.
├── INCAR_bader
├── INCAR_dimer_frequency_calculation
├── INCAR_frequency_calculation
├── INCAR_NEB
├── INCAR_structural_relaxation
└── INCAR_structural_relaxation_fine

而vaspkit_command文件夹记录了vaspkit需要使用的参数还有一个KPOINTS文件


(VASP) [zhaoxb@ibcln01 vaspkit_command]$ tree -L 1
.
├── KPOINTS
├── vaspkit_402
├── vaspkit_501
└── vaspkit_502

2.2 输入文件准备

2.2.1 new_miss_dir.py

在Material Studio上建模并生成POSCAR，此时是一系列POSCAR.txt的文件，将这些文件复制到服务器上。运行new_miss_dir.py脚本。（注意此时POSCAR第一行是在MS中设置的文件名也是脚本会创建的文件夹名称）


#! /usr/bin/env python3
# -*- coding: UTF-8 -*-
import os
import shutil
from pathlib import Path

def find_folder_in_ancestors(folder_name, depth=3):
    current_path = Path.cwd()
    for _ in range(depth):
        folder_path = current_path / folder_name
        if folder_path.is_dir():
            return folder_path
        current_path = current_path.parent
    return None

def main():
    for file in os.listdir('.'):
        if file.startswith('POSCAR') and file.endswith('.txt'):
            with open(file, 'r') as f:
                folder_name = f.readline().strip()
            new_folder_path = Path(folder_name)
            new_folder_path.mkdir(exist_ok=True)

            # Rename and move POSCAR.txt
            new_poscar_path = new_folder_path / 'POSCAR'
            shutil.move(file, new_poscar_path)

            # Copy INCAR
            incar_folder = find_folder_in_ancestors('incar')
            if incar_folder:
                shutil.copy(incar_folder / 'INCAR_structural_relaxation', new_folder_path / 'INCAR')

            # Copy KPOINTS
            vaspkit_command_folder = find_folder_in_ancestors('vaspkit_command')
            if vaspkit_command_folder:
                shutil.copy(vaspkit_command_folder / 'KPOINTS', new_folder_path / 'KPOINTS')

            # Run vaspkit command
            os.chdir(new_folder_path)
            os.system('(echo 103)|vaspkit')
            os.chdir('..')

if __name__ == '__main__':
    main()

这个Python脚本执行了以下任务：

遍历当前工作目录，寻找所有以“POSCAR”开头并以“.txt”结尾的文件。

对于找到的每个这样的文件，读取其第一行内容，这一行内容被用作新建文件夹的名称。

在当前目录中创建一个新的文件夹，其名称从步骤2中获得。

将原始的“.txt”文件重命名为“POSCAR”（无后缀名），并将其移动到步骤3中创建的新文件夹内。

查找名为“incar”的文件夹，该文件夹位于当前目录、上一级目录、上上一级目录或上上上一级目录中。

从找到的“incar”文件夹中复制“INCAR_structural_relaxation”文件到新建文件夹，并将其重命名为“INCAR”。INCAR_structural_relaxation是粗收敛的参数。

查找名为“vaspkit_command”的文件夹，该文件夹位于与“incar”文件夹相同的层级目录中。

从找到的“vaspkit_command”文件夹中复制“KPOINTS”文件到新建文件夹，并保持名称不变。

在新建文件夹内运行“vaspkit”的103功能，准备POTCAR。

完成上述操作后，脚本将返回到原始的工作目录。

脚本的目的是自动化创建VASP计算所需文件夹的过程，并准备好进行计算的所有必要文件。这可以大大简化设置多个VASP计算的工作流程。

运行脚本后得到多个计算文件，可以直接运行计算。

2.2.2 auto_vibm.py

在结构优化计算之后，需要进行频率计算以进行零点能校正和热力学计算。其中需要固定基底原子，放开吸附原子。这个命令对于不同路径不一样，于是将其写在vaspkit_command文件夹的fixatom中。可以按照不同体系更改

这里的fixatom


(VASP) [zhaoxb@ibcln01 vaspkit_command]$ cat fixatom 
O 1 24
Ag 1 48

表明固定前24个O 和前48个Ag 把1改成-1则是固定自后向前的原子

在结构优化后的文件夹内运行auto_vibm.py


#!/usr/bin/env python3
# -*- coding: UTF-8 -*-
import time
import os
import shutil
from pathlib import Path
import subprocess


# 查找 vaspkit_command 文件夹
def find_folder(folder_name):
    current_path = Path(os.getcwd())
    for parent in [current_path.parent, current_path.parent.parent, current_path.parent.parent.parent]:
        folder_path = parent / folder_name
        if folder_path.is_dir():
            return folder_path
    return None


# 执行 auto_fixatom.py 脚本来固定原子
def run_auto_fixatom():
    try:
        subprocess.run(['auto_fixatom.py'], check=True)
        print("auto_fixatom.py ran successfully.")
    except subprocess.CalledProcessError as e:
        print(f"An error occurred while running auto_fixatom.py: {e}")



# 复制文件并提供反馈
def copy_file(src, dst):
    try:
        shutil.copy(src, dst)
        print(f"Copied {src} to {dst}")
    except Exception as e:
        print(f"Error copying {src} to {dst}: {e}")


# 主函数
def main():
    # 首先运行 auto_fixatom.py 脚本来创建 CONTCAR_FIX.vasp
    run_auto_fixatom()

    # 查找必要的文件夹
    vaspkit_command_path = find_folder('vaspkit_command')
    incar_path = find_folder('incar')

    if not vaspkit_command_path:
        print("Error: 'vaspkit_command' folder not found.")
    if not incar_path:
        print("Error: 'incar' folder not found.")

    # 创建 fc 文件夹并复制文件
    cdir_str = os.getcwd()
    cur_filename = Path(cdir_str).name
    vib_dir = Path(cdir_str) / (cur_filename + "fc")
    vib_dir.mkdir(exist_ok=True)

    copy_file(Path(cdir_str) / 'POTCAR', vib_dir / 'POTCAR')
    copy_file(Path(cdir_str) / 'KPOINTS', vib_dir / 'KPOINTS')
    if incar_path:
        copy_file(incar_path / 'INCAR_frequency_calculation', vib_dir / 'INCAR')

    # 现在使用auto_fixatom.py脚本生成的CONTCAR_FIX.vasp
    contcar_fix_path = Path(cdir_str) / 'CONTCAR_FIX.vasp'
    if contcar_fix_path.is_file():
        shutil.copy(contcar_fix_path, vib_dir / 'POSCAR')
        print(f"Copied {contcar_fix_path} to {vib_dir / 'POSCAR'}")
    else:
        print(f"Error: {contcar_fix_path} does not exist.")

    print(f"已建立文件夹及输入文件 {vib_dir}")


if __name__ == "__main__":
    main()

其中auto_fixatom.py为


#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from pymatgen.io.vasp import Poscar
from pymatgen.core.structure import Structure
from pathlib import Path

# Function to find 'fixatom' file in current and up to three parent directories
def find_fixatom_file():
    current_path = Path.cwd()
    for i in range(4):  # Check current and three parent levels
        vaspkit_command_path = current_path / 'vaspkit_command'
        fixatom_path = vaspkit_command_path / 'fixatom'
        if fixatom_path.is_file():
            return fixatom_path
        current_path = current_path.parent
    return None

# Function to parse the fixatom file and get the list of fixed atoms
def parse_fixatom_file(fixatom_path, structure):
    fixed_atoms_info = []
    with open(fixatom_path, 'r') as file:
        for line in file:
            parts = line.split()
            if len(parts) == 3:
                element, start, count = parts
                start, count = int(start), int(count)
                if start > 0:  # Positive start index
                    indices_to_fix = list(range(start - 1, start - 1 + count))
                elif start < 0:  # Negative start index means counting from the end
                    indices_to_fix = list(range(-count, 0))
                fixed_atoms_info.append((element, indices_to_fix))
    return fixed_atoms_info

# Function to fix atoms based on fixatom file
def fix_atoms(structure, fixed_atoms_info):
    # Create a mapping from element to its indices
    element_indices_map = {}
    for idx, site in enumerate(structure):
        element = site.species_string
        if element not in element_indices_map:
            element_indices_map[element] = []
        element_indices_map[element].append(idx)

    selective_dynamics = [[True, True, True]] * len(structure)  # Default to all movable

    for element, indices_to_fix in fixed_atoms_info:
        element_indices = element_indices_map[element]
        if indices_to_fix[0] >= 0:  # Positive indices are directly used
            for idx in indices_to_fix:
                selective_dynamics[element_indices[idx]] = [False, False, False]
        else:  # Negative indices mean counting from the end
            for idx in indices_to_fix:
                selective_dynamics[element_indices[idx]] = [False, False, False]

    # Update the selective dynamics property of the structure
    structure.add_site_property('selective_dynamics', selective_dynamics)
    return structure

# Main function
def main():
    fixatom_path = find_fixatom_file()
    if not fixatom_path:
        print("Error: 'fixatom' file not found.")
        return

    structure = Poscar.from_file("CONTCAR").structure
    fixed_atoms_info = parse_fixatom_file(fixatom_path, structure)
    fixed_structure = fix_atoms(structure, fixed_atoms_info)

    new_structure_filename = "CONTCAR_FIX.vasp"
    Poscar(fixed_structure).write_file(new_structure_filename, vasp4_compatible=False)
    print(f"Fixed atoms as per 'fixatom' and saved new structure to {new_structure_filename}")

if __name__ == "__main__":
    main()

这个脚本为频率计算创建必要的文件和文件夹。脚本的操作步骤如下：

查找必要的文件夹：函数 find_folder 在当前目录的父级、祖父级和曾祖父级目录中查找名为 vaspkit_command 和 incar 的文件夹。

读取fixatom 文件并固定原子：调用auto_fixatom.py，得到固定了原子的CONTCAR_FIX.vasp文件

创建频率计算文件夹：脚本会在当前目录下创建一个名为 当前文件夹名fc 的新文件夹，用于存放频率计算相关的文件。

复制文件到新文件夹：

脚本会从当前目录复制 POTCAR 和 KPOINTS 文件到新文件夹。

如果找到 incar 文件夹，它会复制 INCAR_frequency_calculation 文件到新文件夹，并重命名为 INCAR。

同样，它还会复制 CONTCAR_FIX.vasp 文件到新文件夹，并重命名为 POSCAR。

脚本最后打印出创建的文件夹和输入文件的信息，提供给用户确认。脚本通过自动化这些重复性的步骤，帮助用户节省时间，减少手动设置VASP计算的错误。

2.3 进行自动化工作流的脚本

2.3.1 提交任务脚本svasp与job.vasp

首先介绍在THU的探索1000上提交任务的脚本：

目前我在服务器上使用的提交任务的命令是svasp，这个脚本的内容是：


#!/bin/bash

current_dir=$(basename $(pwd))
sbatch --job-name=$current_dir job.vasp

其中job.vasp的内容是：


#!/bin/bash
#SBATCH -p cnall
#SBATCH -N 1
#SBATCH -o stdout.%j
#SBATCH -e stderr.%j
#SBATCH --no-requeue
#SBATCH --ntasks-per-node=56

module load compilers/intel/oneapi-2023/config
module load soft/vasp/vasp.6.3.2
mpirun -np 56 vasp_std

2.3.2 实现自动化的脚本

这部分脚本有三个

subvasp


#!/bin/bash

# 第一阶段：提交VASP任务
svasp

# 在后台启动收敛检查脚本
nohup bash check_convergence.sh &

check_convergence.sh


#!/bin/bash

# 等待20秒
sleep 20
# 循环检查第一次任务是否收敛
while true; do
    if grep -q 'reached required accuracy - stopping structural energy minimisation' OUTCAR; then
        echo "First phase converged."
        break
    fi
    sleep 60 # 每60秒检查一次
done

# 将CONTCAR复制为新的POSCAR
cp CONTCAR POSCAR

# 调用Python脚本处理INCAR
handle_INCAR_structural_relaxation_fine.py

# 提交第二次VASP任务
svasp

# 等待20秒
sleep 20

# 循环检查第二次任务是否收敛
while true; do
    if grep -q 'reached required accuracy - stopping structural energy minimisation' OUTCAR; then
        echo "Second phase converged."
        break
    fi
    sleep 60
done

# 第二次收敛后，运行 vibm.py 创建频率计算文件夹
vibm.py

# 获取新建立的频率计算文件夹的名称
cdir_str=$(pwd)
cur_filename=$(basename "$cdir_str")
vib_dir="${cdir_str}/${cur_filename}fc"

# 进入频率计算文件夹并提交计算任务
if [ -d "$vib_dir" ]; then
    cd "$vib_dir"
    # 提交频率计算任务
    svasp
    echo "Frequency calculation task submitted in $vib_dir"
else
    echo "Error: Frequency calculation folder not found."
fi

handle_INCAR_structural_relaxation_fine.py


#!/usr/bin/env python3

import os
import shutil

def find_incar_file():
    # 检查各级目录
    for i in range(4):
        path = "../" * i
        incar_path = os.path.join(path, 'incar', 'INCAR_structural_relaxation_fine')
        if os.path.isfile(incar_path):
            return incar_path
    return None

def main():
    incar_file = find_incar_file()
    if incar_file:
        shutil.copy(incar_file, 'INCAR')
        print(f"Copied INCAR file from {incar_file}")
    else:
        print("Error: INCAR file not found.")

if __name__ == "__main__":
    main()

使用时只需要在计算文件夹使用subvasp提交任务即可。

这几个脚本组成了一个自动化流程，用于提交VASP（Vienna Ab initio Simulation Package）计算任务，并根据计算结果自动启动后续的计算步骤。下面是它们各自的功能介绍：