# 
# Water dimer
#

from chemsh       import *
from chemsh.utils import testutils

water = Fragment(coords='../h2o_dimer.cjson')

ff = '''#
# molmec removes Coulomb terms between bonded atoms
# and between atoms bonded to a common third atom
keyword molmec
# Flexible TIP3P model
# last number in LJ term is a cutoff
lennard 12 6 kcal x13
O O 582000.0 595.0 99.9
# intramolecular force-field parameters for MM
# from QM calculation on water monomer
# Note GULP k is twice that of ChemShell format
harmonic bond kcal
H O 1059.162 0.96
# Note GULP k is twice that of ChemShell format
# note different ordering of atoms
three bond kcal
O H H  68.087  104.5
'''


qm = PySCF(method="dft", functional="blyp", basis="3-21g", threshold=1e-9)
mm = GULP(ff=ff, molecule=True)

qmmm = QMMM(frag=water, qm_region=range(3), qm=qm, mm=mm, embedding='electrostatic', coupling='covalent')

sp = Opt(theory=qmmm)

sp.run()

ecalc = sp.result.energy

print(ecalc)
print(qm._result.gradients)
print(sp.result.gradients)
print(qm.frag.bqs.gradients)

# Tcl-ChemShell 3.7.1/ORCA 4.2.1 reference result
# convergence='tightscf'
#eref = -75.59438189
# convergence='normalscf'
#eref = -75.59438188
# ORCA 5.0.2
#eref = -75.59438189
#tolerance = 1.0e-7

#orca reference
#eref = -75.962356680
#pyscf reference 
eref = -75.962387547303

testutils.validate(ecalc, eref)
