本文首发于朋友的公众号计算运维鸟，如需转载请注明出处。
进行一个过渡态计算，使用NEB方法时，初始路径插值方式好不好会很大程度影响计算效率。恰好今天科音论坛上有坛友问到如何使用idpp的方式进行初始路径插值，也提到李湃和ponychen的脚本，但使用时遇到了些许问题。虽然问题部分被解决了，但仔细考虑后脑袋里蹦出来这几个问题：- POSCAR中的原子顺序不匹配或元素名称被省略是否会报错？
- 是否能支持更多格式，例如Materials Studio导出的xtd或arc？
- idpp的选项能否有更多的选择，如极小化算法可选择BFGS或CG？
- 是否能够在Windows甚至是Mac上也能完成产生和预览？
- 有没有现成的玩意能一次性解决上面几个问题并基本上再也不用改脚本？
本想以前四个问题为导向将原来的脚本进行修改，但后来发现通过ase直接就能够解决，根本没必要重新造这个轮子，遂分享一下程序设计的经验，如果只是想使用这个脚本的伙伴可以直接看本文最后程序使用的方式。这个python程序我们暂时起名为makeneb.py，接下来设计思路可以分为三步走：可以假想这个程序输入参数的方式，例如用makeneb.py -i POSCAR_ini POSCAR_fin -n 5来读入初始结构和最终结构并指定插入3帧结构，那么我们就可以借助python的argparse模块来完成这个功能。示例代码如下：# 读取数据import argparsep = argparse.ArgumentParser()# 用nargs来让-i参数获取1个或多个结构的文件名。p.add_argument('-i', '--images', nargs='+', type=str, help='Images defining path from initial to final state.')# 默认为6，因为大多服务器核心数是4的倍数。p.add_argument('-n', '--nimage', type=int, default=6, help='Number of images in a band. Default: 6')
最终将包含运行参数的对象p传回，如果再将p中的参数解析出来就可以得到这些参数。这里使用argparse只是比较适合初学者了解，并非唯一方案，若要实现更好的体验可以使用PyQT或curses来实现更复杂的用户交互。接下来我们希望程序能读取POSCAR的格式，恰好ase的io模块能够胜任这个问题，甚至还能读取一段轨迹，例如：# 通过ase.io读取始末结构。from ase import io# 读取一个结构，返回的是ase.Atoms的数据类型。initial = io.read('POSCAR_ini')final = io.read('POSCAR_fin')# 读取一个轨迹，注意返回的是生成器，使用list可以被转换为ase.Atoms的数组。images = io.iread('XDATCAR')
接下来就可以使用ase.neb模块来完成插点的工作了，neb模块的相关参数说明可以参考https://wiki.fysik.dtu.dk/ase/ase/neb.html的说法，本文不再赘述，仅需要注意进行NEB插点前需要先初始化一个与NEB需要计算的帧数相同的数组，将初始结构和最终结构分别放在数组首尾，当中可以填入初始结构或最终结构，之后转换为ase.neb.NEB的类型并进行插点：# 初始化路径，共nimage个结构。images = [initial]+[initial.copy() for i in range(nimage-2)]+[final]from ase.neb import NEB# 转换为ase.neb.NEB类型，并且设置弹簧常数为0.1。neb = NEB(images, k=0.1)# 对于idpp插点方法设置插点参数，使用BFGS作为优化算法，力收敛阈值为0.1，最大步数为100。from ase.optimize import BFGSneb.idpp_interpolate( fmax=0.1, optimizer=BFGS, steps=100)
此时neb中的images发生更新，我们将这个迭代完成后得到的images再按照需要的方式写入即可。
这里以加入VTST的VASP为例子，需要建立与NEB路径上点数相同的文件夹，例如6帧就要建立00, 01, ... , 05六个文件夹，然后将路径上各个结构的POSCAR格式分别拷贝到这些文件夹下。通过python的os模块直接就可以实现上述功能：
import os# 建立00, 01, ..., 0X一系列文件夹，如果这些文件夹存在就跳过。[os.mkdir('%02d' % i)for i in range(len(images)) if not os.path.exists('%02d' % i)]# 到对应目录去写入POSCAR文件。from ase import io[io.write('%02d/POSCAR' % i, images[i]) for i in range(len(images))]# 如果需要将images输出为XDATCAR，可以这样。io.write('XDATCAR', images)
最终将上述的三部分进行组合适当调整，构造出如下完整程序(https://github.com/DYang90/ase_extend/blob/main/makeneb.py)。程序中内置了简单的帮助，可以通过：python3 makeneb.py --help
要使用上述程序产生《Learn VASP The Hard Way》中提到的乙烷旋转的例子，此处以使用Matertials Studio作为辅助为例，步骤如下：- 构建一个乙烷结构，保存成两份，命名为r.xsd，另一份命名为p.xsd。
- 在菜单Tools->Reaction Preview选项卡，将反应物设置为r.xsd，产物设置为p.xsd，Number of frames设置为2，将乙烷中的氢原子按旋转的情况手工进行匹配，并点击Preview按钮，此时会产生r-p.xtd。
- 上述r-p.xtd文件实际上只保留了显示样式，实际的轨迹是被保存在相同文件夹下一个对应的r-p.arc的隐藏文件，简单起见我们可以将r-p.arc先转换为r-p.xyz轨迹：
ase convert -i dmol-arc r-p.xyz
需要注意的是，此处xyz其实是扩展的xyz，允许保存晶格信息。
- 为了做对照，我们分别采用线性插点和IDPP两种方式来得到轨迹，
python3 makeneb.py -i r-p.xyz -n 10 -o
python3 makeneb.py -i r-p.xyz -n 10 --method idpp --optimizer BFGS -o
最终得到两组结果，通过ase+Povray进行渲染，预览效果如上图所示，可见这样的体系线性插点确实有严重的缺陷，也同时也体现出ase强大之处，大家如果对程序有什么使用问题可以在本文下方留言。