'''
This sample is only for demonstration purposes and may not be chemically meaningful
'''

from chemsh       import *
from chemsh.utils import testutils

# read in the partitioned fragment from punch file
mgo_cluster = Fragment(coords='../../02b_partitioning_explicit/mgo_shells_qmmm_explicit_larger_qm.pun', connect_mode=None)

# use the total charge for QM
qm_charge = int(mgo_cluster.totalcharge)

# separate_ecps=True is mandatory for NWChem (already chosen by NWChem default)
qm_theory = NWChem(method='dft',
                   functional='blyp',
                   # use external filename
                   basis='mgo_nwchem.bas',
                   # you'd better well name the filename
                   ecp='mgo_nwchem_explicit_separate_fragmovecs.ecp',
                   charge=qm_charge,
                   fragmovecs='elements',
                   maxiter=500,
                   tol_pre_hf=1E-2,
                   threshold=1E-2,
                   tol_energy=1E-2,
                  )

mm_theory = GULP(ff='mgo_2body.ff',
                 # conjugate gradients
                 conjugate=True)

# NWChem only supports explicit interface so we choose region "1" as the QM region
qm_region = mgo_cluster.getRegion(1)
print(" qm_region =", qm_region)
# alternatively
qm_region = mgo_cluster.getRegion("1")
print(" qm_region =", qm_region)

# fragment is assigned to QMMM instance
qmmm = QMMM(qm=qm_theory,
            mm=mm_theory,
            frag=mgo_cluster,
            qm_region=qm_region,
            # choose coupling='ionic' to switch off link atom and charge shift (those required by coupling='covalent')
            coupling='ionic',
            # normally shell relaxation converges robustly within 20 cycles
            shl_maxcycles=20,
            # default: 0.0001 although higher criterion can be reached
            shl_tolerance=0.01)

# create a single-point task instance and run
# alternative: sp = SP(theory=qmmm).run()
sp = SP(theory=qmmm)
sp.run(dryrun=False)

print(__doc__)
