vasp_neb_idpp-使用方法

[TOC]

Nudged Elastic Band（NEB）介绍

搬运链接：https://mp.weixin.qq.com/s/pQ-GNPx94W0DndkIoRZ74A

Nudged Elastic Band（NEB）是一种常用的用于搜寻过渡态和反应路径的(chain-of-states)方法。将VTST的代码集成到VASP的源代码后重新编译（http://hmli.ustc.edu.cn/doc/app/vasp.5.4.1-vtst.htm），就可以做`CI-NEB（Climbing Image）方法, 相比于常规的NEB方法，它能够更好地定位过渡态。同时VTST`提供了大量的脚本用于辅助过渡态的计算。（https://www.bigbrosci.com/2018/11/10/ex72/）

chain-of-states这类方法主要好处是只需要提供反应物和产物结构就能得到准确的反应路径和过渡态。首先在二者结构之间以类似LST的方式线性、均匀地插入一批新的结构（使用内坐标更为适宜），一般为5~40个，每个结构就是势能面上的一个点(称为image)，并将相邻的点以某种势函数相连，这样它们在势能面上就如同组成了一条链子。对这些点在某些限制条件下优化后，在势能面上的分布描述的就是MEP，能量最高的结构就是近似的过渡态位置。

CI-NEB与NEB的关键区别是能量最高的点受力的定义，在CI-NEB中这个点不会受到相邻点的弹簧力，避免位置被拉离过渡态，而且将此点平行于路径方向的势能力分量的符号反转，促使此点沿着路径往能量升高的方向上爬到过渡态。这个方法只需要很少的点，比如包含初、末态总共5个甚至3个点就能准确定位过渡态，是最有效率的寻找过渡态的方法之一。

（http://sobereva.com/44）

线性插点和IDPP插点方法

如何生成合理的images呢？VTST提供了一种线性插值的生成初始结构的方法(nebmake.pl)。简单来说就是将初态和末态的原子坐标的差值均分来获得一系列处于初态和末态之间的images。这种线性插值的方法对于近乎于直线的原子扩散问题的研究比较适用。如果反应路径并非直线，通过线性插值的点可能距离真正的能量最小路径甚远，此时可能就需要花费更高的时间用来优化。更严重的是涉及到原子旋转的反应路径，线性插值得到的images很可能发生原子重叠的情况，此时优化任务就进行不下去了。

因此Hannes Jónsson等人在2014年发表了一篇文章：Improved initial guess for minimum energy path calculations专门研究了非线性插点的问题。文章里面讲到一种名为image dependent pair potential (IDPP)的方法来简单优化初始的线性插值的结构。其主要思想是利用路径上每个离散点（image）处的成对距离插值和Nudged Elastic Band构造目标函数曲面（IDPP），然后在新构造的势能面上找到最佳的路径。一旦构造了势能面，即可通过现有的优化算法获得最小能量路径。博客中的代码没有考虑在images之间加弹簧力限制，因此可能会出现原子飞走的情况。（http://blog.sina.com.cn/s/blog_b364ab230102wj0g.html）

幸运的是ASE和pymatgen-diffusion这两个Python第三方库已经实现了IDPP插点方法。测试发现ASE的IDPP插点方法不完善，原子可能会出现错位的情况。因此本人基于pymatgen-diffusion和pymatgen两个库开发了一个python脚本idpp.py，实现了与nebmake.pl一样的插点用法，因此熟悉nebmake.pl的用户可以很快熟悉IDPP的插点方法。

idpp.py所需环境安装

idpp.py使用到了pymatgen-diffusion和pymatgen两个Python库，因此需要安装这两个库。使用pip工具可以很快地进行安装。

先更新pip工具
1
python -m pip install --upgrade pip

安装pymatgen和pymatgen-diffusion

pip3 install pymatgen pymatgen-diffusion --user 
#or pip install pymatgen pymatgen-diffusion --user

pip3 --default-timeout=1000 install -U pymatgen pymatgen-diffusion

安装好之后测试pymatgen和pymatgen-diffusion有没有安装好：

1 2	import pymatgen import pymatgen_diffusion

在优化好初态（./ini）和末态结构（./fin）后，使用idpp.py进行插点,在初态和末态之间插6个images:

1	python idpp.py ./ini/CONTCAR ./fin/CONTCAR 6

短暂等待之后，就会得到00~07一个8个文件夹。接下来就设置INCAR，POTCAR和KPOINTS运行NEB任务叭。

下图比较了使用线性插值方法和idpp.py两种插点方法插点的路径的区别。插点后使用vaspkit 504功能将images合并后生成PDB轨迹文件，并使用VMD软件可视化。可以看到使用线性插值的方法得到的路径存在有原子过近的问题，而idpp方法得到了更合理的插点结果。

线性插值方法得到的路径

IDPP方法得到的路径

#!/public/software/apps/Anaconda3/bin/python

import numpy as np
import os
import sys
import warnings
from pymatgen.core import Structure
from pymatgen_diffusion.neb.pathfinder import IDPPSolver

def warn(*args, **kwargs):
    pass
warnings.warn = warn
sys.stdout = open(os.devnull, 'w')

if len(sys.argv) <4:
    raise SystemError('Sytax Error! Run as python idpp ini/POSCAR fin/POSCAR 4')


warnings.filterwarnings("ignore")
init_struct = Structure.from_file(sys.argv[1], False)
final_struct = Structure.from_file(sys.argv[2], False)

obj = IDPPSolver.from_endpoints(endpoints=[init_struct, final_struct], nimages=int(sys.argv[3]),
                                sort_tol=1.0)
new_path = obj.run(maxiter=5000, tol=1e-5, gtol=1e-3,step_size=0.05,\
                         max_disp=0.05, spring_const=5.0)


for i in range(len(new_path)):
    image_file='{0:02d}'.format(i)
    if not os.path.exists(image_file):
        os.makedirs(image_file)
    POSCAR_file=image_file+'/POSCAR'
    new_path[i].to(fmt="poscar", filename=POSCAR_file)

sys.stdout = sys.__stdout__
#Image Dependent Pair Potential for improved interpolation of NEB initial guess
#Reference: S. Smidstrup, A. Pedersen, K. Stokbro and H. Jonsson, Improved initial guess for minimum energy path calculations, J. Chem. Phys. 140, 214106 (2014).
print("Improved interpolation of NEB initial guess has been generated. BYE.")