#
# (MgO)25 nanocluster optimisation using separate boundary ECPs
# NWChem/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_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 = NWChem(method='dft',
                   functional='blyp',
                   basis='mgo.basis',
                   ecp='mgo_sepecp.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=5)

#
# Run QM/MM optimisation
#
opt = Opt(theory=qmmm,
          algorithm="lbfgs", tolerance=0.0003,
          maxcycle=400, maxene=300)

opt.run()

mgo.save('mgo25_sepecp_opt.pun', 'pun')

ecalc = opt.result.energy

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

testutils.validate(ecalc, eref)
