from chemsh import Fragment, GULP, Layer, NLayerSubtractive, SP, Opt, DL_POLY, ORCA
from chemsh.utils import testutils
import numpy as np

water_dimer = 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
'''

method_inner = ORCA(functional="blyp", method="dft", basis="6-31G", threshold=1.e-8)
method_outer = GULP(ff=ff)#, molecule=True)

layer_inner = Layer(method_inner, [0, 1, 2])
layer_outer = Layer(method_outer, [3, 4, 5])

oniom = NLayerSubtractive(frag=water_dimer, layers=(layer_inner, layer_outer), embedding="mechanical", coupling="covalent")

opt = Opt(theory=oniom)

opt.run()

# Tcl-Chemshell reference file ./tcl_tests/h2o_dimer_subtractive_mm_orca.chm
eref = -76.377892545745

testutils.validate(opt.result.energy, eref)

