'''MgO shell model single point'''

from chemsh       import *
from chemsh.utils import testutils

mgo = Fragment(coords="mgo.cjson", shellatoms=range(8))
#mgo = Fragment(coords="mgo.cjson", shellatoms=['Mg','O'])

ff = '''species
Mg core  1.580
Mg shel  0.420
O  core  0.513
O  shel -2.513
buckingham
Mg shel O shel 2457.243 0.2610  0.00 0.0 10.0
O  shel O shel   25.410 0.6937 32.32 0.0 12.0
spring
Mg 349.95
O   20.53
'''

print("### mgo.nshells =", mgo.nshells)
print("### mgo.shells.coords =", mgo.shells.coords)

mm = GULP(frag=mgo, ff=ff)

sp = SP(theory=mm, gradients=True)
sp.run()

ecalc = sp.result.energy

# Test result
eref = -5.98354835


testutils.validate(ecalc, eref)
