# 
# Water dimer
#

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
O1 core O2 shel 582000.0 595.0 99.9
O2 shel O2 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
H2 core O2 core 1059.162 0.96
# Note GULP k is twice that of ChemShell format
# note different ordering of atoms
three bond kcal
O2 core H2 core H2 core 68.087  104.5
# core-shell spring
spring
O2 209.45
'''

qm = FHIaims(functional='pbe', basis='fhiaims.basis.shells',
             ghost_species='Emptium') # settings='fhiaims.settings')

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

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

sp = SP(theory=qmmm, gradients=True)

sp.run()

ecalc = sp.result.energy

# Test result from ChemShell 3
eref = -76.389453044790 

testutils.validate(ecalc, eref)
