Setting up the QM/MM Cluster

Note

The scripts for this step of the tutorial can be found in mgo_surface/1_cut_cluster.

Preparation of a surface cluster

This section of the tutorial will outline how to set up the QM/MM cluster using the input script construct_cluster.py.

The periodic crystalline structure of interest is first read in to ChemShell, which is stored as a 2D slab in the file mgo_shells_periodic.pun in ChemShell’s punch file format. This is imported by construct_cluster.py as follows:

# Import periodic MgO slab geometry
mgo = Fragment(coords='mgo_shells_periodic.pun', connect_mode='ionic')

where mgo is a Fragment containing the periodic input structure in mgo_shells_periodic.pun. The connectivity (denoted as connect_mode) should be set to ionic, as MgO is an ionic system.

Note

Periodic inputs can also be read in from CIF files using the ChemShell/ASE interface, with an example given in mgo_cluster_cif.py.

The cluster is then built using the construct_cluster() method:

# Cut out cluster from periodic fragment
# Note distances are in bohr and python counts atoms from zero
cluster = mgo.construct_cluster(origin=40, radius_cluster=40.0,
                                adjust_charge='coordination_scaled',
                                radius_active=24.0, bq_margin=2.0, 
                                bq_density=5, bq_layer=12.0)

When this command is run ChemShell will first cut the cluster from the crystal, which is controlled by the following options:

  • origin dictates where the cluster is centred on, normally one of the atoms at the surface (note that in Py-ChemShell atoms are counted from zero). The atoms in the resulting cluster will be ordered by increasing distance from the origin.

  • radius_cluster determines the overall size of the cluster in atomic units.

Next, the ionic charges on the outermost layer of the generated cluster are adjusted to account for the loss of their near neighbours when they were cut from the slab. As MgO is a straightforward two species system, the charges can simply be scaled by their coordination number, which is selected by the adjust_charge option.

Finally, additional point charges are fitted around the outside of the cluster, in order to mimic as closely as possible the periodic electrostatic interactions that are lost when the cluster is cut out of the slab. A number of options control this fit:

  • radius_active sets the size of the electrostatic sampling region from which potential and electric field information is taken for the fit. This should always be set to be at least as large as the largest geometry optimisation region you intend you use.

  • bq_density relates to how many point charge are placed in the shell of charges immediately surrounding the cluster. This should be chosen to ensure a good fit, but not so high that the point charges are overfitted (see below).

  • bq_margin sets the distance of the shell of fitted point charges fitted point charges around the edge of the cluster. It is OK for the point charges to be quite close as the outer edge of ions in the cluster is not part of the active (potential sampling) region.

  • For the surface case, bq_layer adds additional point charge layers above the surface level, to ensure that the electrostatic environment is well represented at the level at which reactions will take place.

There are no hard and fast rules for determining these settings, but the fit procedure will output a number of measures that indicate how well the resulting potential and fields match the original periodic structure. You should examine the fitted charges in the output (under the header “BQ charges”) to check the fit has produced sensible results. The first six are located at a significant distance from the cluster in order to account for long range dipoles, and so can be large. The remainder should, in general, be reasonably small if the fit is of good quality.

A good strategy is to try the default values first, and then to check how sensitive the results are to changes in the input parameters.

Once the cluster has been constructed, the construct_cluster.py script will save it to a new punch file, mgo_shells_cluster.pun, for use in the next stage of the process:

# Save the resulting cluster in ChemShell punch format
cluster.save('mgo_shells_cluster.pun', 'pun')

Note that by default the point charges are denoted (arbitrarily) as F, so they can be understood by external MM software in subsequent calculations. For validation, compare against the pre-calculated example mgo_shells_cluster.pun.REF.

Partitioning the cluster

Next, the cluster needs to be partitioned into regions. The ionic embedding model consists of five numbered regions corresponding to the QM region, QM/MM boundary, active MM, inactive MM, and fitted point charges.

To proceed, it is required that the user identifies the QM region they wish to use. This can either be done manually by relabelling the relevant atoms in the punch file with a suffix of “1”, e.g. Mg becomes Mg1, or by providing a list of atoms to the partitioning command using the option qm_region. In our example we take the latter approach, taking the 14 atoms closest to the origin as our QM region:

# Define the atoms to have in the QM region. 
# Here we choose atoms 0-13 inclusive.
qm_region=range(14)

The embedded cluster is then split into regions using the partition() method:

# Partition the cluster
regions = mgo.partition(qm_region=qm_region, cutoff_boundary=6.0, 
                        radius_active=24.0, origin=origin,
                        qmmm_interface='explicit', interface_exclude=["O"])

Beyond specifying the QM region, the following options should be set to control how the other regions are selected:

  • cutoff_boundary determines which ions are contained within the QM/MM boundary region (region 2). Here we are using the so-called “explicit” interface (specified by qmmm_interface='explicit'), in which the interface ions will be explicitly represented in the MM calculation by charged particles. The interface ions will enter the QM calculation as point charges with effective core potentials attached. The alternative approach, known as an “implicit” interface, is not supported by NWChem and so we do not consider it in this tutorial.

  • interface_exclude is used to exclude certain ion types from the boundary region. In our example we exclude oxygen ions, which then automatically go into the active MM region (region 3). We do this because effective core potentials for the oxygen anions are not readily available, and it is sufficient to have only cations in the boundary region.

  • radius_active defines the region that you want to optimise, i.e. the active MM region (region 3). This can be the same size or smaller than the sampling region used to fit charges around the cluster, but no bigger.

Finally the script will write a new punch file where each atom is labelled with a region number:

# Save the partitioned cluster
regions.save('mgo_shells_regions.pun', 'pun')

As well as relabelling, the partitioning routine also deletes the superfluous classical charges on the QM region atoms, as well as any corresponding shells, and instead provides an overall charge to be used in QM calculations (in this case -8).

Note

As MM partial charges are floating point numbers the calculated QM charge will also be a floating point number, and should be rounded to an integer before use in a QM calculation.