Development and Application of QM/MM Methods with Advanced Polarizable Potentials

Jorge Nochebuena Sehr Naseem-Khan Andrés Cisneros

Abstract

Quantum Mechanics/Molecular Mechanics (QM/MM) simulations are a popular approach to study various features of large systems. A common application of QM/MM calculations is in the investigation of reaction mechanisms in condensed–phase and biological systems. The combination of QM and MM methods to represent a system gives rise to several challenges that need to be addressed. The increase in computational speed has allowed the expanded use of more complicated and accurate methods for both QM and MM simulations. Here, we review some approaches that address several common challenges encountered in QM/MM simulations with advanced polarizable potentials, from methods to account for boundary across covalent bonds and long–range effects, to polarization and advanced embedding potentials.

1 Introduction

The simulation of reaction processes in complex biological systems such as enzymes, or condensed–phase systems, requires the description of changes in electronic and nuclear degrees of freedom. Although quantum methods can adequately describe electronic effects, the modeling of the systems of interest is limited by the size and the time scale that can be carried out in reasonable computational times. By contrast, molecular mechanics methods allow the study of systems with millions atoms but are limited to the type of system for which the force field was designed, and most force fields used for these simulations are non–reactive [123].

Quantum mechanics/molecular mechanics (QM/MM) methods emerged in the mid-1970s as a proposal to study enzymatic reactions [119]. Since then, QM/MM methods have been widely applied to study various chemical and biochemical systems [79, 108, 114, 16, 32, 3, 87]. QM/MM methods aim to leverage the advantages of both quantum and classical approaches by combining two levels of theory. The quantum model focuses on the region of the system where the reactive processes of interest take place, while the classical model describes the rest of the system. For example, if the process involves a reaction in the active site of an enzyme (see Figure 1), the residues and any other molecules (or fragments) involved in the reaction will be assigned as the QM subsystem, and the rest of the enzyme, solvent, counterions, etc. will be represented by the MM force field.

Refer to caption — Figure 1: Main protease of SARS-CoV2 in a box of water [73] The whole protein (left) and a close-up in the active site (right) are shown. Histidine 41 and cysteine 145 amino acid residues are represented by balls and cylinders. PDB: 6LU7 [66].

There are a wide variety of combinations for the treatment of the QM and MM subsystem. The QM subsystem may be represented by empirical valence bond (EVB) [120, 7, 72], semi–empirical [56, 95] or ab initio Hamiltonians [110, 65, 63, 41]. For the MM environment many QM/MM simulations employ non–polarizable point–charge force fields (npFFs) [15, 70, 19]. More recently, increase in computational speed and algorithmic developments have allowed the use of more accurate potentials which may include explicit representation of electronic polarization, more elaborate electrostatics, and in some cases, explicit inclusion of other quantum effects [58, 116, 51, 52].

The coupling of QM and MM methods to represent a system gives rise to several challenges such as how to couple these two levels of theory across covalent boundaries, how to treat long–range effects in the context of QM subsystems, coupling of the quantum and classical Hamiltonians and how to treat the explicit couplings, among others. In this review, we present approaches that tackle several of these issues.

This article is organized as follows. Section 2 provides an overview of the QM/MM approach. It describes what force fields are, the different types of embedding, as well as some methods to include polarization. Section 3 describes the challenges involved in separating the QM and MM regions across covalent bonds, and solutions that have been proposed. The differences between the link atom, double link atom, as well as LSCF and GHO methods will be described. In subsection 3.2 the role played by long-range interactions and methods that have been developed to include them efficiently will be described. In section 4 we will describe some recent development for QM/MM simulations using advanced polarizable force fields. Finally, in section 5 we will show some of the methods implemented in LICHEM, a code developed in the group to perform QM/MM simulations, followed by concluding remarks.

2 Classical MM environments in QM/MM

Many QM/MM simulations rely on the combination of QM methods with the use of potential-energy functions, called force fields [94]. Many force fields approximate the energy and forces of the system by separating the contributions into bonded and non–bonded interactions. In several cases, the non–bonded interactions are approximated by Coulomb and Van der Waals interactions. Some examples are: AMBER [19], CHARMM [15], GROMOS [107], MMFF [59], MM3 [5], MM4 [4], OPLS [70] and UFF [103].

A common expression for fixed–charge, non–polarizable force fields is show in Eq. (1)

	$\displaystyle E_{\textnormal{int}}=$	$\displaystyle\sum_{\textnormal{bonds}}K_{b}(b-b_{0})^{2}+\sum_{\textnormal{angles}}K_{\theta}(\theta-\theta_{0})^{2}+\sum_{\textnormal{dihedrals}}K_{\phi}[1+\cos(n\phi+\delta)]$
		$\displaystyle+\sum_{i<j}\frac{q_{i}q_{j}}{r_{ij}}+\sum_{i<j}4\varepsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right],$		(1)

where the first three terms are the bonded contributions and the last two represent non-bonded interactions. $K_{b}$ , $K_{\theta}$ and $K_{\phi}$ are force constants for bonds, angles and dihedral terms; $b_{0}$ and $\theta_{0}$ are equilibrium values for bond length and angles between atoms; $n$ is the dihedral multiplicity; $\delta$ is the dihedral angle phase; $q$ are partial (generally atomic) charges; $\varepsilon$ are the well depths, and $\sigma$ are the van der Waals radii. In the above equation, and throughout the rest of the review all equations use atomic units.

Given the functional form of the potential, the total energy for a QM/MM system can be separated in 3 contributions: the contribution from the QM subsystem, the contribution from the MM environment, and the interaction between both regions. The first two contributions are straightforward to evaluate. In principle, any combination of QM and MM methods can be chosen. The interaction between the QM and MM regions poses some challenges as mentioned above. The combination of the levels of theory can be achieved by two main approaches termed the subtractive and additive coupling schemes.

In the subtractive method the QM/MM energy is obtained by:

E_{\textnormal{QM/MM}}^{\textnormal{sub}}=E_{\textnormal{MM}}(\textnormal{MM}+\textnormal{QM})+E_{\textnormal{QM}}(\textnormal{QM})-E_{\textnormal{MM}}(\textnormal{QM}),

(2)

that is, first the energy of the complete system (QM and MM regions) is evaluated at the MM level, then the energy of the QM region calculated in the QM level is added, and finally, the energy of the QM region evaluated at the MM level is subtracted. This method is straightforward but has the disadvantage that it requires several calculations at different levels of theory of the same set of atoms.

In the additive scheme, the QM/MM energy is calculated by:

E_{\textnormal{QM/MM}}^{\textnormal{add}}=E_{\textnormal{QM}}(\textnormal{QM})+V_{\textnormal{MM}}(\textnormal{MM})+E_{\textnormal{QM/MM}}(\textnormal{QM}+\textnormal{MM}),

(3)

that is, the energy QM/MM is the sum of the contributions QM and MM regions evaluated at its own level plus the contribution due to coupling. The complexity of the QM/MM methods depends on how the coupling term is expressed.

The last term in Eq. (3) can be further subdivided depending on the individual terms of the classical force field. If the force field uses a functional form as in Eq. (1), the QM/MM interaction can be divided as:

E_{\textnormal{QM/MM}}^{\textnormal{total}}=E_{\textnormal{QM/MM}}^{\textnormal{bonds}}+E_{\textnormal{QM/MM}}^{\textnormal{angles}}+E_{\textnormal{QM/MM}}^{\textnormal{torsions}}+E_{\textnormal{QM/MM}}^{\textnormal{Coulomb}}+E_{\textnormal{QM/MM}}^{\textnormal{VDW}}.

(4)

In both substractive and attractive schemes, it is possible to further subdivide the approaches depending on how each of the terms in the QM/MM interaction is calculated. In the mechanical embedding approach, the interaction between the QM and MM regions are handled at the force field level. The QM subsystem is thus kept in place by MM interactions. Similarly to Eq. (1), all bonded terms involving QM atoms are represented with the respective MM functions, and all non–bonded interactions involving QM and MM atoms are represented by point charges and Lennard–Jones potential energy terms. In this approach, the calculation of the QM system is performed isolated from the MM region, and therefore, the wave function is not explicitly polarized by the external field from the MM environment.

In the electrostatic embedding scheme, the electrostatic interactions between the two regions, $E_{\textnormal{QM/MM}}^{\textnormal{Coulomb}}$ , are included in the QM calculation. In this case, the MM atoms polarize the electrons in the QM subsystem via an effective Hamiltonian,

H_{\textnormal{{eff}}}=H_{\textnormal{{QM}}}-\sum_{j}^{M}\frac{q_{j}}{|\textbf{r}_{i}-\textbf{R}_{j}|},

(5)

where $H_{\textnormal{QM}}$ is the electronic Hamiltonian for the QM subsystem, $M$ is the number of MM atoms that have a partial charge $q_{j}$ , and $\textbf{r}_{i}$ and $\textbf{R}_{j}$ are the positions of electron $i$ and MM atom $j$ . In this approach, the QM subsystem interacts with the MM region as a set of external point charges (or higher order multipoles). However, mutual polarization effects or other non-electrostatic interactions are not taken into account.

2.1 Polarizable QM/MM

The force fields mentioned so far are extremely efficient and provide a convenient approach to investigate myriad systems. These force fields in some cases can exhibit limited accuracy in the reproduction of certain regions of the potential energy surface. One major reason for the reduced accuracy involves the approximations employed to represent the bonded and non–bonded interactions. In many cases, the non–bonded interactions are represented by two separate terms assuming pair–wise additivity. These terms aim to represent the underlying physical interactions, however there are several approximations that are employed, which lead to the neglect of some physical interactions. In general, the electrostatic interactions between particles are approximated by fixed atomic partial charges. Dispersion and repulsion effects are approximated by a Lennard–Jones (or similar) function [92, 115, 13].

The use of fixed atomic point charges may be unable to capture or accurately represent several important effects such as charge density anisotropy, and charge penetration. In addition, the approximation of the charge density by fixed point charges neglects certain many-body effects such as explicitly accounting for electronic polarization or charge–transfer. Polarization effects can be included by explicitly taking into account the effect of the electric field on the charge distribution of the constituent molecules. Several polarizable force fields have been developed including AMOEBA [101], CHARMM [85], EFP [55], GEM [24, 99], GROMOS [48], NEMO [61] and SIBFA [57], which employ one of the methods below to describe polarization effects.

In the case where the force field has an explicit term to represent the electronic polarization of the classical environment, it is necessary to account for the interaction between this effect in the MM and the QM subsystem. In this case, the polarization embedding scheme is such that the QM and MM regions are mutually polarized. Thus, not only the QM region is polarized by the QM charges, but also the MM region responds to the chemical environment in the QM region. Three basic types of polarizable models have been developed: fluctuating charge model [80], Drude oscillator model [14, 14], and induced dipole model [84, 12].

The fluctuating charge (FQ) model is one of the most straightforward polarizable models. It is based on the principle of electronegativity equalization. That is, a charge oscillates between the atoms until the electronegativities of the atoms are equalized. This method allows the modification of the value of the partial charges in response to the electric field, thus altering molecular polarization. This can be done by coupling the charges to the environment using electronegativity or chemical potential equalization schemes. FQ method uses a Taylor expansion truncated up to second-order terms of the atomic chemical potential to express the energy required to create a charge, Q,

E(Q)=E^{0}+\left(\frac{\partial E}{\partial Q}\right)_{0}Q+\frac{1}{2}\left(\frac{\partial^{2}E}{\partial Q^{2}}\right)_{0}Q^{2}

(6)

The energies required to create +1 and -1 charges on an atom, obtained from Eq. (6)

	$\displaystyle E(+1)$	$\displaystyle=E^{0}+\left(\frac{\partial E}{\partial Q}\right)_{0}+\frac{1}{2}\left(\frac{\partial^{2}E}{\partial Q^{2}}\right)_{0}$		(7)
	$\displaystyle E(-1)$	$\displaystyle=E^{0}-\left(\frac{\partial E}{\partial Q}\right)_{0}+\frac{1}{2}\left(\frac{\partial^{2}E}{\partial Q^{2}}\right)_{0}$		(8)

using the definition of ionization potential as $IP=E(+1)-E(0)$ and electron affinity as $EA=E(0)-E(-1)$ we obtain:

	$\displaystyle\left(\frac{\partial E}{\partial Q}\right)_{0}$	$\displaystyle=\frac{1}{2}(IP+EA)=\chi^{0}$		(9)
	$\displaystyle\left(\frac{\partial^{2}E}{\partial Q^{2}}\right)_{0}$	$\displaystyle=IP-EA=J^{0}=2\eta^{0}$		(10)

where the $\chi^{0}$ is the electronegativity and $\eta^{0}$ is the chemical hardness. Equation (6) can be rewritten as:

E(Q)=E^{0}+\chi^{0}Q+\eta^{0}Q^{2}

(11)

The previous equation can be generalized for a set of $N$ atoms

E(Q_{1}...Q_{N})=\sum_{i}\left(E_{i0}+\chi_{i}^{0}Q_{i}+n_{ii}^{0}\right)+\sum_{i,j>i}2\eta_{ij}Q_{i}Q_{j}

(12)

The optimal charge distribution is achieved by minimizing the energy with respect to the charges on each atom:

\frac{\partial E}{\partial Q_{i}}=0,\hskip 19.91684pti=1,N_{i}

(13)

The constraints used in solving the previous equation include total charge conservation:

Q_{\textnormal{tot}}=\sum_{i=1}^{N}Q_{i}

(14)

and electronegativity equalization at convergence:

\chi_{1}=\chi_{2}...=\chi_{N}

(15)

The Drude oscillator model incorporates electronic polarizability by representing an atom as a two-particle system: a charged core with charge $q_{0}$ and a charged shell, also called a Drude particle, with charge $q_{D}$ . The core and shell are linked by a harmonic spring. In the absence of an electric field, the Drude particle is located at the atomic core. Otherwise, the Drude particle will be at a distance $d$ from the atomic core in the presence of an electric field:

d=\frac{q_{D}\textnormal{E}}{k_{D}}

(16)

The atomic induced dipole then is treated as:

\mu=q_{D}d=q_{D}\times\frac{q_{D}E}{k_{D}}=\frac{q_{D}^{2}E}{k_{D}}

(17)

The atomic polarizability, $\alpha$ , is related to the force constant $k$ of the harmonic spring connecting the core and shell and is determined by

\alpha=\frac{q_{D}^{2}}{k_{D}}

(18)

In principle, the Drude oscillator, $q_{D}$ , the force constant of the harmonic spring, $k_{D}$ , or both can be tuned to achieve an appropriate polarization response.

In the induced dipole model the polarization response is represented by a set of inducible dipoles that arise due to the external permanent and induced electric field. The induced dipole moment at the site $i$ is:

\mu_{i}=\alpha_{i}\left[E_{i}-\sum_{j\neq i}^{N}T_{ij}\mu_{j}\right]

(19)

where $T_{ij}$ is the dipole–dipole interaction tensor defined by:

T_{ij}=\frac{1}{r_{ij}^{3}}I-\frac{3}{r_{ij}^{5}}\left[\begin{matrix}x^{2}&xy&xz\\ yx&y^{2}&yz\\ zx&zy&z^{2}\end{matrix}\right]

(20)

where $I$ is the identity matrix and $x$ , $y$ and $z$ are Cartesian components along the vector between atoms $i$ and $j$ at distance $r_{ij}$ .

When polarizable potentials are employed in QM/MM calculations, an additional energy term arises in the QM/MM interaction: $E_{\textnormal{QM/MM}}^{\textnormal{Pol}}$ , which is added to the total energy:

E_{\textnormal{QM/MM}}^{\textnormal{total}}=E_{\textnormal{QM/MM}}^{\textnormal{Bonded}}+E_{\textnormal{QM/MM}}^{\textnormal{Coulomb}}+E_{\textnormal{QM/MM}}^{\textnormal{VanderWaals}}+E_{\textnormal{QM/MM}}^{\textnormal{Pol}}

(21)

where $E_{\textnormal{QM/MM}}^{\textnormal{Pol}}$ describes the explicit polarization, in the case of dipolar force fields, arising from the dipolar interaction between the induced dipoles and the permanent and induced electric field. This term is calculated via the following formula:

E_{\textnormal{QM/MM}}^{\textnormal{Pol}}=-\frac{1}{2}\sum_{i}\alpha_{i}E_{i}^{0}E_{i}

(22)

where $E_{i}^{0}$ represents the electrostatic field on atom $i$ arising from the permanent charges, and higher permanent moments (if present), and $E_{i}$ represents the sum of $E_{i}^{0}$ and the electrostatic field on atom $i$ due to the induced dipoles on other sites.

AMOEBA is an example of a polarizable force field that employs the induced dipole formalism and has been implemented for QM/MM simulations. [74, 37, 88, 82, 54, 36] AMOEBA uses the following functional form:

$\displaystyle E_{\textnormal{int}}^{\textnormal{AMOEBA}}(\mathbf{r}^{N})$	$\displaystyle=$	$\displaystyle\sum_{\textnormal{bonds}}k_{i}^{b}(l_{i}-l_{i,0})^{2}[1+2.55(l_{i}-l_{i,0})+3.793125(l_{i}-l_{i,0})^{2}]$
	$\displaystyle+$	$\displaystyle\sum_{\textnormal{angles}}k_{i}^{\theta}(\theta_{i}-\theta_{i,0})^{2}[1+0.014(\theta_{i}-\theta_{i,0})+5.6\times 10^{-5}(\theta_{i}-\theta_{i,0})^{2}$
	$\displaystyle+$	$\displaystyle 7.0\times 10^{-7}(\theta_{i}-\theta_{i,0})^{3}+2.2\times 10^{-8}(\theta_{i}-\theta_{i,0})^{4}]$
	$\displaystyle+$	$\displaystyle\sum_{\textnormal{torsions}}\frac{V_{n}}{2}(1+\cos(n\omega-\gamma))+\sum_{\textnormal{oop}}0.02191414k_{\chi}\chi^{2}$
	$\displaystyle+$	$\displaystyle\sum_{\textnormal{str-bend}}k_{sb}(b_{i}-b_{i,0})(\theta_{i}-\theta_{i,0})+\sum_{\textnormal{PI-tor}}U_{\textnormal{PI-tor},i}+\sum_{\textnormal{tor-tor}}U_{\textnormal{tor-tor},i}$
	$\displaystyle+$	$\displaystyle\frac{1}{2}\left[\sum_{i<j}^{N}M_{i}T_{ij}M_{j}+\sum_{i<j}^{N}\epsilon_{ij}\left(\frac{1+\delta}{\rho_{ij}+\delta}\right)^{n-m}\left(\frac{1+\gamma}{\rho_{ij}^{m}+\gamma}-2\right)\right]+\frac{1}{2}\sum_{i}^{N}\mu_{i}^{ind}E_{i}$

Given that the induced dipoles depend on the electric fields, Eq. (22) can be re–written as:

E_{\textnormal{QM/MM}}^{\textnormal{Pol}}=\frac{1}{2}\sum_{i}\sum_{j\neq i}\mu_{i}T_{ij}\mu_{j}-\sum_{i}\mu_{i}(E_{i}^{\textnormal{MM}}+E_{i}^{\textnormal{QM}})

(23)

where the first term corresponds to the induced dipole interaction, and the second corresponds to the interaction of the induced dipoles at site $i$ with the electric field at site $I$ arising from the MM ( $E_{i}^{\textnormal{MM}}$ ) and QM ( $E_{i}^{\textnormal{QM}}$ ) subsystems.

It is important to note that the equations for the induced dipoles and the QM wave function are coupled. These equations can be solved iteratively using an SCF approach by computing the induced dipoles at each SCF cycle, resulting in a fully self–consistent (fsc) approach [84]. Alternatively, it is possible to de–couple the solution of the iterative dipoles by approximating the MM external field interacting with the QM wavefunction by only including the permanent field arising from the MM environment [74], resulting in a partially self–consistent (psc) QM/MM polarization approach.

Force fields that explicitly consider polarization changes in response to the chemical environment offer a more accurate description of electrostatic interactions. In fact, it is known that polarization effects become more relevant in highly charged systems [78]. The dipole moment in a molecule is known to change significantly according to the state of aggregation. A polarizable force field has the versatility to reproduce both gas phase and liquid phase properties simultaneously. For example, the non-polarizable TIP3P [69] additive model calculates a hydrogen bond energy for a water dimer greater than the accepted value. In contrast, polarizable models produce energy closer to the values calculated with ab initio methods [86]. In addition, distributed multipoles introduce charge density anisotropy that is lost in isotropic point charge distributions. Table 1 shows the mean absolute error for 10 water dimers calculated with TIP3P and AMOEBA [74]. The results show that the AMOEBA polarizable force field describes hydrogen bonds better than the TIP3P model.

Table 1: Mean absolute error of the binding energy for 10 water dimers compared with BSSE corrected binding energies from PBE0/aug- cc-pVTZ (PBE0/6-311++G(d,p)). Adapted with permission from E. G. Kratz, A. R. Walker, L. Lagardère, F. Lipparini, J.-P. Piquemal, G. Andrés Cisneros J. Comput. Chem. 2016, 37, 1019–1029. Copyright 2016 American Chemical Society

		MAE, $E_{\textnormal{bind}}$
QM	MM	meV	kcal/mol
PBE0	AMOEBA	19.02 (25.54)	0.4385 (0.5889)
PBE0	TIP3P	32.89 (27.75)	0.7584 (0.6398)
-	AMOEBA	9.26 (36.62)	0.2136 (0.8446)
-	TIP3P	44.84 (27.19)	1.0340 (0.6270)

Another recent example of the importance of polarizable environments in QM/MM is shown in the work of Loco et al. [83]. In this work the QM/AMOEBA method, which interfaces the GAUSSIAN [46] and Tinker/Tinker-HP [76] packages has been used to compute excitations properties on the Thiazole Orange dye (TO) intercalated in DNA solvated in water (Figure 2). Several different types of QM/MM molecular dynamics (MD) simulations have been conducted, here we will focus on two of them. First, the QM subsystem is only composed of the TO dye (denoted QM/MM) as previously done by Loco et al. [81]. Second, the QM subsystem is composed of the TO dye and two pairs of nucleobases (denoted QM/MM (PB)). The calculated results indicate charge-transfer intruder states for the larger QM subsystem (QM/MM (PB)), leading to the delocalization of the TO $\rm\pi\pi^{*}$ transition on the QM nucleobases. Though, we observe a small difference of the excitation energy value between QM/MM and QM/MM (PB) when using the AMOEBA polarizable force field. However, complementary studies, using the non-polarizable AMBER99sb force field [62] and the TIP3P water model, have shown a non-negligible difference of the excitation energy between QM/MM and QM/MM (PB). Therefore, computed excitation energies from polarizable QM/MM are less sensitive to the choice of the QM subsystem than classical QM/MM simulations. Embedding polarization effects are a key factor not only for QM/MM calculations but also to compute gas and condensed phase properties.

3 QM/MM Boundary Considerations

The combination of quantum and classical methods results in short– and long–range boundaries that require special consideration. In the short–range regime, the junction between the QM subsystem and the MM environment requires care in cases where covalent bonds are present across the QM/MM boundary. There are scenarios where, even if no chemical bonds are cut, it is necessary to pay special attention to this boundary, for example if molecular exchange between the QM and MM regions is possible or likely. For long–range interactions, it is necessary to consider the role of these effects on the QM/MM energies and forces. In this section we provide a brief review of approaches that have been developed to address both types of scenarios.

3.1 Boundary Between QM and MM subsystems

One of the challenges in QM/MM is the description of atoms at the boundary between QM and MM subsystems involving covalent bonds. In general, atoms that are important in the reactive process are described at the QM level of theory, while the remaining atoms are described at the MM level of theory. In the case where only non-covalent interactions are involved, the separation between the QM and MM subsystems is straightforward. In the cases where a covalent bond (or more) need to be cut to create the two subsystems, the atoms involved in the covalent bond, and atoms near the boundary, need then to be treated with care. Indeed, the two different levels of theories used in QM/MM must be coupled to describe atoms and accurately compute the energy of the system. The challenging description of QM/MM boundary atoms have been addressed in different manners leading to the development of various approaches, some of which are mentioned below.

The link atom is the most commonly used method due to its simplicity and ease of implementation [44, 109]. The idea is to add a supplementary atom (i.e link atom), generally an hydrogen atom, to fill the empty valence in the QM subsystem arising for each cut covalent bond (Figure 3a). The link atom is then only described at the QM level, resulting in the addition of three extra degrees of freedom. This can result in modifications of the electronic structure and the chemical properties of the QM subsystem. Moreover, the closeness between the link atom and the charge of the MM atoms may induce an over-polarization from the MM towards the QM subsystem. In order to overcome those issues, improvements have been suggested, such as the double link atom (DLA) [28] among others. In the DLA method, two supplementary atoms, usually hydrogen atoms, are linked to both atoms of the broken bond to create the QM and MM subsystems. The idea here is to overlap both QM and MM link atoms (QMH and MMH) to ensure an overall neutral net charge, avoiding then an over-polarization from the MM towards the QM (Figure 3a). Also, it is possible to use the double link-atom method with a point charge or a single s–type Gaussian function on the boundary atoms. In the latter, the charges are delocalized, providing a better description of the charge density distribution of the MM boundary atoms.

Methods using frozen localized, hybrid orbitals have been developed to describe QM/MM boundary atoms, following pioneering work from Warshel and Levitt [119]. Later, Rivail and coworkers [90, 42, 111, 8, 91, 43] have developed the Localized Self-Consistent Field (LSCF) method. The idea is to describe the QM atom of the broken bond using a set of hybrid orbitals. The orbital pointing towards the MM atom (shaded gray orbital in Figure 3b) is kept frozen during the QM SCF optimization, while other orbitals (plain white orbitals in Figure 3b) are optimized with all QM orbitals. In a related approach, the Generalized Hybrid Orbital (GHO) method developed by Gao et al. [47, 6], the set of orbitals is placed on the MM atom of the broken bond. In this case, the hybrid orbital pointing towards the QM atom is optimized with all the QM orbitals (plain white orbital in Figure 3c) while the other orbitals on the MM atom are kept frozen during the SCF optimization procedure. Overall, the use of frozen orbitals does not alter the electronic and chemical properties of the system. However, it requires additional parametrization steps of the orbitals depending on the nature of the bond to be cut.

The pseudo-bond approach involves the use of "taylor–made" pseudopotentials. In this method the QM atom of the broken bond is described using effective core potentials (ECP) specifically parametrized to reproduce the environment of the covalent bond. The QM atom is denoted as C_ps and is located in the QM region (Figure 4a). Pseudo-bonds should mimic the original chemical and structural properties of the broken bond between the QM and MM subsystems. Following Zhang et al. previous works [126, 125], Parks et al. [97] have developed new parameters for pseudo-bonds for common bonds found in biological systems such as: C_ps( $sp^{3}$ )-C( $sp^{3}$ ), C_ps( $sp^{3}$ )-C( $sp^{2}$ ,carbonyl), and C_ps( $sp^{3}$ )-N( $sp^{3}$ ). The parametrization of the ECPs are independent of the force fields and the latest pseudobonds also include a minimal basis-set, which overcomes previous dependence of the pseudobond to specific basis sets. In parallel, Dilabio et al. [31] have developed quantum capping potentials, which include additional terms compare to ECP as spherical screening and Pauli repulsion, which contribute then to the refinement of pseudo-bonds.

In the case of polarizable QM/MM, special care needs to be taken to address the polarization interaction across the QM/MM boundary. The pseudobond approach has been extended to encompass polarization effects [74]. In addition to using pseudobonds, boundary atoms (purple shaded surface in Figure 4) need to be considered to compute a QM/MM optimization, which is performed by means of four energy calculations such as : QM energy in the MM static field (Figure 4b), MM energy with no charge on the QM atoms and boundary atoms (Figure 4c), MM polarization energy with charges on the QM atoms (Figure 4d), MM optimization with no polarization from the QM subsystem, pseudo-bond and boundary atoms (Figure 4e). This allow then to avoid the over-polarization from the MM towards the QM, Kratz et al. [74] have demonstrated the accuracy of using this approach on oligopeptide chains in the gas phase with LICHEM [74].

In another recent example, Das et al. [29] studied the polarization effect of uracil DNA glycosylase (UDG) on substrates leading to a ”tautomeric strain” phenomena. UDG belongs to the excision repair enzyme family which its role is to cleave the N₁-C ${}_{1}^{\prime}$ to remove mutagenic uracil from DNA. The activation of the N₁-C ${}_{1}^{\prime}$ bond has been investigated actively both the focus of experimental and computational studies because uracil is not considered as the best of leaving group, giving rise to questions about the reaction mechanism.

Das et al. performed QM calculations on a cluster model based on the active site as well as QM/MM calculations on the full UDG system including eight pseudobonds to investigate the effect of the environment on the aromaticity of uracil, and how this affects the reactivity [29]. Based on the results from both approaches, a step–wise dissociative mechanism for UDG was proposed (Figure 5). Interestingly, the active site of UDG includes several nearby residues interacting with UDG by means of hydrogen bonds. Calculated aromaticity indicators (NICS and HOMED) for uracil in the cluster model suggest that the uracil experiences an increase in aromaticity compared with uracil in the gas phase, to values similar to higher energy uracil tautomers. Including the full polarized protein environment increases the aromaticity of uracil even more. This ”tautomeric strain” enables the reduction of the activation barrier to cleave uracil from UDG thanks to the polarization effects [29].

3.2 Long-Range Electrostatics in QM/MM

Long-range interactions in biomolecular and condensed–phase simulations play an important role that needs to be accounted for in order to properly represent a bulk system. One approach to simulate a bulk system is by the use of periodic boundary conditions (PBC). A popular approach to include long–range effects under PBC in classical simulations is the smooth particle mesh Ewald (sPME) approach [38]. The sPME approach relies on an approximation of the Ewald summation method [39] that uses the fast Fourier transform approach, resulting in $\mathcal{O}(N\log N)$ scaling.

Briefly, for a system with $N$ point charges, the electrostatic energy of a system under PBC is given by:

E_{\textnormal{elst}}=\frac{1}{2}\sum_{n}^{*}\sum_{i,j=1}^{N}\frac{q_{i}q_{j}}{|\textbf{r}_{i}-\textbf{r}_{j}+\textbf{n}|^{2}}

(24)

where the asterisk over the sum denotes that the terms with $\textbf{n}=0$ and either $j=i$ are omitted.

The direct calculation of the electrostatic energy in Eq. (24) is $\mathcal{O}(N^{2})$ and additionally, is only conditionally convergent. The Ewald method [39] divides the electrostatic interaction into a sum of three contributions: direct, reciprocal, and correction terms,

E_{\textnormal{elst}}=E_{\textnormal{direct}}+E_{\textnormal{recip}}+E_{\textnormal{corr}}

(25)

The direct term is given by

E_{\textnormal{dir}}=\frac{1}{2}\sum_{n}^{*}\sum_{i,j=1}^{N}q_{i}q_{j}\left(\frac{\textnormal{erfc}(\beta|\textbf{r}_{j}-\textbf{r}_{j}+\textbf{n}|)}{|\textbf{r}_{j}-\textbf{r}_{i}+\textbf{n}|}\right)

(26)

where $\beta$ is the Ewald splitting parameter. The reciprocal term is defined as

E_{\textnormal{recip}}=\frac{1}{2\pi V}\sum_{m\neq 0}\frac{\textnormal{exp}(-\pi^{2}\textbf{m}^{2}/\beta^{2})}{\textbf{m}^{2}}S(\textbf{m})S(\textbf{-m})

(27)

where $V=\textbf{a}_{1}\cdot\textbf{a}_{2}\times\textbf{a}_{3}$ is the volume of the box, $\textbf{m}=m_{1}\textbf{a}_{1}^{*}+m_{2}\textbf{a}_{2}^{*}+m_{3}\textnormal{a}_{3}^{*}$ is the lattice in reciprocal space and $S(\textbf{m})$ is the structure factor defined by

S(\textbf{m})=\sum_{j=1}^{N}q_{j}\textnormal{exp}(2\pi im\cdot r_{j})

(28)

The correction term in Eq. (25) is composed by the self energy ( $E_{\textnormal{self}})$ .

E_{\textnormal{self}}=-\frac{\beta}{\sqrt{\pi}}\sum_{i=1}^{N}q_{i}^{2}

(29)

The above expressions for the correction can be modified to take into account other corrections such as non–neutral unit cells, polarization response, surface terms, etc [27].

The Ewald method avoids the conditional convergence in Eq. (24), albeit the algorithm is $\mathcal{O}(N^{2})$ . The sPME and related approaches approximate the reciprocal sum by using a discrete convolution on an interpolating grid, resulting in a $\mathcal{O}(NlogN)$ algorithm. This approach is extremely efficient when a discrete representation, i.e. point charges or point multipoles, is employed, because the discrete charge distribution can be interpolated on a grid in a relatively straightforward fashion [27].

The interpolation on a regular grid (or set of grids) of the continuous charge density, such as that from the QM subsystem in QM/MM calculations, presents a challenge. One way to overcome this challenge is to approximate the continuous QM charge density by a discrete representation (e.g. point charges) as was originally proposed by Nam et al. [93] for QM/MM–Ewald, and by Walker et al. [118]. Giese and York developed the ambient potential composite Ewald approach, which avoids the use of dense Fourier grids for the QM electronic density, and performs a direct interaction of the QM charge density with the reciprocal potential [49]. The use of an accurate approach to interact the continuous charge density of the QM subsystem with the charges of the MM environment was shown to avoid SCF instabilities and artifacts in PMF calculations. [49]

Long–range effects can also be approximated by boundary potential models. In the case of QM/MM calculations several methods have been proposed, including the empirical valence bond (EVB/MM) framework [120], the spherical solvent boundary potential (SSBP) [10], the generalized solvent boundary potential (GSBP) [64], the solvated macromolecule boundary potential (SMBP) [11]. One more approach to approximate long–range effects is to use a cutoff method coupled with switching or shifting functions, without taking the long–range effects explicitly into account. [89, 1, 2]

Another method that has been developed for both point–charge, and multipolar/polarizable QM/MM calculations is the QM/MM Long-Range Electrostatic Correction (QM/MM–LREC). The QM/MM–LREC scheme uses a minimum image convention and a switching and smoothing function to avoid edge effects and incorporate long–range electrostatic interactions, while at the same time matching the energies and forces with sPME levels of accuracy [40, 74].

The minimum image convention method considers the interactions between each particle $i$ with all particles $j$ within a cutoff radius $R_{\textnormal{cut}}$ . If the distance $r_{ij}$ between the particles $i$ and $j$ is less than half the length of the primary unit cell then the same image is considered, on the other hand, if the distance between the particles $i$ and $j$ is larger, the neighboring image is used.

In the QM/MM-LREC scheme, each particle $i$ interacts with every particle $j$ within a certain cutoff distance $r_{ij}$ . This leads to a modified Coulomb interaction such that

E(q_{i},q_{j},r_{ij})=\frac{q_{i}q_{j}}{r_{ij}}f(r_{ij}^{\prime},s)

(30)

where

f(r_{ij}^{\prime},s)=\left[1-\left(2r_{ij}^{\prime 3}-3r_{ij}^{\prime 2}+1\right)^{s}\right]

(31)

and

r_{ij}^{\prime}=\left(\frac{R_{\textnormal{cut}}-r_{ij}}{R_{\textnormal{cut}}}\right)

(32)

Here $s$ is an adjustable integer exponent controlling the decay of the damping function (see Figure 6). The LREC scheme avoids the need to approximate the QM charge density for reciprocal space by considering only the first replicas around the primary cell. It also has the advantage that analytical Hessians can be calculated analytically, and the switching function $f(r^{\prime}_{ij},s)$ can be employed for any discrete embedding field including monopoles and higher order multipoles. This approach has been recently extended to a hybrid approach allowing a reduction in the cutoff distance [96].

4 QM/MM with Advanced Potentials

The QM/MM approaches described so far only involve electrostatic and polarization embedding of the QM wavefunction (22). The inclusion of explicit polarization improves the description of the MM environment, however, other effects are approximated by calculating them via the classical terms, e.g. Van der Waals, [105] or are neglected altogether, e.g. charge transfer.

Several approaches have been developed to include a more accurate description of the MM environment. One recently proposed approach is to employ potentials based on the many body decomposition. One such family of potentials is the MB-pol and MB-DFT potentials [104]. We have recently described a new QM/MM implementation involving MB-pol and MB-DFT within LICHEM and within Gaussian to compare partially self–consistent (psc) and fully self consistent (fsc) polarizable QM/MM approaches [77]. These implementations include explicit Coulomb and polarization embedding of the QM wavefunction and a many-body potential representation of the MM environment. The use of MB-pol and MB-DFT allows for a very accurate representation of the environment. This results in various advantages such as no longer needing additional functions or approximations to ensure continuous energy and forces in adaptive QM/MM calculations [35, 121].

Another approach to improve the description of the environment in QM/MM simulations is to use advanced potentials that explicitly include one or more of the missing components in the embedding environment. One approach that allows an explicit term–by–term embedding is by employing potentials that include separate terms for each physical interaction, such as the effective fragment potential (EFP) which can include exchange and charge transfer [117], the fragment exchange potential [21], the QMFF [50], Exchange fragment potential [20], among others. See the recent review of QM:MM and QM:QM embedding approaches by Jones et al. for a thorough review [67].

We have developed a QM/MM implementation within LICHEM that interfaces PSI4 with GEM where the MM subsystem is represented by GEM and thus all MM and QM/MM interactions are calculated using the GEM interaction terms. The Gaussian Electrostatic Model (GEM) is a polarizable potential that uses Gaussian auxiliary basis sets (ABSs) to build the molecular electron density [24]. The philosophy for GEM is to have a density–based force field with errors below 0.2 kcal/mol for total energy and forces compared with high level QM reference. This is achieved by using separate contributions for each of the terms obtained from QM energy decomposition analysis (EDA): Coulomb, Exchange, polarization, charge transfer and dispersion. The GEM densities are obtained by fitting QM relaxed densities to a set of auxiliary Hermite Gaussian basis sets.

Two procedures can be used to perform the adjustment of molecular electron densities through Hermite Gaussians: the conventional analytical variational Coulomb adjustment or the numerical adjustment of the molecular electrostatic potential. The analytical density fitting method relies on the use of Gaussian auxiliary bas functions to expand molecular electron density

\tilde{\rho}(r)=\sum_{k}c_{k}K(r)

(33)

where $\tilde{\rho}$ is the approximate density and $K(r)$ are Hermite Gaussians. The expansion coefficients $c_{k}$ can be obtained by minimizing the self-energy of the error in the density according to some metric

E_{\textnormal{self}}=\langle\rho(r)-\tilde{\rho}(r)|\hat{O}|\rho(r)-\tilde{\rho}(r)\rangle

(34)

Different operators $\hat{O}$ can be used, including the coulomb operator $\hat{O}=1/r$ or the damped Coulomb operator $\hat{O}=\textnormal{erfc}(\beta r)/r$ [71]. The most common is to use the overlap operator $\hat{O}=1$ . The minimization of eq. (34) with respect to the expansion coefficients $c_{k}$ origin a linear system of equations:

\frac{\partial E_{\textnormal{self}}}{\partial c_{l}}=-\sum_{\mu,\nu}P_{\mu\nu}\langle\mu\nu|\hat{O}|l\rangle+\sum_{k}c_{k}\langle k|\hat{O}|l\rangle

(35)

where $P_{\mu\nu}$ is the density matrix. The solution of Eq. (35) requires the inversion of a the ABS matrix $\textbf{G}=\langle k|\hat{O}|l\rangle$ . Although in principle G must be positive definite and symmetric, in practice, it is almost singular and, therefore, the diagonalization to obtain its inverse must be done carefully. To this end, both analytical and numerical approaches have been developed to perform the fitting of the molecular densities [25, 22].

The latest GEM water model has been fitted to match the total energy of the dimer surface and hexamer clusters reported by Paesani et al., [9] as well as the individual GEM terms (Coulomb, exchange, polarization and Van der Waals) to reference symmetry adapted perturbation theory (SAPT) with aug–cc–pVTZ [34]. The fitted density for GEM was calculated at the CCSD(T)/aug–cc–pVTZ level. [34, 33] The auxiliary basis sets employed involve 70 primitives, resulting in approximately 280,000 functions for a system of 512 waters. T he integrals required for the calculation of the GEM terms are evaluated using an extension of the sPME approach for Gaussian functions, [34, 33] which has been implemented in the gem.pmemd program, which is available with the AmberTools release [18].

The first implementation of the QM/GEM approach involved only the Coulomb embedding term [26]. The use of GEM provides an accurate polarization environment for the QM wavefuntion and avoids penetration error effects because of the description of the MM environment with explicit (frozen) densities compared with conventional point charge force fields as shown in Figure 7.

For the current implementation, the total QM/MM interaction energy between the subsystems, $E^{\textnormal{QM/GEM}}$ , involves four terms: Coulomb, Exchange–repulsion, polarization, and dispersion: $E_{\textnormal{Coul}}^{\textnormal{QM/GEM}}+E_{\textnormal{exch}}^{\textnormal{QM/GEM}}+E_{\textnormal{pol}}^{\textnormal{QM/GEM}}+E_{\textnormal{disp}}^{\textnormal{QM/GEM}}$ [54].

The Coulomb and exchange–repulsion terms involve 3 center integrals between the QM density and the fitted GEM density:

E_{\textnormal{Coul}}^{\textnormal{QM/GEM}}=\int\int\frac{\rho(r_{1})\tilde{\rho}(r_{2})}{r_{12}}dr_{1}dr_{2}

(36)

and

E_{\textnormal{exch}}^{\textnormal{QM/GEM}}=K_{\textnormal{exch}}\int\int\rho(r_{1})\tilde{\rho}(r_{2})dr_{1}dr_{2},

(37)

where $K_{\textnormal{exch}}$ is a proportionality constant [122, 99, 25]. The polarization interaction is calculated using the same approximation described above for LICHEM with AMOEBA. The dispersion term is approximated by a multipolar expansion taking the 6, 8 and 12 terms into consideration. The exchange proportionality parameter, as well as the coefficients for the dispersion term have been parametrized by linear least squares to match the SAPT2+3/aug–cc–pVTZ components for the ten water dimers [112].

In this implementation we considered two approaches for the calculation of the total interaction depending on the calculation of the QM/GEM exchange–repulsion term. In both cases the QM/GEM dispersion terms are added a posteriori to the total energy (and forces). The Coulomb term is calculated by including the GEM density in the core Hamiltonian in both alternative approaches. Thus, the difference in the two approaches involves only the exchange–repulsion component. In one approach, the exchange–repulsion inter–molecular interaction is calculated after the SCF has completed The second approach involves the inclusion of the exchange–repulsion term explicitly in the SCF, i.e. explicit exchange-repulsion embedding. This is achieved by including the GEM density in the core Hamiltonian:

H_{\textnormal{eff}}=H_{\textnormal{core}}+V_{\textnormal{GEM}}

(38)

in which $V_{GEM}$ involves the introduction of the 3–center integrals of the QM MOs with the GEM density in the core Hamiltonian by:

V_{\textnormal{GEM}}=\sum_{l}x_{l}\sum_{\mu\nu}\langle\mu\nu\parallel l\rangle+K^{\prime}_{exch}\sum_{l}x_{l}\sum_{\mu\nu}\langle\mu\nu\mid l\rangle,

(39)

where the first term corresponds to the Coulomb integrals and the second term to the overlap integrals multiplied by the exchange–repulsion proportionality constant. In this second case the exchange–repulsion proportionality constant is different than the former case since the QM wavefunction experiences a different external potential (see below). In all cases the required integrals have been programmed into a modified version of Psi4 [98].

The QM/GEM implementation in LICHEM was assesed by comparing the total intermolecular interaction Energies, as well as the individual interaction terms for a sub–set of the water dimer potential reported by Babin et al. [9] as described in Ref. [54]. The maximum unsigned error (MUE), and root mean squared error (RMSE) for the clustered dimers is shown in Figure 8. In the original comparison, the term–by–term comparison was done with respect to SAPT2+3/aug–cc–pVTZ, which does not provide a straightforward separation of the induction and polarization components. Therefore, the dispersion+polarization terms were analyzed together. As can be seen from Figure 8, the use of GEM for the Coulomb, polarization and exchange–repulsion embedding provides an accurate environment for the hybrid calculation.

5 LICHEM

Various methods described above including the pseudobond, and QM/MM–LREC, have been implemented in the LICHEM software (Layered Interacting CHEmical Models) package in the context of both non–polarizable and polarizable potentials including AMOEBA, GEM, and MB–Pol. [75, 54]. At its core, LICHEM is an interface that can couple QM codes such as Gaussian [46], NWChem [113] and Psi4 [98] and MM codes such as TINKER [102], TINKER-HP [76] and LAMMPS [100].

When AMOEBA (or another multipolar force field) is employed, LICHEM uses an approximation for the classical multipolar environment to allow the seamless integration with multiple QM codes [30, 74] LICHEM provides a variety of approaches for non–polarizable and polarizable QM/MM simulations including single point energies, Monte Carlo (Grand Canonical and Isothermal–Isobaric ensembles), and path integral Monte Carlo. Geometry optimizations can be performed in an iterative procedure, where QM and MM atoms are optimized separately. QM atoms are optimized including MM charges or multipoles. Then, MM atoms are optimized including point charges comes from electron density from QM atoms. This procedure continue until convergence criteria are met.

Geometry optimization algorithms for single structures are limited to critical points along a surface. In many cases, it is useful to determine the entire reaction path that connects two, or several critical points to investigate a reaction mechanism. An entire reaction path can be optimized by performing an optimization with multiple replicas between the reactant and product states. Several algorithms have been developed for this end. One such family is called the chain–of–replica methods. These methods represent a path as a string of beads connected to each other. LICHEM has two chain–of–replica algorithms for path optimization, the Nudged Elastic Band (NEB) [68, 60], and the Quadratic String Method (QSM) [17].

The nudged elastic band (NEB) method [68] is an efficient method for finding the minimum energy path (MEP) between a given initial and final state of a transition. To sample regions between stationary points (minima or transition states), spring forces can be added to ensure continuity of the path. In the NEB method, the modified forces ( $F$ ’) of the elastic band on the bead $p$ and the step $n$ are given by:

F_{p,n}^{\prime}=F_{p,n}+K\Delta R_{p,n}

(40)

where $K$ is the spring constant and $\Delta R$ is the displacement between replicas defined as:

\Delta R_{p,n}=|R_{p+1,n}-R_{p,n}|-|R_{p,n}-R_{p-1,n}|

(41)

Spring forces contribute to the elastic band method in locating the minimum energy path. However, to prevent the beads from slipping to the minimum, it is necessary to eliminate the forces parallel to the reaction. The forces are then given by

F_{p,n}^{\prime}=F_{p,n}-F_{p,n}\cdot\tau_{p,n}\tau_{p,n}+K\Delta R_{p,n}\tau{p,n}

(42)

where $\tau$ is a tangent unit vector indicating the direction of the reaction path.

The climbing image (CI) NEB method is a modification to the NEB method. The IC travels upward along the tangent to find an approximate transition state structure. The transition state forces in the IC (replica c) are expressed as:

F_{c,n}^{\prime}=F_{c,n}-2F_{c,n}\cdot\tau_{c,n}\tau_{c,n}

where the forces are inverted along the tangent [60].

The Quadratic String Method (QSM) [17] does not require a predetermined step size and eliminates the need for an artificial spring constant. Instead, the approximate Hessian confidence radius and surface are updated at the end of each iteration. First, the QSM algorithm calculates the energy and gradients of all user-defined beads. The approximate Hessians ( $H_{k}$ ) constructed from the gradients in the initial step are updated with the damped Broyden-Fletcher-Goldfarb-Shanno (DBFGS) algorithm [45]

H^{i+1}=H^{i}-\frac{H^{i}\delta^{i}(\delta^{i})^{T}H^{i}}{(\delta^{i})^{T}H^{i}\delta^{i}}+\frac{r^{i}(r^{i})^{T}}{(r^{i})^{T}\delta^{i}}

(43)

where $i$ indicates to the iteration number and

\delta^{i}=x^{i+1}-x^{i}

(44)

r^{i}=\theta^{i}\gamma^{i}+(1-\theta^{i})H^{i}\delta^{i}

(45)

with

y^{i}=g^{i+1}-g^{i}

(46)

\theta^{i}=\begin{cases}1,\quad\text{if}\quad(\delta^{i})^{T}\gamma^{i}>0.2(\delta^{i})^{T}H^{i}\delta^{i}\\ \displaystyle\frac{0.8(\delta^{i})^{T}H^{i}\delta^{i}}{(\delta^{i})^{T}H^{k}\delta^{k}-(\delta^{i})^{T}},\quad\text{otherwise}\end{cases}

(47)

The update of $H_{k}$ is followed by the update of the trust radii using the energy as a merit function;

\rho=\frac{E_{k}^{i+1}-E_{k}^{i}}{dx_{k}^{T}g_{k}^{i}+\displaystyle\frac{1}{2}dx_{k}^{T}H_{k}^{i}dx_{k}}

(48)

Both NEB and QSM require optimized structures for the end–points (e.g. reactant and product) as well as an initial guess of the intermediate images (beads). One approach to generate this initial guess is to use a linear interpolation between the end–point coordinates. In QM/MM simulations this is usually accomplished using only the coordinates of the QM subsystem, and use the MM environment from one of the ends for all the beads. However, this can leads to issues during the path optimization such as discontinuities along the trajectory. One approach to avoid this is to employ a restraints MM optimization approach [124, 23].

In this case, the initial stages of the iterative QM/MM path optimization are performed with restraints on the MM atoms during the MM part of the optimization. These restraints are reduced at each subsequent iteration. Figure 9 shows a flowchart that describes these approaches in LICHEM.

One more algorithmic advantage that can be applied in the case of "chain–of–replica" approaches is the fact that each bead can be calculated independently. Therefore, it is possible to employ hybrid-parallelization paradigms such that the QM or MM parts of the calculation employ the parallelism inherent in the respective codes (i.e. MPI, OpenMP, cuda, etc) and the path optimizer used MPI to run all beads in parallel. This approach has been implemented in LICHEM for QSM [53]. Figure 10 shows performance metrics for a test reaction comparing NEB, QSM and QSM with restrained MM optimization.

One example of the use of the methods mentioned above for a QM/MM calculation with AMOEBA was reported by Vázquez-Montelongo et al. [116], who performed polarizable QM/MM simulations to investigate the N-tert-Butyloxycarbonylation (Boc) of aniline in ionic liquid mixture composed of water/ [EMIm][BF₄]. Boc is a common group used to protect molecules in organic synthesis [106].

Vázquez-Montelongo et al. have assessed the role of the solvent in four different possible reaction mechanisms of N-tert-Boc protection of aniline using LICHEM. Two of them require an activation of the Boc group by a water molecule or a cation ([EMIm]⁺), while the other two do not (denoted configuration C2). Here we will focus on the two possible mechanisms of configuration C2 which are shown in Figures 11a and 11d. Reaction (1) (Fig. 11a, denoted scheme 1c in the original study) represents a sequential mechanism: first, the nucleophic attack of the aniline to the carbonyl of one of the Boc groups occurs, followed by the formation of CO₂. Reaction (2) corresponds to the concerted mechanism (Fig. 11d, denoted scheme 1d in the original study) where the nucleophic attack of the aniline to the carbonyl of one of the Boc groups and the formation of CO₂ occur concurrently.

The path optimization methods described above have been applied to obtain minimum energy paths (MEP) for both reaction schemes (Figures 11b and 11e). The calculated rate–limiting barrier for reaction (1) is 21.38 kcal/mol compared with 20.42 kcal/mol for reaction (2). In order to assess the effect of polarization on reaction paths, the MEPs were re-computed with the AMOEBA polarization term set to zero during the calculation (Figures 11c and 11f). In this case, MEP with or without the polarization terms give similar barriers for reaction (1), although the shape of the reaction path without polarization is drastically altered, and the intermediate and TS2 observed in the original path are no longer observed. In the case of reaction (2) the differences between the calculated paths with and without polarization are more significant. The energy barrier computed between the reactant and product increases significantly when the polarization is not taken into account, and the shape of the path is also affected. This example illustrates the necessity to include explicit polarization within QM/MM calculations in highly charged environments.

6 Conclusions

The computational study of reactions in biological and condensed-phase systems has been advanced the combination of QM and MM methods. Here, we have reported recent advances in the development of QM/MM including polarizable force fields. The use of polarizable potentials such as AMOEBA, gives rise to novel issues that need to be considered, including challenges at the internal and external boundaries (long–range effects) and coupling to the QM software. This review focused on the methods and approaches that have been developed and implemented in LICHEM to address these, and other challenges, such as the extension of the pseudobond approach, the development and implementation of QM/MM–LREc, the implementation of restrained MM path optimizations coupled with the quadratic string method, implementation of advanced force fields for QM/MM such as MB–Pol and GEM, etc. Several examples have been provided to underscore the utility of polarizable potentials and their usefulness in the study of complex biological and condensed phase systems.

References

[1] Orlando Acevedo. Simulating Chemical Reactions in Ionic Liquids Using QM/MM Methodology. The Journal of Physical Chemistry A, 118(50):11653–11666, dec 2014.
[2] Orlando Acevedo and Wiliiam L Jorgensen. Quantum and molecular mechanical Monte Carlo techniques for modeling condensed-phase reactions. WIREs Computational Molecular Science, 4(5):422–435, sep 2014.
[3] Shideh Ahmadi, Lizandra Barrios Herrera, Morteza Chehelamirani, Jiří Hostaš, Said Jalife, and Dennis R Salahub. Multiscale modeling of enzymes: QM-cluster, QM/MM, and QM/MM/MD: A tutorial review. International Journal of Quantum Chemistry, 118(9):e25558, 2018.
[4] Norman L Allinger, Kuohsiang Chen, and Jenn-Huei Lii. An improved force field (MM4) for saturated hydrocarbons. Journal of Computational Chemistry, 17(5-6):642–668, apr 1996.
[5] Norman L Allinger, Young H Yuh, and Jenn Huei Lii. Molecular mechanics. The MM3 force field for hydrocarbons. 1. Journal of the American Chemical Society, 111(23):8551–8566, nov 1989.
[6] Patricia Amara, Martin J. Field, Cristobal Alhambra, and Jiali Gao. The generalized hybrid orbital method for combined quantum mechanical/molecular mechanical calculations: formulation and tests of the analytical derivatives. Theoretical Chemistry Accounts: Theory, Computation, and Modeling (Theoretica Chimica Acta), 104(5):336–343, aug 2000.
[7] Johan Åqvist and Arieh Warshel. Simulation of enzyme reactions using valence bond force fields and other hybrid quantum/classical approaches. Chemical Reviews, 93(7):2523–2544, nov 1993.
[8] Xavier Assfeld and Jean-Louis Rivail. Quantum chemical computations on parts of large molecules: the ab initio local self consistent field method. Chemical Physics Letters, 263(1-2):100–106, dec 1996.
[9] Volodymyr Babin, Claude Leforestier, and Francesco Paesani. Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient. Journal of Chemical Theory and Computation, 9(12):5395–5403, 2013.
[10] Dmitrii Beglov and Benoît Roux. Finite representation of an infinite bulk system: Solvent boundary potential for computer simulations. The Journal of Chemical Physics, 100(12):9050–9063, 1994.
[11] Tobias Benighaus and Walter Thiel. Long-Range Electrostatic Effects in QM/MM Studies of Enzymatic Reactions: Application of the Solvated Macromolecule Boundary Potential. Journal of Chemical Theory and Computation, 7(1):238–249, jan 2011.
[12] Mattia Bondanza, Michele Nottoli, Lorenzo Cupellini, Filippo Lipparini, and Benedetta Mennucci. Polarizable embedding QM/MM: the future gold standard for complex (bio)systems? Phys. Chem. Chem. Phys., 22(26):14433–14448, 2020.
[13] Eliot Boulanger, Lei Huang, Chetan Rupakheti, Alexander D MacKerell Jr, and Benoît Roux. Optimized Lennard-Jones Parameters for Druglike Small Molecules. Journal of chemical theory and computation, 14(6):3121–3131, jun 2018.
[14] Eliot Boulanger and Walter Thiel. Solvent Boundary Potentials for Hybrid QM/MM Computations Using Classical Drude Oscillators: A Fully Polarizable Model. Journal of Chemical Theory and Computation, 8(11):4527–4538, nov 2012.
[15] Bernard R Brooks, Robert E Bruccoleri, Barry D Olafson, David J States, S Swaminathan, and Martin Karplus. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. Journal of Computational Chemistry, 4(2):187–217, 1983.
[16] Elizabeth Brunk and Ursula Rothlisberger. Mixed Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States. Chemical Reviews, 115(12):6217–6263, jun 2015.
[17] Steven K Burger and Weitao Yang. Quadratic string method for determining the minimum-energy path based on multiobjective optimization. The Journal of Chemical Physics, 124(5):54109, 2006.
[18] D.A. Case, K. Belfon, I.Y. Ben-Shalom, S.R. Brozell, D.S. Cerutti, T.E. Cheatham, V.W.D. Cruzeiro, T.A. Darden, R.E. Duke, G. Giambasu, M.K. Gilson, H. Gohlke, A.W. Goetz, R. Harris, S. Izadi, S.A. Izmailov, K. Kasavajhala, A. Kovalenko, R. Krasny, T. Kurtzman, T.S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, V. Man, K.M. Merz, Y. Miao, O. Mikhailovskii, G. Monard, H. Nguyen, A. Onufriev, F. Pan, S. Pantano, R. Qi, D.R. Roe, A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, C.L. Simmerling, N.R.Skrynnikov, J. Smith, J. Swails, R.C. Walker, J. Wang, L. Wilson, R.M. Wolf, X. Wu, Y. Xiong, Y. Xue, D.M. York, and P.A. Kollman. AMBER 2020. University of California, San Francisco., 2020.
[19] David A Case, Thomas E Cheatham III, Tom Darden, Holger Gohlke, Ray Luo, Kenneth M Merz Jr., Alexey Onufriev, Carlos Simmerling, Bing Wang, and Robert J Woods. The Amber biomolecular simulation programs. Journal of Computational Chemistry, 26(16):1668–1688, 2005.
[20] Xin Chen and Jiali Gao. Fragment Exchange Potential for Realizing Pauli Deformation of Interfragment Interactions. The Journal of Physical Chemistry Letters, 11(10):4008–4016, may 2020.
[21] Xin Chen, Zexing Qu, Bingbing Suo, and Jiali Gao. A self-consistent coulomb bath model using density fitting. Journal of Computational Chemistry, 41(18):1698–1708, 2020.
[22] G Andrés Cisneros, Dennis Elking, Jean-Philip Piquemal, and Thomas A Darden. Numerical fitting of molecular properties to Hermite Gaussians. The journal of physical chemistry. A, 111(47):12049–12056, nov 2007.
[23] G Andrés Cisneros, Haiyan Liu, Zhenyu Lu, and Weitao Yang. Reaction path determination for quantum mechanical/molecular mechanical modeling of enzyme reactions by combining first order and second order "chain-of-replicas" methods. The Journal of chemical physics, 122(11):114502, mar 2005.
[24] G Andrés Cisneros, Jean-Philip Piquemal, and Thomas A Darden. Intermolecular electrostatic energies using density fitting. The Journal of Chemical Physics, 123(4):44109, 2005.
[25] G Andrés Cisneros, Jean-Philip Piquemal, and Thomas A Darden. Generalization of the Gaussian electrostatic model: Extension to arbitrary angular momentum, distributed multipoles, and speedup with reciprocal space methods. The Journal of Chemical Physics, 125(18):184101, 2006.
[26] G Andrés Cisneros, Jean-Philip Piquemal, and Thomas A Darden. Quantum Mechanics/Molecular Mechanics Electrostatic Embedding with Continuous and Discrete Functions. The Journal of Physical Chemistry B, 110(28):13682–13684, jul 2006.
[27] Thomas A Darden. Extensions of the Ewald method for Coulomb interactions in crystals. International Tables for Crystallography, pages 458–481, jun 2010.
[28] Debananda Das, Kirsten P. Eurenius, Eric M. Billings, Paul Sherwood, David C. Chatfield, Milan Hodošček, and Bernard R. Brooks. Optimization of quantum mechanical molecular mechanical partitioning schemes: Gaussian delocalization of molecular mechanical charges and the double link atom method. The Journal of Chemical Physics, 117(23):10534–10547, dec 2002.
[29] Ranjita Das, Erik A Vázquezvázquez-Montelongo, G Andrés, Andrés Cisneros, and Judy I Wu. Ground State Destabilization in Uracil DNA Glycosylase: Let’s Not Forget "Tautomeric Strain" in Substrates. Journal of the American Chemical Society, 2019.
[30] Mike Devereux, Shampa Raghunathan, Dmitri G Fedorov, and Markus Meuwly. A Novel, Computationally Efficient Multipolar Model Employing Distributed Charges for Molecular Dynamics Simulations. Journal of Chemical Theory and Computation, 10(10):4229–4241, oct 2014.
[31] Gino A. DiLabio, Margaret M. Hurley, and Phillip A. Christiansen. Simple one-electron quantum capping potentials for use in hybrid QM/MM studies of biological molecules. The Journal of Chemical Physics, 116(22):9578–9584, jun 2002.
[32] Fernanda Duarte, Beat A Amrein, David Blaha-Nelson, and Shina C L Kamerlin. Recent advances in QM/MM free energy calculations using reference potentials. Biochimica et biophysica acta, 1850(5):954–965, may 2015.
[33] Robert E Duke and G Andrés Cisneros. Ewald-based methods for Gaussian integral evaluation: application to a new parameterization of GEM. Journal of molecular modeling, 25(10):307, sep 2019.
[34] Robert E Duke, Oleg N Starovoytov, Jean-Philip Piquemal, and G Andrés Cisneros. GEM*: A Molecular Electronic Density-Based Force Field for Molecular Dynamics Simulations. Journal of Chemical Theory and Computation, 10(4):1361–1365, apr 2014.
[35] Adam Duster, Chun-Hung Wang, and Hai Lin. Adaptive QM/MM for Molecular Dynamics Simulations: 5. On the Energy-Conserved Permuted Adaptive-Partitioning Schemes. Molecules, 23(9):2170, aug 2018.
[36] Jacek Dziedzic, Teresa Head-Gordon, Martin Head-Gordon, and Chris-Kriton Skylaris. Mutually polarizable QM/MM model with in situ optimized localized basis functions. The Journal of Chemical Physics, 150(7):074103, feb 2019.
[37] Jacek Dziedzic, Yuezhi Mao, Yihan Shao, Jay Ponder, Teresa Head-Gordon, Martin Head-Gordon, and Chris-Kriton Skylaris. TINKTEP: A fully self-consistent, mutually polarizable QM/MM approach based on the AMOEBA force field. The Journal of Chemical Physics, 145(12):124106, sep 2016.
[38] Ulrich Essmann, Lalith Perera, Max L Berkowitz, Tom Darden, Hsing Lee, and Lee G Pedersen. A smooth particle mesh Ewald method. The Journal of Chemical Physics, 103(19):8577–8593, nov 1995.
[39] P P Ewald. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annalen der Physik, 369(3):253–287, jan 1921.
[40] Dong Fang, Robert E Duke, and G Andrés Cisneros. A new smoothing function to introduce long-range electrostatic effects in QM/MM calculations. The Journal of Chemical Physics, 143(4):44103, 2015.
[41] Dong Fang, Richard L Lord, and G Andrés Cisneros. Ab Initio QM/MM Calculations Show an Intersystem Crossing in the Hydrogen Abstraction Step in Dealkylation Catalyzed by AlkB. The Journal of Physical Chemistry B, 117(21):6410–6420, may 2013.
[42] György G. Ferenczy, Jean-Louis Rivail, Péter R. Surján, and Gábor Náray-Szabó. NDDO fragment self-consistent field approximation for large electronic systems. Journal of Computational Chemistry, 13(7):830–837, sep 1992.
[43] Nicolas Ferré, Xavier Assfeld, and Jean-Louis Rivail. Specific force field parameters determination for the hybrid ab initio QM/MM LSCF method. Journal of Computational Chemistry, 23(6):610–624, apr 2002.
[44] Martin J Field, Paul A Bash, and Martin Karplus. A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations. Journal of Computational Chemistry, 11(6):700–733, 1990.
[45] R Fletcher. Practical Methods of Optimization; (2nd Ed.). Wiley-Interscience, USA, 1987.
[46] M J Frisch, G W Trucks, H B Schlegel, G E Scuseria, M A Robb, J R Cheeseman, G Scalmani, V Barone, G A Petersson, H Nakatsuji, X Li, M Caricato, A V Marenich, J Bloino, B G Janesko, R Gomperts, B Mennucci, H P Hratchian, J V Ortiz, A F Izmaylov, J L Sonnenberg, D Williams-Young, F Ding, F Lipparini, F Egidi, J Goings, B Peng, A Petrone, T Henderson, D Ranasinghe, V G Zakrzewski, J Gao, N Rega, G Zheng, W Liang, M Hada, M Ehara, K Toyota, R Fukuda, J Hasegawa, M Ishida, T Nakajima, Y Honda, O Kitao, H Nakai, T Vreven, K Throssell, J A Montgomery Jr., J E Peralta, F Ogliaro, M J Bearpark, J J Heyd, E N Brothers, K N Kudin, V N Staroverov, T A Keith, R Kobayashi, J Normand, K Raghavachari, A P Rendell, J C Burant, S S Iyengar, J Tomasi, M Cossi, J M Millam, M Klene, C Adamo, R Cammi, J W Ochterski, R L Martin, K Morokuma, O Farkas, J B Foresman, and D J Fox. Gaussian 16 Revision C.01. Gaussian Inc. Wallingford CT, 2016.
[47] Jiali Gao, Patricia Amara, Cristobal Alhambra, and Martin J. Field. A Generalized Hybrid Orbital (GHO) Method for the Treatment of Boundary Atoms in Combined QM/MM Calculations. The Journal of Physical Chemistry A, 1998.
[48] Daan P Geerke and Wilfred F van Gunsteren. On the Calculation of Atomic Forces in Classical Simulation Using the Charge-on-Spring Method To Explicitly Treat Electronic Polarization. Journal of Chemical Theory and Computation, 3(6):2128–2137, nov 2007.
[49] Timothy J Giese and Darrin M York. Ambient-Potential Composite Ewald Method for ab Initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulation. Journal of Chemical Theory and Computation, 12(6):2611–2632, jun 2016.
[50] Timothy J Giese and Darrin M York. Quantum mechanical force fields for condensed phase molecular simulations. Journal of physics. Condensed matter : an Institute of Physics journal, 29(38):383002, sep 2017.
[51] Tommaso Giovannini, Alessandra Puglisi, Matteo Ambrosetti, and Chiara Cappelli. Polarizable QM/MM Approach with Fluctuating Charges and Fluctuating Dipoles: The QM/FQF $\mu$ Model. Journal of Chemical Theory and Computation, 15(4):2233–2245, apr 2019.
[52] Tommaso Giovannini, Rosario Roberto Riso, Matteo Ambrosetti, Alessandra Puglisi, and Chiara Cappelli. Electronic transitions for a fully polarizable QM/MM approach based on fluctuating charges and fluctuating dipoles: Linear and corrected linear response regimes. The Journal of Chemical Physics, 151(17):174104, nov 2019.
[53] Hatice Gökcan, Eric Kratz, Thomas A Darden, Jean-Philip Piquemal, and G Andrés Cisneros. QM/MM Simulations with the Gaussian Electrostatic Model: A Density-based Polarizable Potential. The Journal of Physical Chemistry Letters, 9(11):3062–3067, jun 2018.
[54] Hatice Gökcan, Erik Antonio Vázquez-Montelongo, and G. Andrés Cisneros. LICHEM 1.1: Recent Improvements and New Capabilities. Journal of Chemical Theory and Computation, 15(5):3056–3065, may 2019.
[55] Mark S Gordon, Dmitri G Fedorov, Spencer R Pruitt, and Lyudmila V Slipchenko. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chemical Reviews, 112(1):632–672, jan 2012.
[56] Krishna Govender, Jiali Gao, and Kevin J Naidoo. AM1/d-CB1: A Semiempirical Model for QM/MM Simulations of Chemical Glycobiology Systems. Journal of chemical theory and computation, 10:4694–4707, 2014.
[57] Nohad Gresh, G Andrés Cisneros, Thomas A Darden, and Jean-Philip Piquemal. Anisotropic, Polarizable Molecular Mechanics Studies of Inter- and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy. Journal of chemical theory and computation, 3(6):1960–1986, nov 2007.
[58] Muhammad A Hagras and William J Glover. Polarizable Embedding for Excited-State Reactions: Dynamically Weighted Polarizable QM/MM. Journal of Chemical Theory and Computation, 14(4):2137–2144, apr 2018.
[59] Thomas A Halgren. MMFF VI. MMFF94s option for energy minimization studies. Journal of Computational Chemistry, 20(7):720–729, may 1999.
[60] Graeme Henkelman, Blas P Uberuaga, and Hannes Jónsson. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. The Journal of Chemical Physics, 113(22):9901–9904, 2000.
[61] Asbjørn Holt and Gunnar Karlström. Improvement of the NEMO potential by inclusion of intramolecular polarization. International Journal of Quantum Chemistry, 109(6):1255–1266, may 2009.
[62] Viktor Hornak, Robert Abel, Asim Okur, Bentley Strockbine, Adrian Roitberg, and Carlos Simmerling. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Structure, Function, and Bioinformatics, 65(3):712–725, nov 2006.
[63] Hao Hu and Weitao Yang. Development and application of ab initio QM/MM methods for mechanistic simulation of reactions in solution and in enzymes. Theochem, 898(1-3):17–30, mar 2009.
[64] Wonpil Im, Simon Bernèche, and Benoıt Roux. Generalized solvent boundary potential for computer simulations. The Journal of Chemical Physics, 114(7):2924–2937, 2001.
[65] Pathumwadee Intharathep, Anan Tongraar, and Kritsana Sagarik. Ab initio QM/MM dynamics of H3O+ in water. Journal of Computational Chemistry, 27(14):1723–1732, 2006.
[66] Zhenming Jin, Xiaoyu Du, Yechun Xu, Yongqiang Deng, Meiqin Liu, Yao Zhao, Bing Zhang, Xiaofeng Li, Leike Zhang, Chao Peng, Yinkai Duan, Jing Yu, Lin Wang, Kailin Yang, Fengjiang Liu, Rendi Jiang, Xinglou Yang, Tian You, Xiaoce Liu, Xiuna Yang, Fang Bai, Hong Liu, Xiang Liu, Luke W Guddat, Wenqing Xu, Gengfu Xiao, Chengfeng Qin, Zhengli Shi, Hualiang Jiang, Zihe Rao, and Haitao Yang. Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors. Nature, 582(7811):289–293, 2020.
[67] Leighton O Jones, Martín A Mosquera, George C Schatz, and Mark A Ratner. Embedding Methods for Quantum Chemistry: Applications from Materials to Life Sciences. Journal of the American Chemical Society, 142(7):3281–3295, feb 2020.
[68] H. Jónsson, G. Mills, and K. W. Jacobsen. Nudged elastic band method for finding minimum energy paths of transitions. In Classical and Quantum Dynamics in Condensed Phase Simulations, pages 385–404. WORLD SCIENTIFIC, jun 1998.
[69] William L. Jorgensen, Jayaraman Chandrasekhar, Jeffry D. Madura, Roger W. Impey, and Michael L. Klein. Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics, 79(2):926–935, jul 1983.
[70] William L Jorgensen, David S Maxwell, and Julian Tirado-Rives. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. Journal of the American Chemical Society, 118(45):11225–11236, nov 1996.
[71] Yousung Jung, Alex Sodt, Peter M W Gill, and Martin Head-Gordon. Auxiliary basis expansions for large-scale electronic structure calculations. Proceedings of the National Academy of Sciences, 102(19):6692–6697, 2005.
[72] Shina C L Kamerlin and Arieh Warshel. The EVB as a quantitative tool for formulating simulations and analyzing biological and chemical reactions. Faraday Discuss., 145(0):71–106, 2010.
[73] Teruhisa S. Komatsu, Yohei; Koyama, Noriaki; Okimoto, Gentaro; Morimoto, Yousuke; Ohno, and Makoto Taiji. COVID-19 related trajectory data of 10 microseconds all atom molecular dynamics simulation of SARS-CoV-2 dimeric main protease. Mendeley Data, 2020.
[74] Eric G Kratz, Robert E Duke, and G Andrés Cisneros. Long-range electrostatic corrections in multipolar/polarizable QM/MM simulations. Theoretical Chemistry Accounts, 135(7):166, 2016.
[75] Eric G. Kratz, Alice R. Walker, Louis Lagardère, Filippo Lipparini, Jean-Philip Piquemal, and G. Andrés Cisneros. LICHEM: A QM/MM program for simulations with multipolar and polarizable force fields. Journal of Computational Chemistry, 37(11):1019–1029, apr 2016.
[76] Louis Lagardère, Luc-Henri Jolly, Filippo Lipparini, Félix Aviat, Benjamin Stamm, Zhifeng F Jing, Matthew Harger, Hedieh Torabifard, G Andrés Cisneros, Michael J Schnieders, Nohad Gresh, Yvon Maday, Pengyu Y Ren, Jay W Ponder, and Jean-Philip Piquemal. Tinker-HP: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields. Chem. Sci., 9(4):956–972, 2018.
[77] Eleftherios Lambros, Filippo Lipparini, G Andres Cisneros, and Francesco Paesani. A Many-Body, Fully Polarizable Approach to QM/MM Simulations. ChemRxiv, sep 2020.
[78] Pengfei Li, Lin Frank Song, and Kenneth M Merz. Parameterization of Highly Charged Metal Ions Using the 12-6-4 LJ-Type Nonbonded Model in Explicit Water. The Journal of Physical Chemistry B, 119(3):883–895, jan 2015.
[79] Hai Lin and Donald G Truhlar. QM/MM: what have we learned, where are we, and where do we go from here? Theoretical Chemistry Accounts, 117(2):185, 2006.
[80] Filippo Lipparini and Vincenzo Barone. Polarizable Force Fields and Polarizable Continuum Model: A Fluctuating Charges/PCM Approach. 1. Theory and Implementation. Journal of Chemical Theory and Computation, 7(11):3711–3724, nov 2011.
[81] Daniele Loco, Sandro Jurinovich, Lorenzo Cupellini, Maximilian F. S. J. Menger, and Benedetta Mennucci. The modeling of the absorption lineshape for embedded molecules through a polarizable QM/MM approach. Photochemical & Photobiological Sciences, 17(5):552–560, may 2018.
[82] Daniele Loco, Louis Lagardère, Stefano Caprasecca, Filippo Lipparini, Benedetta Mennucci, and Jean-Philip Piquemal. Hybrid QM/MM Molecular Dynamics with AMOEBA Polarizable Embedding. Journal of Chemical Theory and Computation, 13(9):4025–4033, sep 2017.
[83] Daniele Loco, Louis Lagardère, Gérardo A. Cisneros, Giovanni Scalmani, Michael Frisch, Filippo Lipparini, Benedetta Mennucci, and Jean-Philip Piquemal. Towards large scale hybrid QM/MM dynamics of complex systems with advanced point dipole polarizable embeddings. Chemical Science, 10(30):7200–7211, jul 2019.
[84] Daniele Loco, Étienne Polack, Stefano Caprasecca, Louis Lagardère, Filippo Lipparini, Jean-Philip Piquemal, and Benedetta Mennucci. A QM/MM Approach Using the AMOEBA Polarizable Embedding: From Ground State Energies to Electronic Excitations. Journal of Chemical Theory and Computation, 12(8):3654–3661, aug 2016.
[85] Pedro E M Lopes, Guillaume Lamoureux, and Alexander D Mackerell Jr. Polarizable empirical force field for nitrogen-containing heteroaromatic compounds based on the classical Drude oscillator. Journal of Computational Chemistry, 30(12):1821–1838, sep 2009.
[86] Pedro E M Lopes, Benoit Roux, and Alexander D MacKerell. Molecular modeling and dynamics studies with explicit inclusion of electronic polarizability: theory and applications. Theoretical Chemistry Accounts, 124(1):11–28, 2009.
[87] Rita P Magalhães, Henriques S Fernandes, and Sérgio F Sousa. Modelling Enzymatic Mechanisms with QM/MM Approaches: Current Status and Future Challenges. Israel Journal of Chemistry, 60(7):655–666, 2020.
[88] Yuezhi Mao, Yihan Shao, Jacek Dziedzic, Chris-Kriton Skylaris, Teresa Head-Gordon, and Martin Head-Gordon. Performance of the AMOEBA Water Model in the Vicinity of QM Solutes: A Diagnosis Using Energy Decomposition Analysis. Journal of chemical theory and computation, 13(5):1963–1979, may 2017.
[89] Billy W McCann and Orlando Acevedo. Pairwise Alternatives to Ewald Summation for Calculating Long-Range Electrostatics in Ionic Liquids. Journal of chemical theory and computation, 9(2):944–950, feb 2013.
[90] Gérald Monard, Michel Loos, Vincent Théry, Kristofor Baka, and Jean-Louis Rivail. Hybrid classical quantum force field for modeling very large molecules. International Journal of Quantum Chemistry, 58(2):153–159, jan 1996.
[91] Antonio Monari, Jean-Louis Rivail, and Xavier Assfeld. Theoretical Modeling of Large Molecular Systems. Advances in the Local Self Consistent Field Method for Mixed Quantum Mechanics/Molecular Mechanics Calculations. Accounts of Chemical Research, 46(2):596–603, feb 2013.
[92] Y Mauricio Muñoz-Muñoz, Gabriela Guevara-Carrion, Mario Llano-Restrepo, and Jadran Vrabec. Lennard-Jones force field parameters for cyclic alkanes from cyclopropane to cyclohexane. Fluid Phase Equilibria, 404:150–160, 2015.
[93] Kwangho Nam, Jiali Gao, and Darrin M York. An Efficient Linear-Scaling Ewald Method for Long-Range Electrostatic Interactions in Combined QM/MM Calculations. Journal of Chemical Theory and Computation, 1(1):2–13, jan 2005.
[94] Rebecca Notman and Jamshed Anwar. Breaching the skin barrier — Insights from molecular simulation of model membranes. Advanced Drug Delivery Reviews, 65(2):237–250, 2013.
[95] Pedro Ojeda-May and Kwangho Nam. Acceleration of Semiempirical QM/MM Methods through Message Passage Interface (MPI), Hybrid MPI/Open Multiprocessing, and Self-Consistent Field Accelerator Implementations. Journal of Chemical Theory and Computation, 13(8):3525–3536, aug 2017.
[96] Xiaoliang Pan, Edina Rosta, and Yihan Shao. Representation of the QM Subsystem for Long-Range Electrostatic Interaction in Non-Periodic Ab Initio QM/MM Calculations. Molecules, 23(10):2500, 2018.
[97] Jerry M. Parks, Hao Hu, Aron J. Cohen, and Weitao Yang. A pseudobond parametrization for improved electrostatics in quantum mechanical/molecular mechanical simulations of enzymes. The Journal of Chemical Physics, 129(15):154106, oct 2008.
[98] Robert M Parrish, Lori A Burns, Daniel G A Smith, Andrew C Simmonett, A Eugene DePrince, Edward G Hohenstein, Uğur Bozkaya, Alexander Yu. Sokolov, Roberto Di Remigio, Ryan M Richard, Jérôme F Gonthier, Andrew M James, Harley R McAlexander, Ashutosh Kumar, Masaaki Saitow, Xiao Wang, Benjamin P Pritchard, Prakash Verma, Henry F Schaefer, Konrad Patkowski, Rollin A King, Edward F Valeev, Francesco A Evangelista, Justin M Turney, T Daniel Crawford, and C David Sherrill. Psi4 1.1: An Open-Source Electronic Structure Program Emphasizing Automation, Advanced Libraries, and Interoperability. Journal of Chemical Theory and Computation, 13(7):3185–3197, jul 2017.
[99] Jean-Philip Piquemal, G Andrés Cisneros, Peter Reinhardt, Nohad Gresh, and Thomas A Darden. Towards a force field based on density fitting. The Journal of Chemical Physics, 124(10):104101, 2006.
[100] Steve Plimpton. Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of Computational Physics, 117(1):1–19, 1995.
[101] Jay W Ponder, Chuanjie Wu, Pengyu Ren, Vijay S Pande, John D Chodera, Michael J Schnieders, Imran Haque, David L Mobley, Daniel S Lambrecht, Robert A DiStasio, Martin Head-Gordon, Gary N I Clark, Margaret E Johnson, and Teresa Head-Gordon. Current Status of the AMOEBA Polarizable Force Field. The Journal of Physical Chemistry B, 114(8):2549–2564, mar 2010.
[102] Joshua A Rackers, Zhi Wang, Chao Lu, Marie L Laury, Louis Lagardère, Michael J Schnieders, Jean-Philip Piquemal, Pengyu Ren, and Jay W Ponder. Tinker 8: Software Tools for Molecular Design. Journal of chemical theory and computation, 14(10):5273–5289, oct 2018.
[103] A K Rappe, C J Casewit, K S Colwell, W A Goddard, and W M Skiff. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. Journal of the American Chemical Society, 114(25):10024–10035, dec 1992.
[104] Sandeep K Reddy, Shelby C Straight, Pushp Bajaj, C Huy Pham, Marc Riera, Daniel R Moberg, Miguel A Morales, Chris Knight, Andreas W Götz, and Francesco Paesani. On the accuracy of the MB-pol many-body potential for water: Interaction energies, vibrational frequencies, and classical thermodynamic and dynamical properties from clusters to liquid water and ice. The Journal of chemical physics, 145(19):194504, nov 2016.
[105] Demian Riccardi, Guohui Li, and Qiang Cui. Importance of van der Waals Interactions in QM/MM Simulations. The Journal of Physical Chemistry B, 108(20):6467–6478, may 2004.
[106] Anirban Sarkar, Sudipta Raha Roy, Naisargee Parikh, and Asit K. Chakraborti. Nonsolvent Application of Ionic Liquids: Organo-Catalysis by 1-Alkyl-3-methylimidazolium Cation Based Room-Temperature Ionic Liquids for Chemoselective N - tert -Butyloxycarbonylation of Amines and the Influence of the C-2 Hydrogen on Catalytic Efficiency. The Journal of Organic Chemistry, 76(17):7132–7140, sep 2011.
[107] Walter R P Scott, Philippe H Hünenberger, Ilario G Tironi, Alan E Mark, Salomon R Billeter, Jens Fennen, Andrew E Torda, Thomas Huber, Peter Krüger, and Wilfred F van Gunsteren. The GROMOS Biomolecular Simulation Program Package. The Journal of Physical Chemistry A, 103(19):3596–3607, may 1999.
[108] Hans Martin Senn and Walter Thiel. QM/MM Methods for Biomolecular Systems. Angewandte Chemie International Edition, 48(7):1198–1229, feb 2009.
[109] U. Chandra Singh and Peter A. Kollman. A combined ab initio quantum mechanical and molecular mechanical method for carrying out simulations on complex molecular systems: Applications to the CH3Cl + Cl- exchange reaction and gas phase protonation of polyethers. Journal of Computational Chemistry, 7(6):718–730, dec 1986.
[110] Marek Štrajbl, Gongyi Hong, and Arieh Warshel. Ab Initio QM/MM Simulation with Proper Sampling: “First Principle” Calculations of the Free Energy of the Autodissociation of Water in Aqueous Solution. The Journal of Physical Chemistry B, 106(51):13333–13343, dec 2002.
[111] Vincent Théry, Daniel Rinaldi, Jean-Louis Rivail, Bernard Maigret, and György G. Ferenczy. Quantum mechanical computations on very large molecular systems: The local self-consistent field method. Journal of Computational Chemistry, 15(3):269–282, mar 1994.
[112] Gregory S Tschumper, Matthew L Leininger, Brian C Hoffman, Edward F Valeev, Henry F Schaefer, and Martin Quack. Anchoring the water dimer potential energy surface with explicitly correlated computations and focal point analyses. The Journal of Chemical Physics, 116(2):690–701, dec 2002.
[113] M Valiev, E J Bylaska, N Govind, K Kowalski, T P Straatsma, H J J Van Dam, D Wang, J Nieplocha, E Apra, T L Windus, and W A de Jong. NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations. Computer Physics Communications, 181(9):1477–1489, 2010.
[114] Marc W van der Kamp and Adrian J Mulholland. Combined Quantum Mechanics/Molecular Mechanics (QM/MM) Methods in Computational Enzymology. Biochemistry, 52(16):2708–2728, apr 2013.
[115] Mary J Van Vleet, Alston J Misquitta, Anthony J Stone, and J R Schmidt. Beyond Born–Mayer: Improved Models for Short-Range Repulsion in ab Initio Force Fields. Journal of Chemical Theory and Computation, 12(8):3851–3870, aug 2016.
[116] Erik Vázquez-Montelongo, José Vázquez-Cervantes, and G. Cisneros. Polarizable ab initio QM/MM Study of the Reaction Mechanism of N-tert-Butyloxycarbonylation of Aniline in [EMIm][BF4]. Molecules, 23(11):2830, oct 2018.
[117] Claudia I Viquez Rojas and Lyudmila V Slipchenko. Exchange-repulsion in QM/EFP excitation energies – beyond polarizable embedding. Journal of Chemical Theory and Computation, aug 2020.
[118] Ross C Walker, Michael F Crowley, and David A Case. The implementation of a fast and accurate QM/MM potential method in Amber. Journal of computational chemistry, 29(7):1019–1031, may 2008.
[119] A Warshel and M Levitt. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. Journal of Molecular Biology, 103(2):227–249, 1976.
[120] Arieh Warshel and Robert M Weiss. An empirical valence bond approach for comparing reactions in solutions and in enzymes. Journal of the American Chemical Society, 102(20):6218–6226, sep 1980.
[121] Hiroshi C Watanabe and Qiang Cui. Quantitative Analysis of QM/MM Boundary Artifacts and Correction in Adaptive QM/MM Simulations. Journal of Chemical Theory and Computation, 15(7):3917–3928, jul 2019.
[122] Richard J Wheatley and Sarah L Price. An overlap model for estimating the anisotropy of repulsion. Molecular Physics, 69(3):507–533, feb 1990.
[123] Gaobo Xiao, Mingjun Ren, and Haibo Hong. 50 million atoms scale molecular dynamics modelling on a single consumer graphics card. Advances in Engineering Software, 124:66–72, 2018.
[124] Li Xie, Haiyan Liu, and Weitao Yang. Adapting the nudged elastic band method for determining minimum-energy paths of chemical reactions in enzymes. The Journal of Chemical Physics, 120(17):8039–8052, apr 2004.
[125] Yingkai Zhang. Improved pseudobonds for combined ab initio quantum mechanical/molecular mechanical methods. The Journal of Chemical Physics, 122(2):024114, jan 2005.
[126] Yingkai Zhang, Tai-Sung Lee, and Weitao Yang. A pseudobond approach to combining quantum mechanical and molecular mechanical methods. The Journal of Chemical Physics, 110(1):46–54, jan 1999.