# Single point MgO surface cluster energy evaluation

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

# Import MgO cluster geometry from previous stage
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)

# 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=20)

# Calculate the QM/MM energy
sp = SP(theory=qmmm, gradients=False)
sp.run()

ecalc = sp.result.energy

# Validate result against pre-calculated reference energy
eref = -1083.508825142250
tolerance = 1.0e-7

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