#
# (MgO)25 nanocluster optimisation 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=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 = -310.166389367120 
# AJL, 03/02/2024. Updated as believed outdated due to FHI-aims version updates
eref = -310.166388197651

testutils.validate(ecalc, eref)
