NEB方法寻找过渡态（Gaussian与ASE联用）
| Created | |
|---|---|
| Tags | ASEGaussian公众号投稿 |
寻找过渡态是个棘手的工作，如果尝尽了各种办法（软扫描、换计算精度等）都找不到过渡态，那么不妨尝试下NEB方法。
一串弹簧连接珠子放到马鞍上，那么珠子串在重力与之间弹簧拉力共同作用下会分布在马鞍的“最低能量路径”上。NEB方法利用类似原理，生成一串连接反应物与产物的珠子，经过优化，珠子构成的路径就是势能面上的能量最低路径（MEP），能量最高点与过渡态结构相近。
使用NEB方法的应用场景是，用户已经找到反应物与产物稳定结构，且Gaussian的QST2、QST3尝试失败。
ASE简介与安装
ASE 全程 Atomic Simulation Environment，是一个操作原子的python环境。它被设计成一个平台，提供了调用多种量化计算软件的接口，通过调用这些量化软件得到需要的信息，从而完成稍微复杂的组合计算。它提供一些算法优化结构、找过渡态、寻找全局最小值点，也能用进行MD模拟或者NEB计算。ASE用户大多来自材料模拟领域。使用该环境需要用户有一定的python编程基础，因为用户需要编写一个脚本完成一套组合计算。
本文讲解Gaussian与ASE联用进行NEB计算。ASE的安装非常简单，只需要用户执行，
pip install --upgrade --user ase
如果没有网络，用户可以尝试在有网络的环境下将ASE与Gaussian打包成docker镜像，传到服务器上。
运行python，如果import ase成功，则说明ASE安装成功。另外机器里需要安装了Gaussian。
计算脚本的编写
先把需要的类与方法导入
from ase import Atoms,io
from ase.io import iread,read
from ase.build import minimize_rotation_and_translation
from ase.optimize import BFGS,MDMin
from ase.neb import NEB
from ase.parallel import world
from ase.calculators.gaussian import Gaussian
import os
import sys1. 输入反应物、产物结构
将反应物、产物存为xyz格式。本文中反应物为A.xyz，产物是B.xyz。之后削减两个结构之间的RMSD[cite Kerbash算法]。
reac = read('A.xyz') # Reactant saved in A.xyz
prod = read('B.xyz') # Product saved in B.xyz
minimize_rotation_and_translation(reac,prod) # Reduce RMSD between reac and prod这样反应物、产物便被实例化为reac和prod，他们属于ASE中的Atoms类，可以理解为一个分子。
2. 定义与gaussian的接口
与gaussian的接口定义在calculator类中。我们通过以下方法定义calculator，这些calculator决定了Gaussian完成什么样的计算。
gauss_CALC = Gaussian(label='calc/gaussian', # Calculation file name
nprocshared="32", # %nprocshared in Gaussian
mem='5GB', # %mem in Gaussian
xc='m06l', # functional
basis='def2svp', # basis set
ioplist=["2/9=2000"], # Must add, avoid standard rotation
scf='maxcycle=100') # scf options需要注意的是，Gaussian计算会旋转分子到standard orientation，必须加入ioplist=["2/9=2000"]以防止分子旋转。我们实例化了一个gauss_CALC，这是与gaussian的接口。下文会用到。
3. 计算反应物、产物的能量
反应物、产物在计算中不进行结构优化。因此建议计算下这两个结构的能量。
reac.calc = gauss_CALC # define calculator for reactant
reac.calc.label = "calc/gaussianr" # Rename Gaussian File for reactant
E = reac.get_potential_energy() # calculate reactant energy
print("Single-point energy of reactant: ",E)
prod.calc = gauss_CALC # define calculator for product
prod.calc.label = "calc/gaussianp" # Rename Gaussian File for product
E = prod.get_potential_energy() # calculate product energy
print("Single-point energy of product: ",E)4. 定义NEB
在这里，我们在string上定义21个beads（不含反应物与产物）。指定方法为将首（反应物）与尾（产物）之间进行线性插值。
Nbeads=21
images = [reac]
images += [reac.copy() for i in range(Nbeads)]
images += [prod]
i = 0
for image in images:
i = i+1
image.calc = deepcopy(gauss_CALC)
image.calc.label='calc/gaussian'+str(i)
neb = NEB(images)
neb.interpolate() # Interpolate the images Linearly5. 优化NEB
dt是优化步长，建议比默认值小一些。fmax是收敛阈，当这个NEB上每个beads受力很小时候，我们认为优化结束。
optimizer = MDMin(neb, trajectory='A2B4.traj',dt=0.05)
optimizer.run(fmax=0.05)###neb for multi points interpolate
from ase import Atoms,io
from ase.io import iread,read
from ase.build import minimize_rotation_and_translation
from ase.optimize import BFGS,MDMin
from ase.neb import NEB
from ase.parallel import world
from ase.calculators.gaussian import Gaussian
import os
import sys
from ase.io.trajectory import Trajectory
from ase.io import write
import io, os, copy
os.environ["OMP_NUM_THREADS"] = "8"
os.environ["MKL_NUM_THREADS"] = "8"
os.environ["OMP_STACKSIZE"] = "8G"
E = 0
Nbeads=21
TS_guess_num = 1
## Define Structures
reacStruct = "A.xyz"
prodStruct = "B.xyz"
TSs =["TS1.xyz", "TS2.xyz"]
#QM Calculator
CalcTemplate = Gaussian(label='calc/gaussianr',
nprocshared="2",
mem='5GB',
xc='pbepbe',
basis='sto-6g',
ioplist=["2/9=2000"],
scf='maxcycle=100,xqc')
def outputAtomstoXYZ(Stringatoms,fname):
string = 'structure'
#get current working directory and make a scratch
#directory
path = os.getcwd()
path = path + '/scratch'
if not os.path.exists(path): os.makedirs(path)
#output file name
outFileName = fname
if(os.path.exists(fname)):os.remove(fname)
#write each structure from the .traj file in .xyz format
for img in Stringatoms:
string = 'structure%03d' % (i,) +'.xyz'
outStruct = os.path.join(path, string)
write(outStruct,img)
#combines all optimization structures in one trajectory
#file
inFile = open(os.path.join(path, 'structure%03d' %
(i,) +'.xyz'), 'r')
fileStr = inFile.read()
outFile = open(outFileName, 'a')
outFile.write(fileStr)
## 1. Load Structures and generate string.
### 1.0 Load structures in
reac = read(reacStruct) # Reactant saved in 1.xyz
prod = read(prodStruct) # Product saved in 1.xyz
TS_List = []
if(len(TSs)>0):
for tsName in TSs:
TS_List.append(read(tsName))
reac.calc = copy.deepcopy(CalcTemplate)
reac.calc.label = "R"
prod.calc = copy.deepcopy(CalcTemplate)
prod.calc.label = "P"
E = reac.get_potential_energy()
E = prod.get_potential_energy()
print("Single-point energy of reactant: ",E)
print("Single-point energy of product: ",E)
guessPointOnString = [reac] + TS_List + [prod]
### 1.1 Generate string
nbeadsperPath = int(Nbeads/(1+len(TS_List)))
beadsArrange = [nbeadsperPath for i in range(len(TS_List)+1)]
beadsArrange[-1] = int(Nbeads - nbeadsperPath*(1+len(TS_List)))+beadsArrange[-1]
print("Beads on each paths ",beadsArrange)
imagesList = []
for i in range(len(beadsArrange)):
if(i == len(beadsArrange)-1):
StrPoint = beadsArrange[i]-2
LenAdd = StrPoint+2
else:
StrPoint = beadsArrange[i]-1
LenAdd = StrPoint+1
minimize_rotation_and_translation(guessPointOnString[i],guessPointOnString[i+1])
images = [guessPointOnString[i]]
images += [guessPointOnString[i].copy() for j in range(StrPoint)]
images += [guessPointOnString[i+1]]
neb = NEB(images,climb=True)
neb.interpolate()
images = neb.images
for img in range(LenAdd):
imagesList.append(images[img])
outFileName = 'Guess.xyz'
# Save string of guess
outputAtomstoXYZ(imagesList,outFileName)
print(len(imagesList))
## 2. Optimize:
Tag = 1
for img in imagesList:
img.calc = copy.deepcopy(CalcTemplate)
img.calc.label='calc/gaussian'+str(Tag)
Tag += 1
neb = NEB(imagesList,climb=True)
optimizer = MDMin(neb, trajectory='A2B01.traj',dt=0.5)
optimizer.run(fmax=0.5)
outputAtomstoXYZ(neb.images,"opt05.xyz")
optimizer = MDMin(neb, trajectory='A2B01.traj',dt=0.1)
optimizer.run(fmax=0.1)
outputAtomstoXYZ(neb.images,"opt01.xyz")
optimizer = MDMin(neb, trajectory='A2B01.traj',dt=0.1)
optimizer.run(fmax=0.05)
outputAtomstoXYZ(neb.images,"opt005.xyz")
optimizer = MDMin(neb, trajectory='A2B001.traj',dt=0.05)
optimizer.run(fmax=0.01)
outputAtomstoXYZ(neb.images,"opt001.xyz")