.. _qm_mm_spoint: ******************************* 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: .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_nwchem.basis .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_nwchem.ecp 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 :doc:`ffield` tutorial. The force field used in this tutorial, ``mgo_2body.ff``, is shown below: .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_2body.ff 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**. .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py :lines: 7-8 Next, we read in the QM region and charge as defined in the saved fragment. .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py :lines: 10-15 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. .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py :lines: 18-23 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. .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py :lines: 25-26 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. .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py :lines: 28-33 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. .. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py :lines: 35-37 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.