12.3.2 Using the Stand-Alone QM/MM Features

To perform QM/MM calculations, the user must assign MM atom types for each atom in the $molecule section. The format for this specification is modeled upon that used by the Tinker molecular modeling package,⁷⁸³ although the Tinker program is not required to perform QM/MM calculations using Q-Chem. Force field parameters and MM atom type numbers used within Q-Chem are identical to those used Tinker for the AMBER99, CHARMM27, and OPLSAA force fields, and the format of the force field parameters files is also the same.

The $molecule section must use Cartesian coordinates to define the molecular geometry for internal QM/MM calculations; the Z-matrix format is not valid. MM atom types are specified in the $molecule section immediately after the Cartesian coordinates on a line so that the general format for the $molecule section is

$molecule
  <Charge> <Multiplicity>
  <Atom> <X> <Y> <Z> <MM atom type>
  . . .
$end

For example, one can define a TIP3P water molecule using AMBER99 atom types, as follows:

$molecule
   0 1
   O        -0.790909    1.149780    0.907453  2001
   H        -1.628044    1.245320    1.376372  2002
   H        -0.669346    1.913705    0.331002  2002
$end

When the input is specified as above, Q-Chem will determine the MM bond connectivity based on the distances between atoms; if two atoms are sufficiently close, they are considered to be bonded. Occasionally this approach can lead to problems when non-bonded atoms are in close proximity of one another, in which case Q-Chem might classify them as bonded regardless of whether the appropriate MM bond parameters are available. To avoid such a scenario, the user can specify the bonds explicitly by setting the $rem variable USER_CONNECT = TRUE, in which case the $molecule section must have the following format

$molecule
  <Charge> <Multiplicity>
  <Atom> <X> <Y> <Z> <MM atom type> <Bond 1> <Bond 2> <Bond 3> <Bond 4>
  . . .
$end

Each <Bond #> is the index of an atom to which <Atom> is bonded. Four bonds must be specified for each atom, even if that atom is connected to fewer than four other atoms. (For non-existent bonds, use zero as a placeholder.) Currently, Q-Chem supports no more than four MM bonds per atom.

After setting USER_CONNECT = TRUE, a TIP3P water molecule in the AMBER99 force field could be specified as follows:

$molecule
   0 1
   O        -0.790909    1.149780    0.907453  2001  2 3 0 0
   H        -1.628044    1.245320    1.376372  2002  1 0 0 0
   H        -0.669346    1.913705    0.331002  2002  1 0 0 0
$end

Explicitly defining the bonds in this way is highly recommended.

In many cases, all atoms types (within both the QM and MM subsystems) will be defined by a given force field. In certain cases, however, a particular atom type may not be defined in a given force field. For example, a QM/MM calculation on the propoxide anion might consist of a QM subsystem containing an alkoxide functional group, for which MM parameters do not exist. Even though the alkoxide moiety is described using quantum mechanics, van der Waals parameters are nominally required for atoms within the QM subsystem, which interact with the MM atoms via Lennard-Jones-type interactions.

In such cases, there are four possible options, the choice of which is left to the user’s discretion:

1. 1.

   Use a similar MM atom type as a substitute for the missing atom type.

2. 2.

   Ignore the interactions associated with the missing atom type.

3. 3.

   Define a new MM atom type and associated parameters.

4. 4.

   Define a new force field.

These options should be applied with care. Option 1 involves selecting an atom type that closely resembles the undefined MM atom. For example, the oxygen atom of an alkoxide moiety could perhaps use the MM atom type corresponding to the oxygen atom of a neutral hydroxyl group. Alternatively, the atom type could be ignored altogether (option 2) by specifying MM atom type 0 (zero). Setting the atom type to zero should be accompanied with setting all four explicit bond connections to placeholders if USER_CONNECT = TRUE. An atom type of zero will cause all MM energies involving that atom to be zero.

The third option in the list above requires the user to specify a $force_field_params section in the Q-Chem input file. This input section can be used to add new MM atom type definitions to one of Q-Chem’s built-in force fields. At a minimum, the user must specify the atomic charge and two Lennard-Jones parameters (radius and well depth, $ϵ$), for each new MM atom type. Bond, angle, and torsion parameters for stretches, bends, and torsions involving the new atom type may also be specified, if desired. The format for the $force_field_params input section is

$force_field_params
NumAtomTypes <n>
AtomType -1 <Charge> <LJ Radius> <LJ Epsilon>
AtomType -2 <Charge> <LJ Radius> <LJ Epsilon>
. . .
AtomType -n <Charge> <LJ Radius> <LJ Epsilon>
Bond <a> <b> <Force constant> <Equilibrium Distance>
. . .
Angle <a> <b> <c> <Force constant> <Equilibrium Angle>
. . .
Torsion <a> <b> <c> <d> <Force constant> <Phase Angle> <Multiplicity>
. . .
$end

The first line in this input section specifies how many new MM atom types appear in this section (<n>). These are specified on the following lines labeled with the AtomType tag. The atom type numbers are required to be negative and to appear in the order $-1,-2,-3,\mathrm{\dots },-n$. The $molecule section for a water molecule, with user-defined MM parameters for both oxygen and hydrogen, might appear as follows:

$molecule
   0 1
   O        -0.790909    1.149780    0.907453        -1  2 3 0 0
   H        -1.628044    1.245320    1.376372        -2  1 0 0 0
   H        -0.669346    1.913705    0.331002        -2  1 0 0 0
$end

The remainder of each AtomType line in the $force_field_params section consists of a charge (in elementary charge units), a Lennard-Jones radius (in Å), and a Lennard-Jones well depth ($ϵ$, in kcal/mol).

Each (optional) Bond line in the $force_field_params section defines bond-stretching parameters for a bond that contains a new MM atom type. The bond may consist of both atoms <a> and <b> defined an AtomType line, or else <a> may be defined with an AtomType line and <b> defined as a regular atom type for the force field. In the latter case, the label for <b> should be the number of its general van der Waals type. For example, the atom type for a TIP3P oxygen in AMBER99 is 2001, but its van der Waals type is 21, so the latter would be specified in the Bond line. The remaining entries of each Bond line are the harmonic force constant, in kcal/mol/Å${}^2$, and the equilibrium distance, in Å.

Similar to the Bond lines, each (optional) Angle line consists of one or more new atom types along with existing van der Waals types. The central atom of the angle is <b>. The harmonic force constant (in units of kcal/mol/degree) and equilibrium bond angle (in degrees) are the final entries in each Bond line.

Each (optional) Torsion line consists of one or more new MM atom types along with regular van der Waals types. The connectivity of the four atoms that constitute the dihedral angle is <a>–<b>–<c>–<d>, and the torsional potential energy function is

$$E_{\mathrm{torsion}}(\theta )=k_{\mathrm{torsion}}[1+\cos (m\theta -\varphi )]$$

(12.36)

The force constant ($k_{\mathrm{torsion}}$) is specified in kcal/mol and the phase angle ($\varphi$) in degrees. The multiplicity ($m$) is an integer.

Option 4 in the list on page 12.3.2.2 is the most versatile, and allows the user to define a completely new force field. This option is selected by setting FORCE_FIELD = READ, which tells Q-Chem to read force field parameters from a text file whose name is specified in the $force_field_params section as follows:

$force_field_params
    Filename <path/filename>
$end

Here, <path/filename> is the full (absolute) path and name of the file that Q-Chem will attempt to read for the MM force field. For example, if the user has a file named MyForceField.prm that resides in the path /Users/me/parameters/, then this would be specified as

$force_field_params
   Filename /Users/me/parameters/MyForceField.prm
$end

Within the force field file, the user should first declare various rules that the force field will use, including how van der Waals interactions will be treated, scaling of certain interactions, and the type of improper torsion potential. The rules are declared in the file as follows:

RadiusRule <option>
EpsilonRule <option>
RadiusSize <option>
ImptorType <option>
vdw-14-scale <x>
chg-14-scale <x>
torsion-scale <x>

Currently, only a Lennard-Jones potential is available for van der Waals interactions. RadiusRule and EpsilonRule control how to average $\sigma$ and $ϵ$, respectively, between atoms A and B in their Lennard-Jones potential. The options available for both of these rules are Arithmetic [e.g., ${\sigma }_{AB}=({\sigma }_A+{\sigma }_B)/2$] or Geometric [e.g., ${\sigma }_{AB}={({\sigma }_A{\sigma }_B)}^{1/2}$]. RadiusSize has options Radius or Diameter, which specify whether the parameter $\sigma$ is the van der Waals radius or diameter in the Lennard-Jones potential.

ImptorType controls the type of potential to be used for improper torsion (out-of-plane bending) energies, and has two options: Trigonometric or Harmonic. These options are described in more detail below.

The scaling rules takes a floating point argument <x>. The vdw-14-scale and chg-14-scale rules only affect van der Waals and Coulomb interactions, respectively, between atoms that are separated by three consecutive bonds (atoms 1 and 4 in the chain of bonds). These interaction energies will be scaled by <x>. Similarly, torsion-scale scales dihedral angle torsion energies.

After declaring the force field rules, the number of MM atom types and van der Waals types in the force field must be specified using:

NAtom <n>
Nvdw <n>

where <n> is a positive integer.

Next, the atom types, van der Waals types, bonds, angles, dihedral angle torsion, improper torsions, and Urey-Bradley parameters can be declared in the following format:

Atom 1 <Charge> <vdw Type index> <Optional description>
Atom 2 <Charge> <vdw Type index> <Optional description>
. . .
Atom <NAtom> <Charge> <vdw Type index> <Optional description>
. . .
vdw 1 <Sigma> <Epsilon> <Optional description>
vdw 2 <Sigma> <Epsilon> <Optional description>
. . .
vdw <Nvdw> <Sigma> <Epsilon> <Optional description>
. . .
Bond <a> <b> <Force constant> <Equilibrium Distance>
. . .
Angle <a> <b> <c> <Force constant> <Equilibrium Angle>
. . .
Torsion <a> <b> <c> <d> <Force constant 1> <Phase Angle 1> <Multiplicity 1>
. . .
Improper <a> <b> <c> <d> <Force constant> <Equilibrium Angle> <Multiplicity>
. . .
UreyBrad <a> <b> <c>  <Force constant> <Equilibrium Distance>

The parameters provided in the force field parameter file correspond to a basic MM energy functional of the form

$$E_{\mathrm{MM}}=E_{\mathrm{Coul}}+E_{\mathrm{vdW}}+E_{\mathrm{bond}}+E_{\mathrm{angle}}+E_{\mathrm{torsion}}+E_{\mathrm{imptor}}+E_{\mathrm{UreyBrad}}$$

(12.37)

Coulomb and van der Waals interactions are computed for all non-bonded pairs of atoms that are at least three consecutive bonds apart (i.e., 1–4 pairs and more distant pairs). The Coulomb energy between atom types 1 and 2 is simply

$$E_{\mathrm{Coul}}=f_{\mathrm{scale}}\frac{q_1{q}_2}{r_{12}}$$

(12.38)

where $q_1$ and $q_2$ are the respective charges on the atoms (specified with <Charge> in elementary charge units) and $r_{12}$ is the distance between the two atoms. For 1–4 pairs, $f_{scale}$ is defined with chg-14-scale but is unity for all other valid pairs. The van der Waals energy between two atoms with van der Waals types a and b, and separated by distance $r_{ab}$, is given by a “6-12” Lennard-Jones potential:

$$E_{\mathrm{vdW}}(r_{ab})=f_{\mathrm{scale}}{ϵ}_{ab}\left[{\left(\frac{{\sigma }_{ab}}{r_{ab}}\right)}^{12}-2{\left(\frac{{\sigma }_{ab}}{r_{ab}}\right)}^6\right]$$

(12.39)

Here, $f_{\mathrm{scale}}$ is the scaling factor for 1–4 interactions defined with vdw-14-scale and is unity for other valid interactions. The quantities ${ϵ}_{ab}$ and ${\sigma }_{ab}$ are the averages of the parameters of atoms a and b as defined with EpsilonRule and RadiusRule, respectively (see above). The units of <Sigma> are Å , and the units of <Epsilon> are kcal/mol. Hereafter, we refer to atoms’ van der Waals types with a, b, c, ... and atoms’ charges with 1,2, 3, ....

The bond energy is a harmonic potential,

$$E_{\mathrm{bond}}(r_{ab})=k_{\mathrm{bond}}{(r_{ab}-r_{eq})}^2$$

(12.40)

where $k_{\mathrm{bond}}$ is provided by <Force Constant> in kcal/mol/Å${}^2$ and $r_{eq}$ by <Equilibrium Distance> in Å. Note that <a> and <b> in the Bond definition correspond to the van der Waals type indices from the vdw definitions, not the Atom indices.

The bending potential between two adjacent bonds connecting three different atoms (<a>-<b>-<c>) is also taken to be harmonic,

$$E_{\mathrm{angle}}({\theta }_{abc})=k_{\mathrm{angle}}{({\theta }_{abc}-{\theta }_{eq})}^2$$

(12.41)

Here, $k_{\mathrm{angle}}$ is provided by <Force Constant> in kcal/mol/degrees and ${\theta }_{eq}$ by <Equilibrium Angle> in degrees. Again, <a>, <b>, and <c> correspond to van der Waals types defined with vdw.

The energy dependence of the <a>-<b>-<c>-<d> dihedral torsion angle, where <a>, <b>, <c>, and <d> are van der Waals types, is defined by

$$E_{\mathrm{torsion}}({\theta }_{abcd})=f_{\mathrm{scale}}\sum _m{k}_{abcd}[1+\text{cos}(m{\theta }_{abcd}-\varphi )]$$

(12.42)

Here, $f_{scale}$ is the scaling factor defined by torsion-scale. The force constant $k_{abcd}$ is defined with <Force constant> in kcal/mol, and the phase angle $\varphi$ is defined with <Phase Angle> in degrees. The summation is over multiplicities, $m$, and Q-Chem supports up to three different values of $m$ per dihedral angle. The force constants and phase angles may depend on $m$, so if more than one multiplicity is used, then <Force constant> <Phase Angle> <Multiplicity> should be specified for each multiplicity. For example, to specify a dihedral torsion between van der Waals types 2–1–1–2, with multiplicities $m=2$ and $m=3$, we might have:

   Torsion  2  1  1  2  2.500  180.0  2 1.500 60.0 3

Improper torsion angle energies for four atoms <a>-<b>-<c>-<d>, where <c> is the central atom, can be computed in one of two ways, as controlled by the ImptorType rule. If ImptorType is set to Trigonometric, then the improper torsion energy has a functional form similar to that used for dihedral angle torsions:

$$E_{\mathrm{imptor}}({\theta }_{abcd})=\frac{k_{abcd}}{N_{\mathrm{equiv}}}[1+\text{cos}(m{\theta }_{abcd}-\varphi )]$$

(12.43)

Here, ${\theta }_{abcd}$ is the out-of-plane angle of atom <c>, in degrees, and $k_{abcd}$ is the force constant defined with <Force Constant>, in kcal/mol. The phase $\varphi$ and multiplicity $m$ need to be specified in the Improper declaration, although the definition of an improper torsion suggests that these values should be set to $\varphi =0$ and $m=2$. The quantity $N_{equiv}$ accounts for the number of equivalent permutations of atoms <a>, <b>, and <d>, so that the improper torsion angle is only computed once. If ImptorType is set to Harmonic, then in place of Eq. (12.43), the following energy function is used:

$$E_{\mathrm{imptor}}({\theta }_{abcd})=\frac{k_{abcd}}{N_{\mathrm{equiv}}}{\theta }_{abcd}^2$$

(12.44)

The Urey-Bradley energy, which accounts for a non-bonded interaction between atoms <a> and <c> that are separated by two bonds (i.e., a 1-3 interaction through <a>-<b>-<c>), is given by

$$E_{\mathrm{UreyBrad}}(r_{ac})=k_{abc}{(r_{ac}-r_{eq})}^2$$

(12.45)

The distance in Å between atoms <a> and <c> is $r_{ac}$, the equilibrium distance $r_{eq}$ is provided by <Equilibrium Distance> in Å, and the force constant $k_{abc}$ is provided by <Force Constant> in kcal/mol/Å${}^2$.

A short example of a valid text-only file defining a force field for a flexible TIP3P water could be as follows:

//-- Force Field Example --//

// -- Rules -- //
RadiusRule Geometric
RadiusSize Radius
EpsilonRule Geometric
ImptorType Trigonometric
vdw-14-scale 1.0
chg-14-scale 0.8
torsion-scale 0.5

// -- Number of atoms and vdw to expect -- //
NAtom 2
Nvdw 2

// -- Atoms -- //
Atom 1 -0.8340  2  TIP3P Oxygen
Atom 2  0.4170  1  TIP3P Hydrogen

// -- vdw -- //
vdw 1  0.0000  0.0000 H parameters
vdw 2  1.7682  0.1521 O parameters

// -- Bond -- //
Bond 1 2  553.0  0.9572

// -- Angle -- //
Angle 1 2 1  100.0  104.52

Lines that do not begin with one of the keywords will be ignored, and have been used here as comments.

For QM/MM calculations (but not for purely MM calculations) the user must specify the QM subsystem using a $qm_atoms input section, which assumes the following format:

$qm_atoms
   <QM atom 1 index> <QM atom 2 index> . . .
   . . .
   <QM atom n index>
$end

Multiple indices can appear on a single line and the input can be split across multiple lines. Each index is an integer corresponding to one of the atoms in the $molecule section, beginning at 1 for the first atom in the $molecule section. Link atoms for the ONIOM model and YinYang atoms for the Janus model are not specified in the $qm_atoms section, as these are inserted automatically whenever a bond connects a QM atom and an MM atom.

Q-Chem 4.2.2 and later versions also support, for example

$qm_atoms
   18:31 35
$end

which specifies 15 QM atoms (atoms 18 through 31; atom 35).

For Janus QM/MM calculations, there are several ways of dealing with van der Waals interactions between the QM and MM atoms. By default, van der Waals interactions are computed for all QM–MM and MM–MM atom pairs but not for QM–QM atom pairs. In some cases, the user may prefer not to neglect the van der Waals interactions between QM–QM atoms, or the user may prefer to neglect any van der Waals interaction that involves a QM atom. Q-Chem allows the user this control via two options in the $forceman section. To turn on QM–QM atom van der Waals interactions, the user should include the following in their input:

$forceman
   QM-QMvdw
$end

Similarly, to turn off all van der Waals interactions with QM atoms, the following should be included:

$forceman
   NoQM-QMorQM-MMvdw
$end

Periodic boundary conditions (using Ewald summation for the long-range Coulomb interactions) can be used in conjunction with both MM-only calculations and QM/MM calculations. The approach is based off of the work of Nam et al.⁶⁶⁰ and (independently) Riccardi et al.,⁷⁸⁹ as implemented in both the Amber⁹⁶⁶ and Charmm^{789, 122} programs. These approaches use Mulliken charges to represent the periodic images of the wave function, and while suitable for semi-empirical calculations with minimal basis sets, instabilities in the Mulliken charges for extended basis sets lead to SCF convergence failure in the QM/MM-Ewald calculations.³⁹⁰ The implementation in Q-Chem thus allows for the use of ChElPG charges to represent the image wave functions, affording an algorithm that is stable in extended basis sets.^{390, 389}

The efficiency of the Ewald summation is governed by the parameter, $\alpha$, that controls the partition of the Coulomb potential into short- and long-range components, and in the QM/MM-Ewald method there are separate values of $\alpha$ for the QM and MM portions of the calculations. Improper selection of ${\alpha }_{\mathrm{MM}}$ and/or ${\alpha }_{\mathrm{QM}}$ can greatly increase the computational time, and the choices that are optimal for MM calculations need not be optimal for QM/MM calculations.³⁹⁰ The cost of the MM Ewald summation scales as $𝒪(N_{\mathrm{recip}}{N}_{\mathrm{atoms}})$, where $N_{\mathrm{recip}}$ is the number of reciprocal-space lattice vectors that is used for the $k$-space sum. The QM portion of the calculation scales as $𝒪(N_{\mathrm{recip}}{N}_{\mathrm{QM}}{N}_{\mathrm{atoms}})$, where $N_{\mathrm{QM}}$ is the number of QM atoms (whereas $N_{\mathrm{atoms}}=N_{\mathrm{QM}}+N_{\mathrm{MM}}$ is the total number of atoms). The MM Ewald parameter is thus selected to minimize the amount of work that is done in real space. The optimal value, which is typically ${\alpha }_{\mathrm{MM}}\approx 0.5{\text{Å}}^{-1}$, can be found by solving the equation³⁹⁰

$${\alpha }_{\mathrm{MM}}=2C/L$$

(12.46)

where

$$C={\left[-\ln \left({10}^{-\text{SCF\_CONVERGENCE}}\right)\right]}^{1/2}$$

(12.47)

and $L$ is the length of the simulation cell. (Only cubic simulation cells are available at present.) In contrast, the parameter ${\alpha }_{\mathrm{QM}}$ should be selected to minimize the total number of vectors in both real and reciprocal space. The optimal value, which is often ${\alpha }_{\mathrm{QM}}\approx 0.1{\text{Å}}^{-1}$, is determined by solving the equation³⁹⁰

$$\frac{2C{L}^3{\alpha }_{\mathrm{QM}}^3}{{\pi }^{3/2}}+\frac{{\alpha }_{\mathrm{QM}}^2{L}^2}{{\pi }^{1/2}}-{\alpha }_{\mathrm{QM}}L-2C=0.$$

(12.48)

To perform an MM- or QM/MM-Ewald job, one must set MM_SUBTRACTIVE = TRUE and EWALD_ON = TRUE in the $rem section, but otherwise job control is largely done through the $forceman section. The following variables must be set for every type of Ewald calculation.

• •

  The keyword Ewald will turn on Ewald summation.

• •

  The keyword alpha should be followed by a value for the MM Ewald parameter and then the QM Ewald parameter. (The latter must be set even for MM-only jobs.)

• •

  Box_length specifies the side length of the cubic simulation cell, in Å.

The following parameters are optional for further job control.

• •

  Dielectric specifies a dielectric constant for the surrounding medium, which appears in the “dipole term” of Ewald summation ($E_{\mathrm{dipole}}$ in Ref. 390). If no value is set, the dielectric constant is set to infinity, corresponding to “tin foil" boundary conditions.

• •

  The keyword Ewald_SCF_thresh_on, followed by a real number, causes Q-Chem to wait until the DIIS error falls below the specified value before adding the Ewald correction to the Fock matrix, thus obviating the sometimes-costly Ewald correction in early SCF cycles. (The default value is 1.0, which turns on Ewald summation immediately in most cases)

A short example of a $forceman section using Ewald summation could be as follows:

$forceman
   ewald
   alpha                 0.35   0.1
   box_length            15.00
   dielectric            88.0
   mm_read_scratch
   ewald_scf_thresh_on   0.0001
$end