#
# Water dimer with oxygen shells
#

from chemsh       import *
from chemsh.utils import testutils

water = Fragment(coords='h2o_dimer_shells.cjson')

ff = '''#
# molmec removes Coulomb terms between bonded atoms
# and between atoms bonded to a common third atom
keyword molmec
# Flexible TIP3P model with ad-hoc oxygen shell
# (purely for testing purposes)
# last number in LJ term is a cutoff
lennard 12 6 kcal x13
O shel O shel 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 core O core 1059.162 0.96
# Note GULP k is twice that of ChemShell format
# note different ordering of atoms
three bond kcal
O core H core H core 68.087  104.5
# core-shell spring
spring
O 209.45
'''


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

opt = Opt(theory=mm,
          algorithm="lbfgs", tolerance=0.0003, trust_radius="energy",
          maxcycle=400, maxene=300).run()

ecalc = opt.result.energy

# Test result from ChemShell 3
eref = -0.011496334090
tolerance = 1.0e-5


testutils.validate(ecalc, eref, tol=tolerance)
