from chemsh import *
from chemsh.utils import testutils
import numpy as np

# Define the total system
water = Fragment(coords="h2o_dimer.cjson")

# Define the atoms in the low layer (layer_1)
atoms_layer_1  = [3,4,5]
# Define the method to use for the layer_1
method_1 = ORCA(method="hf", basis="3-21g")
# Add to the Layer class
layer_1 = Layer(method_1, atoms_layer_1)

# Define the atoms in the high layer (layer_2)

atoms_layer_2 = [0, 1, 2]
# Define the method to use for the layer_2
method_2 = ORCA(method="dft", functional="blyp", basis="6-31g", maxiter=10000)
# Add to the Layer class
layer_2 = Layer(method_2, atoms_layer_2)

# Use the layers in the two-layer Subtractive class
oniom = NLayerSubtractive(frag=water, layers=(layer_2, layer_1), embedding="mechanical", coupling="covalent")

opt = SP(theory=oniom)

opt.run()

ecalc = opt.result.energy
eref = -151.96389726

testutils.validate(ecalc, eref)
