# MgO surface cluster geometry optimisation

# Initialise ChemShell
from chemsh import *
from chemsh.utils import testutils

# Import initial MgO cluster geometry
mgo = Fragment(coords='mgo_shells_regions.pun', connect_mode=None)

# Get list of QM atoms and charge from input fragment
qm_region = mgo.getRegion(1)
print("QM region is: ", qm_region)

qm_charge = round(mgo.totalcharge)
print("QM charge is: ", qm_charge)

# Get list of active atoms for optimisation
active_region = mgo.getRegion(1,2,3)
print("Active region is: ", active_region)

# QM/MM settings
qm_theory = NWChem(method='dft',
                   functional='b3lyp',
                   basis='mgo_nwchem.basis',
                   ecp='mgo_nwchem.ecp',
                   charge=qm_charge,
		   path="mpirun -np 4 nwchem")

mm_theory = GULP(ff='mgo_2body.ff',
                 conjugate=True)

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

# Run QM/MM optimisation
opt = Opt(theory=qmmm, 
          algorithm="lbfgs", trust_radius="energy",
          maxcycle=200, maxene=300, active=active_region)
opt.run()

ecalc = opt.result.energy

# Save optimised geometry
mgo.save('mgo_shells_opt.pun', 'pun')
mgo.save('mgo_shells_opt.xyz', 'xyz')

# Validate result against pre-calculated reference energy
eref = -1083.536372673275
tolerance = 1.0e-6

testutils.validate(ecalc, eref, tol=tolerance)
