# 
# Water dimer
#

from numpy        import array
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 = NWChem(method="hf", basis="3-21g")
mm = GULP(ff=ff, molecule=True)

qmmm = QMMM(frag=water, qm_region=range(3), qm=qm, mm=mm, embedding='mechanical')

sp = SP(theory=qmmm, gradients='fd', displace=1.0E-6)

sp.run()

ecalc = sp.result.energy
gcalc = sp.result.gradients

print(" >>> gradients:\n")
print(sp.result.gradients)
print("\n\n")

# Test result
eref = -75.5926115099
gref = [[  3.77349762e-03,   2.58820743e-15,   3.04371695e-02],
        [ -7.87457236e-04,   2.70749616e-02,  -1.46636432e-02],
        [ -7.87457236e-04,  -2.70749616e-02,  -1.46636432e-02],
        [  2.86284787e-02,   0.00000000e+00,  -7.47555424e-03],
        [ -2.50691700e-02,   0.00000000e+00,  -7.87380597e-03],
        [ -5.75789184e-03,   0.00000000e+00,   1.42394770e-02]]
gref = array(gref)

tolerance = 1.0e-7


testutils.validate(gcalc, gref)
