# Transition state vibrational analysis calculation using the constructed Li embedded solid-state MgO cluster

# Initialise ChemShell
from chemsh import *   

# Imports the previously optimised localised-cluster model
guess = Fragment(coords="cimg.pun") 

# Extracts the qm 'region' and 'charge' parameters and stores them as variables
qm_region = guess.getRegion(1)
print("The QM region is ", qm_region)
qm_charge = round(guess.totalcharge)
print("The QM charge is ", qm_charge)
active_region = guess.getRegion(1,2,3)
print("The active region is ", active_region)

# Defines QM and MM theory settings
qm_theory = NWChem(method='dft', functional='pbe0', 
                   basis='li-mgo_nwchem.basis', ecp='li-mgo_nwchem.ecp', 
                   charge=qm_charge, path='mpirun -np 4 nwchem', maxiter=500)
 
mm_theory = GULP(ff='mgo_2body.ff', conjugate=True)

qmmm= QMMM(qm=qm_theory, mm=mm_theory, 
           frag=guess, qm_region=qm_region, 
           coupling='ionic', shl_maxcycles=5)

# Run a vibrational analysis on the guess of the transition state
opt_path = Opt(theory=qmmm, coordinates="cartesian", 
               thermal=True, active=active_region)
opt_path.run()
