#
# (MgO)25 nanocluster energy evaluation using separate boundary ECPs
# FHI-aims/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.
#
# Note the result differs substantially from GAMESS-UK and NWChem
# versions because of all-electron Mg1.

from chemsh import *
from chemsh.utils import testutils

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

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

#
# QM/MM settings
#
qm_theory = FHIaims(
                    #settings='fhiaims.settings.mgo', 
                    basis='fhiaims.basis.mgo',
                    ghost_species='Emptium',
                    ghost_min_cutoff=0.4, #This is for accurate validation against TCL ChemShell. More accurate default is 0.001 Angstrom
                    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 = -309.791089376992 

testutils.validate(ecalc, eref)
