#
# DL-FIND Transition-state and path search using the nudged elastic band method
#      in DLC coordinates
#
# Open shell example (CH3O)
#
from chemsh       import *
from chemsh.utils import testutils

#
# Example 4 for a transition-state search taken from
# J. Baker and F. Chan, J. Computational Chem. 17, 888 (1996) [bak96a]
#
znums = [ 8, 6, 1, 1, 1 ]

# unit is a.u., since keyword 'unit' has been removed from <class Fragment>
# Re-opt for dftb+
coords1 = [ [  3.03528829e-01, -5.44340443e-13,  3.48475824e-01 ],
            [  7.57045551e-03, -9.79404548e-12,  2.93611590e+00 ],
            [  2.08368816e+00,  0.00000000e+00, -1.70468592e-01 ],
            [ -7.60930605e-01,  1.77602238e+00,  3.73911305e+00 ],
            [ -7.60930605e-01, -1.77602238e+00,  3.73911305e+00 ] ]

reactant = Fragment(znums=znums, coords=coords1)

# unit is a.u., since keyword 'unit' has been removed from <class Fragment>
# Re-opt for dftb+
coords2 = [ [ -9.70792445e-01, -8.81489908e-12, -1.08308640e-01 ], 
            [  1.32071035e-01, -3.52317579e-12,  2.18060400e+00 ], 
            [  2.23765669e+00, -8.56735837e-13,  1.86468382e+00 ], 
            [ -2.63004522e-01,  1.74193917e+00,  3.32768502e+00 ], 
            [ -2.63004522e-01, -1.74193917e+00,  3.32768502e+00 ] ] 

product = Fragment(znums=znums, coords=coords2)


qm = DFTBplus(frag=reactant, mult=2, skf_path="../../dftbplus/slakos/mio-1-1/")

opt = Opt(theory=qm,
          frag2=product,
          neb='frozen',
          coordinates='dlc',
          nimages=10,
          nebk=0.01,
          tolerance=0.0003,
          maxcycle=40000,
          maxene=30000,
          maxstep=0.9,
          trust_radius='constant'
         )

opt.run()

ecalc = opt.result.energy

# Py-ChemShell reference value
eref = -6.008970276800
tolerance = 1.0e-6


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

