#
# (MgO)25 nanocluster energy evaluation using the implicit interface embedding model
# GAMESS-UK/GULP version
#
# The starting geometry was cut from the Tcl-ChemShell mgo_embed tutorial slab
# to match structure 25a in M. Chen et al, JPCA 2014, 188, 3136.
# Force field, basis set and ECP are from the mgo_embed tutorial.
#

from chemsh import *
from chemsh.utils import testutils

#
# Delete any connectivity to ensure link atoms are not used
#
mgo = Fragment(coords='mgo25_implicit.pun', connect_mode=None)

#
# Get list of QM atoms and charge from input fragment
#
qm_charge = round(mgo.totalcharge)
qm_region = mgo.getRegion(1,2)

#
# QM/MM settings
#
qm_theory = GAMESS_UK(method='dft',
                   functional='blyp',
                   basis='mgo.basis',
                   ecp='mgo_implicit.ecp',
                   charge=qm_charge)

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)

#
# Run single point energy evaluation
#
sp = SP(theory=qmmm, gradients=True)

sp.run()

ecalc = sp.result.energy

# Test result from Tcl-ChemShell
eref = -110.47325105
tolerance = 1.0e-7

testutils.validate(ecalc, eref)
