# MgO ionisation potential

# Initialise ChemShell
from chemsh import *
from chemsh.utils import testutils
from chemsh.data import unitconvert
from chemsh.cluster.jost_corr import jost_corr_surf

print("\n-----------------")
print("Un-ionised system")
print("-----------------\n")

# This is our un-ionised fragment
mgo = Fragment(coords='mgo_shells_opt.pun', connect_mode=None)

# Get list of QM atoms and charge from input fragment
qm_region = mgo.getRegion(1)
print("QM region is: ", qm_region)

qm_charge = round(mgo.totalcharge)
print("Un-ionised QM charge is: ", qm_charge)

# QM/MM settings for un-ionised fragment
# This is a singlet calculation as in previous stages
qm_theory = NWChem(method='dft',
                functional='b3lyp',
                basis='mgo_nwchem.basis',
                ecp='mgo_nwchem.ecp',
                charge=qm_charge,
                mult=1,
		path="mpirun -np 4 nwchem")
	
mm_theory = GULP(ff='mgo_2body.ff',
                 conjugate=True)

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

# Calculate the un-ionised QM/MM energy
sp_unionised = SP(theory=qmmm, gradients=False)
sp_unionised.run()
ecalc_unionised = sp_unionised.result.energy

print("\n--------------")
print("Ionised system")
print("--------------\n")

# This is our ionised fragment
ionised_charge = qm_charge + 1
mgo_ionised = Fragment(coords='mgo_shells_opt.pun', connect_mode=None,
                       totalcharge=ionised_charge)

# The ionised fragment is now open shell, so we alter the multiplicity
qm_theory.charge = ionised_charge
qm_theory.mult = 2
print("Ionised QM charge is: ", ionised_charge)

# This our ionised cluster
qmmm.frag = mgo_ionised

# Calculate the ionised QM/MM energy
sp_ionised = SP(theory=qmmm, gradients=False)
sp_ionised.run()
ecalc_ionised = sp_ionised.result.energy

print("\n--------------------")
print("Ionisation Potential")
print("--------------------\n")

# Calculate the ionisation potential.
# Use ChemShell unit conversion data for Hartree -> eV conversion.
ip = (ecalc_ionised - ecalc_unionised) * unitconvert.Ha2eV
print("Ionisation potential = ", ip, " eV")

# Calculate Jost correction
# Charge in units of e
Q = 1 
# Dielectric constant (unitless)
epsilon = 10 
# Radius of relaxed region (a.u.)
R = 24.0
jost_corr = jost_corr_surf(Q,epsilon,R)

# The IP with the Jost correction applied
ip_corr = ip + jost_corr
print("The Jost correction is", jost_corr, " ev")
print("Ionisation potential with Jost correction applied = ", ip_corr, " eV")

# Validate results against pre-calculated reference energies
tolerance = 1.0e-6
# For mgo
eref_unionised = -1083.536372679964
testutils.validate(ecalc_unionised, eref_unionised, tol=tolerance)

# For mgo_ionised
eref_ionised = -1083.303108307366
testutils.validate(ecalc_ionised, eref_ionised, tol=tolerance)
