Single-Point Energy Calculation
Note
The scripts for this step of the tutorial can be found in mgo_surface/2_mgo_sp.
Now that the embedded cluster model has been set up, it is ready for QM/MM simulations. In this tutorial, we demonstrate how a single-point QM/MM energy evaluation can be performed on the model structure.
1. Quantum Mechanics (QM)
Note
The QM calculations in this example use NWChem.
As with all QM packages, an NWChem calculation requires a basis set to be specified. Basis sets for the MgO system can be obtained from the Basis Set Exchange.
For Mg, we will use the Stuttgart relativistic, large core basis which is designed to be used with an effective core potential (ECP). The ECP is a representative potential which replaces some core electrons, thereby meaning a smaller number of basis functions can be used to represent the system.
Note
Technical note: the Stuttgart basis was modified to remove the most diffuse functions in order to help prevent artificial spreading of charge density outside the QM region.
For O, we will use a standard Ahlrichs pVDZ double zeta basis with polarisation functions. We will also use this basis for the adorbate in later stages of the tutorial.
The basis set an ECP specifications are provided to ChemShell as separate .basis and .ecp files in NWChem format.
Below are the mgo_nwchem.basis and mgo_nwchem.ecp files used:
basis cartesian
#BASIS: Ahlrichs pVDZ
O1 S
2266.1767785 -0.53431809926E-02
340.87010191 -0.39890039230E-01
77.363135167 -0.17853911985
21.479644940 -0.46427684959
6.6589433124 -0.44309745172
O1 S
0.80975975668 1.0000000000
O1 S
0.25530772234 1.0000000000
O1 P
17.721504317 0.43394573193E-01
3.8635505440 0.23094120765
1.0480920883 0.51375311064
O1 P
0.27641544411 1.0000000000
O1 D
1.200000 1.000000
#BASIS: Modified Stuttgart RLC
Mg1 S
2.42571930 0.02676400
0.82262500 -0.22388000
Mg1 S
0.10774900 0.62046400
Mg1 P
0.76904700 -0.03664800
0.18867500 0.24314500
Mg1 P
0.07510100 0.55478400
end
ecp
Mg1 nelec 10
Mg1 ul
2 1.00000000 0.00000000
Mg1 s
2 1.73200000 14.67600000
Mg1 p
2 1.11500000 5.17570000
Mg1 d
2 1.20300000 -1.81600000
bq_Mg2_e nelec 0
bq_Mg2_e ul
2 1.00000000 0.00000000
bq_Mg2_e s
2 1.73200000 14.67600000
bq_Mg2_e p
2 1.11500000 5.17570000
bq_Mg2_e d
2 1.20300000 -1.81600000
end
The ECP fulfills two purposes in the ionic embedded cluster model:
In the inner QM region (Mg1), it is used in a conventional way to replace core electrons (specified by
nelec 10, calculated as 12 Mg electrons minus 2 from the valence orbital)In the QM/MM boundary region (bq_Mg2_e), the ECP acts as a boundary potential to prevent charge density diffusing outside the QM region. In the boundary region, no full QM atoms are present so no core electrons are replaced (as specified by
nelec 0).
With the exception of the number of electrons replaced, the ECP definitions for region 1 and region 2 are identical.
2. Molecular Mechanics (MM)
Note
The MM calculations in this example use GULP. We use the term MM here loosely, as is usual in the QM/MM literature. The interatomic potentials used in the ionic embedding model are quite different in form to molecular mechanics forcefields used for covalently-bonded systems, but they fulfill the same purpose.
For the classical environment, we need to specify a forcefield in GULP format.
This can be sourced in one of two ways: either using a pre-existing forcefield
from the literature or generated manually using the approach described in the
Appendix: Generating a Force Field tutorial. The force field used in this tutorial,
mgo_2body.ff, is shown below:
# region 1-3
buck
Mg1 core O3 shel 798.394448 0.323866 0.000 0.000 12.000
lennard 12 6
Mg1 core O3 shel 182.246312 0.000000 0.000 12.000
buck
O1 core O3 shel 22764.300 0.149000 27.879 0.000 12.000
# region 1-4
buck
Mg1 core O4 shel 798.394448 0.323866 0.000 0.000 12.000
lennard 12 6
Mg1 core O4 shel 182.246312 0.000000 0.000 12.000
buck
O1 core O4 shel 22764.300 0.149000 27.879 0.000 12.000
# region 2-2
# region 2-3
buck
Mg2 core O3 shel 798.394448 0.323866 0.000 0.000 12.000
lennard 12 6
Mg2 core O3 shel 182.246312 0.000000 0.000 12.000
# region 2-4
buck
Mg2 core O4 shel 798.394448 0.323866 0.000 0.000 12.000
lennard 12 6
Mg2 core O4 shel 182.246312 0.000000 0.000 12.000
# region 3-3
buck
Mg3 core O3 shel 798.394448 0.323866 0.000 0.000 12.000
lennard 12 6
Mg3 core O3 shel 182.246312 0.000000 0.000 12.000
buck
O3 shel O3 shel 22764.300 0.149000 27.879 0.000 12.000
spring
O3 46.097205 3589.6947
# region 3-4
buck
Mg3 core O4 shel 798.394448 0.323866 0.000 0.000 12.000
lennard 12 6
Mg3 core O4 shel 182.246312 0.000000 0.000 12.000
buck
Mg4 core O3 shel 798.394448 0.323866 0.000 0.000 12.000
lennard 12 6
Mg4 core O3 shel 182.246312 0.000000 0.000 12.000
buck
O3 shel O4 shel 22764.300 0.149000 27.879 0.000 12.000
# exclusions of interactions in and between regions 4 and 5
coulomb
Mg4 core Mg4 core 0.0 100.0
Mg4 core O4 core 0.0 100.0
Mg4 core O4 shel 0.0 100.0
O4 core O4 core 0.0 100.0
O4 core O4 shel 0.6 100.0
O4 shel O4 shel 0.0 100.0
F5 core F5 core 0.0 100.0
F5 core Mg4 core 0.0 100.0
F5 core O4 core 0.0 100.0
F5 core O4 shel 0.0 100.0
# GULP CONTROL PARAMETERS
xtol opt 6.000000
gtol opt 6.000000
ftol opt 8.000000
maxcyc 5
stepmx 1.00
output xyz current.xyz
accuracy 12.000
If you are familiar with GULP calculations, you will recognise standard
forcefield terms such as Buckingham buck, Lennard-Jones lennard,
and Coulomb coulomb interatomic potentials. There are also Spring spring
terms to control the interaction between classical cores and shells.
More information on setting these terms can be found in the
GULP manual.
The forcefield description also depends on which region the atoms are in. The parameters for each term do not vary by region, but certain terms are not necessary. For example, there are no O ions in region 2, so there are no region 2-2 terms (as there are no Mg-Mg terms in the forcefield). All terms involving O ions are between shells, except for region 1 O ions (where shells have been deleted). Similarly, spring terms are only necessary for region 3. Finally, interactions within regions 4 and 5 (including shell interactions) are excluded as these regions are fixed.
3. Calculating the Single-Point Energy
The input structure used in this calculation is the mgo_shells_region.pun
created previously. The coordination table used in the
cluster creation step is no longer required, and so it can (optionally)
be deleted by setting connect_mode to None.
# Import MgO cluster geometry from previous stage
mgo = Fragment(coords='mgo_shells_regions.pun', connect_mode=None)
Next, we read in the QM region and charge as defined in the saved fragment.
# 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("QM charge is: ", qm_charge)
The getRegion() method is used to obtain a list of atoms with the
specified region label (here region 1, the QM region).
For the QM charge, we round the floating point value in the saved fragment to the nearest integer for use in the QM calculation.
The settings for the QM/MM calculation can then be defined. First, we set the options for the QM (NWChem) calculation.
qm_theory = NWChem(method='dft',
functional='b3lyp',
basis='mgo_nwchem.basis',
ecp='mgo_nwchem.ecp',
charge=qm_charge,
path='mpirun -np 4 nwchem')
Here, we have specified a density functional theory calculation using the well-known B3LYP functional, with the basis set and ECP files we discussed above. The saved QM charge is passed in. As this is a singlet calculation, the spin multiplicitly can be left as the default (1).
Special consideration should be given to the path option, which
we have set here to call mpirun in order to run NWChem in parallel with
4 cores. If you do not have an MPI environment set up on your machine,
you should remove the
path option so that NWChem can be run directly. The calculation
will of course take longer as a result. On the other hand, if you have
more than 4 cores available, you could experiment with increasing the
number of cores passed to mpirun.
Note
These instructions assume you are running ChemShell and NWChem as separate executables. If you are running ChemShell with NWChem directly-linked, however (as will be the case on most HPC platforms), the path option should be removed as the NWChem calculation will use the same MPI environment as ChemShell.
Next, the options for the GULP part of the calculation are given.
mm_theory = GULP(ff='mgo_2body.ff',
conjugate=True)
These are quite straightforward, specifying the file containing the forcefield definition, and the optimisation algorithm for shell relaxation.
Finally, the options for the QM/MM driver are given.
qmmm = QMMM(qm=qm_theory,
mm=mm_theory,
frag=mgo,
qm_region=qm_region,
coupling='ionic',
shl_maxcycles=20)
Here, the QM and MM theory options specified above are passed in, together with the QM region definition obtained above, and the QM/MM coupling model used (in the case of MgO, the ionic embedding model). Finally, we set the maximum number of shell relaxation cycles to the relatively high value of 20, to ensure that the shells are fully converged for the final energy evaluation.
The QM/MM calculation is now ready to run. We set up a single-point
calculation by passing the QM/MM options to the SP() class, and then call
the run method to run the calculation.
# Calculate the QM/MM energy
sp = SP(theory=qmmm, gradients=False)
sp.run()
The remainder of the input file is a validation step to check that the calculation you obtain is the same as our pre-calculated value. This is purely for testing purposes and would not normally form part of a production calculation.
Note
If your calculated single-point energy differs significantly from the pre-calculated value and so fails the test, it may indicate that there is a problem with your ChemShell, NWChem or GULP installations.