Ionisation Potential

Note

Ths scripts for this step in the tutorial can be found in mgo_surface/7_mgo_ip.

This tutorial will run through how to calculate the ionisation potential of a system. As with the previous theme, this tutorial will be delivered by using an MgO cluster as the example.

First, the energy of the non-ionised MgO system is calculated. This is the same as the earlier single point calculation of MgO, except run at the optimised cluster geometry:

# 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

Next, we increase the charge by 1, representing the removal of an electron from MgO. This also changes the environment to an open shell system, with a multiplicity change from a singlet (1) to a doublet (2).

# 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

The ionisation potential of MgO is calculated using the equation below:

\[E_{IP} = E_{[MgO_{surf}]^{+}} - E_{[MgO_{surf}]^0} + E_{corr}\]

Note

The calculation should result in \(E_{IP}= 6.35\) eV.

A more accurate calculation can be obtained if the Jost correction is applied to the script. The Jost correction is applied when dealing with charged species in QM/MM calculations, as the relaxed region has a finite space. The Jost correction uses the following equation:

\[E_{corr} = -\frac{{{Q}^{2}(\varepsilon - 1))}}{2R(\varepsilon + 1)}\]

Where:

Q = total charge of the electrons in the defect, units e.

varepsilon = dielectric constant of MgO, no units.

R = radius of the relaxed region (\(a_(0)\)), units a.u. bohr

Note

The corrected value should be \(E_{IP}= 5.88\) eV.