# GAMESS-UK C2H2 tautomerisation DFT (BLYP) reaction path optimisation
# Step 2: nudged elastic band calculation

# Initialise ChemShell
from chemsh import *

# Read in optimised reactant and product state geometries from disk
rs = Fragment(coords="reactant_opt.pun")
ps = Fragment(coords="product_opt.pun")

# Set up the QM calculation
neb_theory = GAMESS_UK(frag=rs, method="dft", functional="blyp", basis="3-21g")

# Run an NEB optimisation between the optimised endpoints
optPath = Opt(theory=neb_theory, frag2=ps, neb="frozen", nimages=8,
              coordinates="cartesian", nebk=0.01, 
              maxstep=0.9, trust_radius="const")
optPath.run()
