# NEB calculation for the migration of the electron hole from oxygen 7 to oxygen 5

# Initialise ChemShell
from chemsh import *   

# Imports the previously optimised localised-cluster model
rs = Fragment(coords="O7-localised.pun") # reactant structure
ps = Fragment(coords="O5-localised.pun") # product structure

# Extracts the qm 'region' and 'charge' parameters and stores them as variables
qm_region = rs.getRegion(1)
print("The QM region is ", qm_region)
qm_charge = round(rs.totalcharge)
print("The QM charge is ", qm_charge)
active_region = rs.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=rs, qm_region=qm_region, 
           coupling='ionic', shl_maxcycles=5)

# Run an NEB optimisation between the optimised endpoints where the holes are localised on different oxygens
opt_path = Opt(theory=qmmm, frag2=ps, neb="frozen", nimages=10,
               coordinates="cartesian", nebk=0.01, trust_radius="const", 
               active=active_region)
opt_path.run()

