Nudged Elastic Band (NEB) 法で最小エネルギー経路を探索しよう

pip3 install --upgrade --user ase

#!/usr/bin/python3

# ===== 設定 =====

# 始状態の構造（xyz 形式）
geom_of_reactant = "open.xyz"
# 終状態の構造（xyz 形式）
geom_of_product  = "closed.xyz"
# NEB 計算の離散点（images）の個数
numb_of_nodes = 18
# 計算の収束条件
convergence = 0.4

# 計算方法
gaussian_method = "B3LYP"
# 基底関数
gaussian_basis = "6-31G(d)"
# 電荷
gaussian_charge = 0
# スピン多重度
gaussian_mult = 1
 
# ===== 必要な設定はここまで =====

import os

from ase import Atoms
from ase.io import read,write
from ase.calculators.gaussian import Gaussian
from ase.build.rotate import minimize_rotation_and_translation
from ase.optimize import BFGS
from ase.neb import NEB

atoms_reac = read(geom_of_reactant)
atoms_prod = read(geom_of_product)

if not os.path.exists("neb"):
    os.mkdir("neb")

images = []
for n in range(numb_of_nodes):
    if n < (numb_of_nodes // 2):
        image = atoms_reac.copy()
    else:
        image = atoms_prod.copy()
    minimize_rotation_and_translation(atoms_reac, image)
    image.calc = Gaussian(label = 'neb/Gaussian', mem = "19200MB", chk = "Gaussian.chk", save = None,
                          method = gaussian_method, basis = gaussian_basis, charge = gaussian_charge, mult = gaussian_mult)
    images.append(image)

neb = NEB(images)
neb.interpolate()

pote_old = [image.get_potential_energy() for image in images]

path_opt = BFGS(neb, trajectory = 'neb.traj', logfile = 'neb.log')
path_opt.run(fmax = convergence)

pote_new = [image.get_potential_energy() for image in images]

with open("neb.dat", 'w') as fw:
    for n in range(numb_of_nodes):
        fw.write(f"{n:2d}  {pote_old[n]:20.10f}  {pote_new[n]:20.10f}\n")

if os.path.exists("neb.xyz"):
    os.remove("neb.xyz")
for n in range(numb_of_nodes):
    write("neb.xyz", images[n], append = True)

#!/bin/sh
#PBS -l select=1:ncpus=16:mpiprocs=1:ompthreads=16
#PBS -l walltime=2:00:00

NCPUS=16
FILE="neb.py"

if [ ! -z "${PBS_O_WORKDIR}" ]; then
  cd "${PBS_O_WORKDIR}"
  WORK=/lwork/users/${USER}/${PBS_JOBID}/gaussian
else
  WORK=/gwork/users/${USER}/tmp.$$
fi

if [ ! -d ${WORK} ]; then
  mkdir ${WORK}
fi

. /apl/gaussian/16c02/g16/bsd/g16.profile
export LANG=C
export GAUSS_SCRDIR=${WORK}

export GAUSS_CDEF=`/apl/gaussian/16c02/rccs/cpu.pl -n $NCPUS`

python ${FILE}

if [ -d ${WORK} ]; then
  /bin/rm -rf ${WORK}
fi
exit 0