Covalent Bonds Divided by the QM/MM Barrier

Hybrid QM/MM simulations of biomolecular systems often present situations where only a portion of a molecule should be treated quantum mechanically, usually to save computational resources since the cost of simulating QM regions rises rapidly with the number of simulated toms. In order to deal with chemical bonds that are split by the QM/MM division of the biomolecular system, that is, bonds that have one atom in the quantum (QM) region and another in the classical (MM) region (we will call these ``QM/MM bonds"), NAMD makes approximations to the molecular system in order to bridge differences in simulation type (QM vs. MM), and minimize errors involved in the QM/MM division of the system (Figure 17 A and B).

Figure 17: Treatment of QM/MM bonds. A) Illustration of all atoms in the vicinity of the QM/MM bond, colored by element: cyan for carbon, white for hydrogen, blue for nitrogen and red for oxygen. B) To the left, in blue, is the region that will be treated with the chosen QC software. To the right, in red, the region treated classically by NAMD. The bond in gray is the one crossing the QM/MM division. The atom marked as QM1 is the quantum atom directly connected to the classical atom on the other side of the QM/MM division. Analogously, the atom marked as MM1 is the classical atom directly connected to the quantum atom on the other side of the QM/MM division. Atoms marked as MM2 are directly bonded to the MM1 atom, and atoms marked MM3 are directly bonded to MM2 atoms. C) Z1 method. Ignored partial charges are indicated in the image with a gray sphere taking the place of its respective classical atom. Directly adjacent to MM1 is a green sphere representing the link atom that is placed along the QM1-MM1 covalent bond. All remaining partial charges representing classical atoms that are passed on to the QC software are indicated in purple spheres. D) Z2 method. E) Z3 method. F) RCD method. Virtual point charges, are represented in yellow spheres. The text indicates the total charge placed at each position, where q indicates the charge of the MM1 atom and q2 represents the partial charge of the MM2 atom at that position. The yellow spheres at MM2 atom positions indicate their partial charge has been changed from its original value. G) CS method. Since in this case the virtual point charges are placed very close to the MM2 atom position, the yellow spheres representing them show significant overlapping.

Link Atoms

As previously mentioned, the information regarding atoms in the QM region is passed on to the chosen QC software, that is, their respective positions and element types, but in order to maintain (or approximate) the effect of the chemical bond between the QM atom and the MM atom, NAMD creates and places a link atom (usually a hydrogen) along the ``broken" QM/MM bond. The user can fine-tune this process by choosing the method of placement of the link atom and even the element of such atom (keywords ``qmBondDist" and ``qmLinkElement").

The introduction of the link atom will invariably place it very near the classical atom involved in the QM/MM bond, therefore the use and placement of partial charges from classical atoms becomes highly relevant. Under the mechanical embedding scheme, the QC software only receives the atoms in the QM region and the link atoms created to approximate QM/MM bonds, so no manipulation of partial charges is required. On the other hand, usual QM/MM simulations are done under the electrostatic embedding scheme, in which case the partial charges of classical atoms involved in the QM/MM bonds and classical atoms directly connected to them require special treatment.

Point Charge Alterations

Several methods have been proposed to handle this situation, and the QM/MM interface developed here implements the most widely accepted ones. One can be chosen using the ``qmBondScheme" keyword (Figure 17 C to G). In all implemented methods, the classical atom participating in the QM/MM bond (MM1 atom) does not have its partial charge passed on to the QC software, since this would create excessive repulsion (or attraction) on the link atom. This is, in fact, the entirety of the ``Z1" method: ignoring the partial charge of the MM1 atom. Analogously, ``Z2" and ``Z3" ignore all partial charges up to MM2 and MM3 atoms, respectively (Figure 17 C to E).

The Redistributed Charge and Dipole (RCD) method (Figure 17 F) is more elaborate, as it rearranges the partial charge of the MM1 atom (indicated as $q$ ) so that the total charge of the region is maintained as well as the dipole moments of the bonds between MM1 and MM2 atoms. This is done by creating ``virtual" point charges, which are passed on to the QC software as if they represented partial charges of classical atoms. More specifically, the RCD method creates a virtual point charge in the middle of all MM1-MM2 bonds with a charge of $2q/n$ , where $n$ is the number of MM2 atoms connected to MM1, and also subtracts a charge $q/n$ from each of the MM2 atoms, so that the total charge of the region remains constant while approximating the dipole moment of the MM1-MM2 bonds. This way there will be no point charge placed at the position of the MM1 atom, but its partial charge is not simply removed, it is redistributed.

A similar approach is taken by the Charge Shifting (CS) method (Figure 17 G). In this case, the MM1 partial charge is equally distributed across the MM2 atoms, and two virtual point charges are placed along the direction of the MM1-MM2 bond, one before the MM2 atom and one after, each one with a charge of $+q/n$ and $-q/n$ respectively. This method will also keep the total charge of the region constant while trying to preserve the local dipoles formed by all MM1-MM2 bonds.

Link Atom Charge and Charge Groups

Along with the gradient over all QM atoms, NAMD can also use the partial charge derived from the QC calculation to update the charge distribution of atoms in the QM region. When a QM/MM bond exists, however, part of the charge of the region will be placed on the link atom, and in order to keep the charge of the QM region constant, the link atom charge is re-distributed on the QM region. This seemingly simple mechanism can cause problems unless special care is be taken when deciding which bond will mark the division of QM and MM regions.

Many force fields divide the topologies of biomolecules in ``charge groups" (Figure 18 A and B). What this means is that not only will the partial charges of all atoms of a molecule add up to the whole number that represents the charge of the molecule, they will also add up to whole numbers in sub groups of atoms (look for the ``GROUP" statements in http://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/node24.htmlto see an example). Therefore, one needs to make sure that the chosen QM/MM bond(s) sits in between ``charge groups", so the total sum of partial charges of atoms defining a QM region is a whole number. This is especially important in order to keep the total charge of the system constant. Since the QC calculation will always distribute a whole charge over all atoms of a system (QM atoms plus a link atom), if the partial charge of QM atoms is not initially a whole number, it will be forced into a whole number after the first QC step, where the charge of the link atom is distributed over the QM region. This will create a mismatch between QM and MM charges, changing the total charge of the entire system (QM plus MM regions).

Figure 18: Charge Groups and QM/MM Bonds. A) Illustration of aspartate and the distribution of charge over its atoms as in CHARMM36 force field parameters. Circles in red indicate oxygen atoms, blue indicate nitrogen atoms, cyan for carbon atoms, and white for hydrogen atoms. ``Bond 1" indicates the peptide bond, ``Bond 2" indicates the one between the alpha carbon and the peptide bond nitrogen, and ``Bond 3" the bond between the alpha carbon and the peptide bond carbon. B) Charge groups are indicated with dashed squares, along with their total charges. C) Depiction of the atoms in the QM region if Bonds 1 and 3 are used to separate it from the MM region. The total charge of QM region is $-1$ . D) Depiction of QM region if the same is defined by Bonds 2 and 3. In this case, the total charge of QM region is $-1.16$ .

An example can be seen in Figure 18, bonds 1 and 3 are chosen as the QM/MM bonds, the charge distribution seen in Figure 18 C shows a whole charge for the QM region (and consequently for the MM region). Therefore, any charge placed on link atoms can be redistributed to the QM atoms with no change in total system charge. However, if bonds 2 and 3 are chosen for the QM/MM bond (Figure 18 D), the charge of the MM region would be $+1.16$ , while the charge of the QM region would be $-1.16$ . Since the QC calculation would place a pre-determined whole charge on the region ($-1$ , in this case), the updated charge of the QM region after the first simulation step would change the total charge of the system to $+0.16$ , in this example.