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

butanol = Fragment(coords="./butanol_mech_opt.cjson")

atoms_layer_1 = list(range(5, butanol.natoms))
method_1 = ORCA(method="hf", basis="3-21G")
layer_1 = Layer(method_1, atoms_layer_1)

atoms_layer_2 = list(range(5))
method_2 = ORCA(method="dft", functional="blyp", basis="6-31g", maxiter=10000)
layer_2 = Layer(method_2, atoms_layer_2)

job = NLayerSubtractive(frag=butanol, layers=(layer_2, layer_1), embedding="mechanical", coupling="covalent")

opt = Opt(theory=job, tolerance=0.0001)
opt.run()
ecalc = opt.result.energy
# From TCL CHEMSHELL (modified) and ORCA 5.0
eref = -232.100175056407

testutils.validate(ecalc, eref)
