Ionisation potentials & Deprotonation energy
7. Ionisation Potential of silicate
The ionisation potential of the mfi cluster can be calculated by finding the single-point energy of the cluster before and after an electron has been removed from an O atom in the qm_region. The difference between these energies must be the energy supplied to ionise the cluster.
This calculation can be done by running the input file mfi_ip.py in the directory zeol_covalent/3_silicate_ip.
In the code, the single-point energy of the neutral silicate is calculated first, using the same method as in section 3, with a few minor differences.
When creating the neutral mfi cluster fragment, there is an additional totalcharge parameter, which is subject to change when an electron is removed:
mfi = Fragment(coords='optimized_mfi.pun', connect_mode='covalent', totalcharge=0)
The value of totalcharge is accessed from the mfi cluster fragment and saved as an integer variable in charge, meaning that we can access this charge elsewhere:
qm_charge = int(mfi.totalcharge)
When defining qm_theory, the NWChem parameters now include charge as well as mult (multiplicity). This is because it is only the qm_region of the cluster that is affected during the ionisation procedure. For the unionised fragment, charge is neutral, and mult is 1:
qm_theory = NWChem(method='scf', basis='3-21g', maxiter=100, charge=qm_charge, mult=1)
After calculating the sp energy of the neutral mfi cluster fragment, it is saved as a variable using the following line.
ecalc_neutral = sp.result.energy
Next, the mfi_ionised cluster fragment is created, and the totalcharge parameter has changed:
mfi_ionised = Fragment(coords='optimized_mfi.pun', connect_mode='covalent',
totalcharge=1)
The qm_theory NWChem parameters charge and mult have also changed to 1 and 2 respectively now that an electron has been removed from the qm_region of the cluster. The charge parameter is accessed and altered using qm_theory.charge, and the mult parameter is accessed and altered using qm_theory.mult.
qm_theory.charge = 1
qm_theory.mult = 2
The qmmm cluster fragment frag is accessed and altered to mfi_ionised using qmmm.frag
qmmm.frag = mfi_ionised
After calculating the sp energy of the ionised mfi cluster fragment, it is saved as a variable using the following line.
ecalc_ionised = sp_ionised.result.energy
Using both of the single-point energies we have aquired before and after ionisation saved as ecalc_neutral and ecalc_ionised, the ionisation potential is calculated using the following equation:
The correction term represents the Jost correction. A function, jost_corr_bulk(), is used to calculate this.
The ionisation potential is first calculated without the Jost correction, and stored in a variable ionisation_potential. A multiplication factor of 27.21138602 is required to display the energy in eV instead of a.u.
ionisation_potential = ecalc_ionised - ecalc_neutral
print("Ionisation potential = ", ionisation_potential*27.21138602, " eV")
After this, the parameters are defined that are required for the jost_corr_bulk() function.
Q = 1 # Total charge of the defect (electrons).
epsilon = 1.71 # This is the relative response to how the atoms in the cluster respond to the change in charge
R = 10 # Radius of relaxed region. This is in a.u. Bohr. Equal to radius of active region.
These parameters are used to find the Jost correction using the jost_corr_bulk() function. The correction is saved in the variable jost_corr.
jost_corr = jost_corr_bulk(Q,epsilon,R)
Finally, the ionisation potential is calculated and displayed with the Jost correction applied using the variables ionisation_potential and jost_corr.
print("Using this correction, the more accurate value of the ionisation potential = ", (ionisation_potential*27.21138602) + jost_corr, " eV")
8. Ionisation Potential of zeolite
The ionisation potential of zeolite can be calculated by running the input file zeol_ip.py in the directory zeol_covalent/8_zeol_ip.
mod:zsm5 and zsm5_ionised are the cluster fragments used instead of mfi and mfi_ionised. The rest of the python script is the same as the one described in section 7.
The Jost correction for zsm5 uses the same constants as mfi also.
9. Deprotonation Energy of zeolite
Due to the H+ cation in the zeolite cluster, the material becomes very acidic, and so zsm5 can be used for acid-catalysed reactions. By calculating the deprotonation energy of the zeolite cluster, this will give us an indication as to how acidic the cluster is.
The deprotonation energy of the zsm5 cluster can be calculated by finding the single-point energy of the cluster before and after deprotonation energy. The difference between these energies must be the energy supplied to deprotonate the cluster.
The deprotonation energy can be calculated by running the input file deprotonation_energy.py in the directory zeol_covalent/9_deproton_energy.
To begin with, the single-point energy of a zsm5 fragment is calculated in the same way as section 3, with a few minor differences.
To create a neutral zsm5 cluster fragment, we include an additional totalcharge parameter, which is subject to change when a proton is removed:
zsm5 = Fragment(coords='zsm5.pun', connect_mode='covalent', totalcharge=0)
The value of totalcharge is accessed from the zsm5 cluster fragment and saved as an integer variable in charge, meaning that we can access this charge elsewhere:
qm_charge = int(zsm5.totalcharge)
When defining qm_theory, the NWChem parameters now include charge as well as mult. This is because it is only the qm_region of the cluster that is affected during the deprotonation procedure. For the regular cluster fragment, charge is neutral, and mult is 1:
qm_theory = NWChem(method='scf', basis='3-21g', maxiter=100, charge=qm_charge, mult=1)
The rest of the sp calculation is carried out as normal, and the value of the sp energy is saved in the variable ecalc_regular.
ecalc_neutral = sp.result.energy
As we are going to continue using this zsm5 cluster fragment to produce the deprotonated cluster, we must revert back the charges changed by the dipole-adjust by changing the bond_modifiers parameter in the qmmm object to reverse_silicate_modifiers.
qmmm.boundary_charge_adjust(reverse_silicate_modifiers)
Now, we must deprotonate the zsm5 cluster. This is achieved using the delete() function on the zsm5 cluster fragment. The function uses positional arguments, and as the H atom was the final atom added during the zeolite construction procedure, it is the 551st atom and therefore has position 550 in the fragment.
zsm5.delete(550)
We must amend to totalcharge parameter in the zsm5 fragment to reflect the fact that the structure is now proton-deficient.
zsm5.totalcharge=-1
After these two changes, a deprotonated zeolite now exists within the zsm5 cluster fragment.
The value of totalcharge is accessed from the zsm5 cluster fragment and saved as an integer variable in charge, meaning that we can access this charge elsewhere:
qm_theory.charge = int(zsm5.totalcharge)
The qm_region is updated to reflect the fact that there is no longer a H atom bound to an O atom in the region anymore. This is again carried out by using the getRegion() function, which adds any atoms with a ‘1’ suffix to the qm_region.
qm_region = zsm5.getRegion(1)
The qmmm object and the atoms frozen during optimization are redefined, which are inserted as parameters to the opt object.
The optimized energy of the deprotonated zsm5 cluster fragment is saved as a variable using the following line.
ecalc_deprotonated = opt.result.energy
Now that we have the energy of the regular and deprotonated structure, ecalc_regular and ecalc_deprotonated, we can find the deprotonation energy. The equation is:
The correction term is the same as the ionisation potential, the Jost correction. The function, jost_corr_bulk(), is used to calculate this.
The deprotonation energy is first calculated without the Jost correction, and stored in a variable deprotonation_energy. A multiplication factor of 27.21138602 is required to display the energy in eV instead of a.u.
deprotonation_energy = ecalc_deprotonated - ecalc_regular
print("The deprotonation energy is: ", deprotonation_energy*27.21138602, " eV")
After this, the parameters are defined that are required for the jost_corr_bulk() function. As you can see, the dielectric constant and the radius of the relaxed region parameters are unchanged when compared to the ionisation potential Jost correction parameters, though the charge has changed from 1 to -1.
Q = -1 #Charge. This is in units e.
epsilon = 1.71 #Dielectric Constant. This is unitless.
R = 10 #Radius of relaxed region. This is in a.u. Bohr. Equal to radius of active region
These parameters are used to find the Jost correction using the jost_corr_bulk() function. Again, the correction is saved in the variable jost_corr.
jost_corr = jost_corr_bulk(Q,epsilon,R)
Finally, the deprotonation energy is calculated and displayed with the Jost correction applied using the variables deprotonation_energy and jost_corr.
print("Using this correction, the more accurate value of the deprotonation energy = ", (deprotonation_energy*27.21138602) + jost_corr, " eV")
A perhaps more accurate value can be recited from the paper Modelling metal centres, acid sites and reaction mechanisms in microporous catalysts. Any difference between the value from the paper (~1100kJ/mol) and the value calculated in the script (~1150kJ/mol) may be due to the size of the default qm region.