#
# 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
'''


mm = GULP(frag=water, ff=ff, molecule=True)

sp = SP(theory=mm, gradients='fd')

sp.run()

ecalc = sp.result.energy
gcalc = sp.result.gradients
print(" >>> Gradients:\n")
print(sp.result.gradients)
print("\n\n")

# analytical gradients
gref = [[ 0.0037735 ,  0.        ,  0.02740732],
        [-0.00078746,  0.01281475, -0.01314872],
        [-0.00078746, -0.01281475, -0.01314872],
        [ 0.02862848,  0.        , -0.00747555],
        [-0.02506917,  0.        , -0.00787381],
        [-0.00575789,  0.        ,  0.01423948]]

gref = array(gref)

testutils.validate(gcalc, gref, tol=1.0e-6)
