Solid State Embedding with a Shell Model Potential
This is an illustration of the setup required to perform embedded cluster calculations on a solid surface.
Starting Point: a 2D Periodic Fragment
mgo_surf.chm should be run first to generate the 2D periodic slab from which the hemispherical cluster will be cut out. The slab must be at least as deep as the intended cluster radius, otherwise part of the cluster will be missing. Also, it should be thick enough to give a sensible electrostatic potential and electric field in the active region.
Note that connect ionic is specified to generate the connectivity table using ionic radii (using the standard covalent radii would result in spurious connections).
The periodic fragment created is mgo_periodic.pun.
Preparation of a Hemispherical Cluster
The cluster is prepared using construct_cluster.chm. It is built with the following command:
construct_cluster coords=mgo_periodic.pun cluster=mgo_model.pun crystal_type=ionic \ origin_atom=41 radius_cluster=40.0 radius_active=24.0 \ bq_margin=2.0 bq_density=5 bq_layer=12.0 \ adjust_charge=coordination_scaled
The input and output fragments are specified by coords and cluster respectively.
The MgO system is of the ionic crystal type. The first stage of the command is to cut out the cluster from the crystal. The resulting cluster is centred on the atom specified by origin_atom, which would normally be one of the atoms at the surface. The size of the cluster is specified by radius_cluster (in atomic units).
In the second stage ionic charges on the boundary are adjusted to reflect the fact that the ions there have fewer neighbours. In this case we simply scale by the coordination number using adjust_charge=coordination_scaled.
Finally point charges are added around the edge of the cluster in order to approximate the missing periodic electrostatic interactions. The number and position of charges is set by bq_density, bq_margin and bq_layer, and the size of the sampling region for the fit by radius_active. radius_active should be at least as large as the largest optimisation region you intend to use.
Note that by default the symbol 'F' is (arbitrarily) used to denote the point charges surrounding the cluster.
Full documentation of these options can be found in the Cutting Clusters section of the manual.
Partitioning the Cluster
The next stage is to specify which parts of the hemispherical model cluster are to be treated as QM, interface, active MM or frozen MM regions.
First copy the mgo_model.pun fragment to mgo_model1.pun. The quantum mechanics region (region 1) is then defined by relabelling at least one atom with a 1 suffix (e.g. Mg -> Mg1). You can see an example of a minimal two-atom QM region in mgo_model1.pun.sample.
The other regions are then assigned by running cluster_partition.chm:
cluster_partition coords=mgo_model1.pun output=mgo_regions.pun cluster_origin=cluster_origin \
cutoff_boundary=5.0 radius_active=12.0 qmmm_interface=implicit interface_exclude= { O }
If any Mg ions are within cutoff_boundary a.u. of region 1, they will be assigned to the interface region (region 2). Here we are using the standard 'implicit' interface, in which the interface ions will be represented in the QM calculation through effective core potentials placed on QM atoms. We choose to exclude the oxygen ions from region 2 because the anionic ECPs required are not readily available. Instead the oxygen ions will be placed in region 3.
The active region (region 3) is defined using the distance radius_active from the cluster origin. Note here we use a smaller active region than previously specified just to speed up the example optimisation. The remaining atoms are assigned to the frozen region (4) and correcting charges to region 5.
In the resulting output mgo_regions.pun the atoms are re-ordered by region and the region number is appended to the element symbol for every atom. This can be used to control forcefield and basis specifications for the QM/MM calculation.
In addition to the relabelling, cluster_partition also removes the shells from atoms in the quantum and interface regions, as these are not used in subsequent QM/MM calculations.
Definition of an Embedded Cluster Calculation
To use the embedded cluster we have set up you need to consider both QM and MM code input as a function of the different regions.
GULP specification
The GULP parameter file is shown below. Note in particular that GULP treats the point charges as if they were normal atoms, using the label F5.
NB: The shell optimisation tolerances are set much higher than normal simply to speed up the example optimisation.
# 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
GAMESS-UK specification
The files mgo.basis and mgo.ecp contain GAMESS-UK directives to specify the basis and ECP details for the atom types in the model. Note that region 1 atoms have both a basis set and pseudopotential, while region 2 atoms only have pseudopotentials. GAMESS-UK may print a warning about this which can be safely ignored.
Running the calculation
The energy calculation is executed in example_ener.chm. The first step is to delete all connectivity, achieved by setting the connectivity control parameters to 0. This is necessary so that the QM/MM code does not try to terminate any bonds between the QM and MM regions using link atoms.
global chemsh_default_connectivity_toler global chemsh_default_connectivity_scale set chemsh_default_connectivity_toler 0.0 set chemsh_default_connectivity_scale 0.0 connect coords=$job.pun conn=$job.pun
Next the QM region atom list and the charge to be set in GAMESS-UK are extracted from the input fragment (mgo_regions.pun). The get_cluster_region command can be used to query any of the regions in the fragment, and we use it again in the optimisation example to get the active atoms list.
set qm_region [ get_cluster_region coords=$job.pun region= {1 2} ]
set qm_charge [ get_molecule_charge coords=$job.pun ]
set qm_charge [ expr round($qm_charge) ]
The QM/MM arguments are then set including the options to be supplied to GAMESS-UK and GULP. In the arguments to GULP, include_qm_force is needed to switch on electrostatic coupling of the QM potential to the shells.
set gulp_arguments " mm_defs=mgo_2body.ff \
include_qm_force=yes \
use_second_derivatives=no "
The arguments to GAMESS-UK include specification of the charge (stored in qm_charge), DFT functional, and the location of the pre-prepared basis set and ECP files:
set gamess_arguments " basisfile=mgo.basis \
ecpfile=mgo.ecp \
maxcyc=100 \
functional=b3lyp \
direct=yes \
echo_input=no print_bqs=no \
unique_listing=yes \
charge=$qm_charge \
scftype=rhf \
mult=1 \
memory=16000000 "
Symmetry should be turned off to ensure that the quantum region is not reorientated by GAMESS-UK as the calculation proceeds. Direct SCF is generally used for all medium to large clusters.
The arguments to hybrid are as follows:
set args "qm_region= [ list $qm_region ] \
mm_theory = gulp : [ list $gulp_arguments ] \
qm_theory = gamess : [ list $gamess_arguments ] \
electrostatics_option=polarised \
name_bqs=yes \
shell_restart=no \
calculate_connectivity=no \
max_shell_cycles=20 "
The QM region is specified using the list stored in $qm_region. The name_bqs setting is optional; it labels each bq in the generated GAMESS-UK input file which helps if you need to check the input. The value of 20 for max_shell_cycles will ensure full convergence in the iterative optimisation of shells and QM charge density.
Finally, the calculation is run to obtain the QM/MM energy:
energy energy=energy coords=$job.pun theory=hybrid : [list $args]
In the tutorial directory there is also an example QM/MM optimisation (example_opti.chm) using DL-FIND. For geometry optimisations a smaller value of max_shell_cycles can be used to accelerate the calculation, as each point will start from a nearly converged set of shell positions.
Further embedding options
The steps above describe the default embedding method which is suitable for GAMESS-UK and GULP. In this section we discuss refinements to embedding in ChemShell v3.6 to support more flexible models.
QM/MM interfaces types: implicit and explicit
In the 'implicit' QM/MM interface used above, the interface (region 2) is represented by atoms in the QM calculation to which effective core potentials are attached but no basis set. In the MM code, the region 2 centres, like the region 1 QM centres, have no charge associated with them.
An alternative - but equivalent - approach is the 'explicit' interface, where the interface ions are explicitly represented in the MM code by charged particles. In this case the interface centres enter the QM calculation as point charges (bqs), albeit with effective core potentials attached. In this case the nuclear charge in the effective core potentials of the region 2 ions should be set to zero.
The explicit interface is selected at the cluster partitioning stage, as in the example cluster_partition_explicit.chm. If you compare the output mgo_regions_explicit.pun with the original mgo_regions.pun, you can see that the MM charges on the Mg2 are retained, with a corresponding change in the overall QM charge from +8 to zero.
Next, before performing a calculation the MgO ECP file needs to be modified as in mgo_explicit.ecp so that the nuclear charge on the Mg2 ECP is zero (note also the change in label to 'bq_Mg2', reflecting the change in the way the centre is handled in the QM calculation).
An example of an explicit interface QM/MM calculation is given in example_ener_explicit.chm. It is instructive to compare this with the implicit version (example_ener.chm). The QM region is now specified as region 1 only, reflecting the change in the nature of the boundary region:
set qm_region [ get_cluster_region coords=$job.pun region= {1} ]
The change of the QM charge is automatically handled, being extracted from the input fragment mgo_regions_explicit.pun. The only other change is to set ecpfile= mgo_explicit.ecp. The output energy should be in good agreement with the implicit case (-156.212 a.u.).
Separable boundary ECPs
When using the explicit interface, it is possible to separate the boundary region point charges and embedding ECPs onto different logical centres using the option separate_boundary_ecps in the hybrid interface.
In the GAMESS-UK case, the effective core potential centres are represented as further bqs (point charges), with the suffix '_e', i.e. bq_Mg2_e for our example. A nominal charge of 10-8 is assigned to the centre and the same amount subtracted from the interface ion point charge. Both centres share the same physical location, so the physical model has not changed.
To assign the embedding ECPs to the correct centres, the ECP definition file should be changed again so that the label reads 'bq_Mg2_e' instead of 'bq_Mg2', as in mgo_sepecp.ecp. An example calculation is given in example_ener_sepecp.chm. The only changes are the addition of:
separate_boundary_ecps=yes \
in the list of hybrid arguments and the change in ECP file name. The output energy should be in good agreement with the standard explicit case.
For the example calculation considered here there is clearly no special advantage to using the explicit interface and separate boundary ECPs with GAMESS-UK, as we get the same energy no matter which approach we take. However, these options do allow for more flexible boundaries, in particular the attachment of embedding ECPs to anion centres (which in the example above we explicitly exclude from the boundary at the partitioning stage using the option interface_exclude). The separate_boundary_ecps option gets around the limitations of the QM code with respect to anions because the ECP is actually assigned to a very small but positive point charge that shares the same location as the anion.
Embedded cluster calculations with NWChem
The new options are particularly important for NWChem because for technical reasons NWChem only supports embedded cluster calculations using the explicit interface and separate boundary ECPs. You should ensure that the cluster is set up accordingly, as the NWChem interface cannot detect an incorrect setup.
The NWChem example, example_ener_nwchem.chm therefore takes a similar form to example_ener_sepecp.chm. The basis set and ECP files, mgo_nwchem.basis and mgo_nwchem.ecp are in a different format to GAMESS-UK but are otherwise identical. The resulting energy should again be in good agreement (-156.212 a.u.), with minor differences due to different QM code defaults.