# MgO shell model optimisation

from chemsh import *

# Read in MgO crystal data from disk
mgo = Fragment(coords="mgo.cjson", shellatoms=[0,1,2,3,4,5,6,7])

# Define GULP forcefield for MgO
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
'''

# Set up forcefield calculation
my_theory = GULP(frag=mgo, ff=ff)

# Request a DL-FIND optimisation
opt = Opt(theory=my_theory)
opt.run()
