# QM/MM Geometry optimization calculation using the constructed solid-state Li_MgO cluster

# Initialise ChemShell
from chemsh import *   

# Import the partioned Li_MgO cluster
li_mgo = Fragment(coords='Li_MgO-cluster.pun', connect_mode=None)

# Extracts the qm 'region' and 'charge' parameters and stores them as variables
qm_region = li_mgo.getRegion(1)
print("The QM region is ", qm_region)
qm_charge = round(li_mgo.totalcharge)
print("The QM charge is ", qm_charge)
active_region = li_mgo.getRegion(1,2,3)
print("The active region is ", active_region)

# Defines QM and MM theory settings
qm_theory = NWChem(method='dft', functional='pbe0', 
                   basis='li-mgo_nwchem.basis', ecp='li-mgo_nwchem.ecp', 
                   charge=qm_charge, path='mpirun -np 4 nwchem', pop_analysis='mulliken', maxiter=500)
 
mm_theory = GULP(ff='mgo_2body.ff', conjugate=True)

qmmm= QMMM(qm=qm_theory, mm=mm_theory, 
           frag=li_mgo, qm_region=qm_region, 
           coupling='ionic', shl_maxcycles=5)

# Defining the QM/MM optimization parameters 
opt= Opt(theory=qmmm, algorithm="lbfgs",
         trust_radius="energy", maxcycle=200, 
         maxene=300, active=active_region)

# Running QM/MM geometry optimization
opt.run()

# Saves the optimized Li_MgO cluster in punch and xyz formats to disc
li_mgo.save('Li_MgO-cluster_optimized.pun', 'pun')
li_mgo.save('Li_MgO-cluster_optimized.xyz', 'xyz')




