Chapter 1: QM/MM with Auxiliary DFT in deMon2k¹

Over the last two decades, many quantum mechanical/molecular mechanical (QM/MM) implementations have emerged based on the pioneering works of Warshel and Karplus¹  and Warshel and Levitt²  published in the seventies. Common to these approaches is the possibility of describing a relatively small part of a complex system at a more or less sophisticated QM level of theory, whereas a computationally efficient MM force-field is employed for the remaining part of the system. The fundamental idea behind these QM/MM methodologies is to balance the accuracy and reliability of the simulation with the size of the underlying model. Naturally, the interface between the QM and MM regions is most delicate and defines to a large extent the achievable accuracy and reliability of the overall simulation. Therefore, it is common practice to extend the QM region as much as possible in order to diminish these interface effects. To this end, a promising first principles QM approach is Kohn–Sham density functional theory (DFT).^3,4  Although its computational demand is similar to the well-known Hartree–Fock method, it incorporates electron correlation through the exchange–correlation functional. In the linear combination of Gaussian-type orbital (LCGTO) approximation, Kohn–Sham DFT possesses a formal N$bas4$ scaling, with N_bas being the number of basis functions in the calculation, due to the expansion of the two-electron Coulomb repulsion into four-center electron repulsion integrals (ERIs). Already in very early molecular LCGTO-DFT Kohn–Sham implementations, the calculation of the four-center ERIs was identified as a critical computational bottleneck. To overcome this bottleneck, Dunlap and co-workers^5,6  introduced the variational fitting of the Coulomb potential. Due to its variational nature, this approximation also yields a variational energy expression free of four-center ERIs. Note, however, that the variational bound shifts below the corresponding four-center ERI Kohn–Sham bound. For relative energies, this shift in the variational bound has no effect. Thus, Kohn–Sham LCGTO-DFT approaches incorporating the variational fitting of the Coulomb potential yield the same relative energies as their four-center ERI counterparts at a considerably reduced computational cost. The same holds for so-called resolution-of-the-identity (RI) approaches⁷  because RI is a direct result of the variational fitting of the Coulomb potential and, therefore, possesses the same variational bound.⁸  To avoid four-center ERI calculations for hybrid functionals, the variational fitting of the Fock potentials^9–15  was introduced into deMon2k,¹⁶  too. The resulting four-center ERI-free Kohn–Sham energy model is one of the first principles QM approaches available in deMon2k for QM/MM simulations. It is named density fitted DFT (DF-DFT).

The computational bottleneck in medium-sized (50 to 500 atoms) DF-DFT calculations is the numerical integration of the exchange–correlation energy and its functional derivatives. Although this computational bottleneck was identified at the same time as that for the two-electron Coulomb repulsion, it has received much less attention. Early LCGTO-DFT implementations like deMon-KS¹⁷  or DGAUSS¹⁸  used the semi-numerical fit of the exchange–correlation potential with auxiliary functions as proposed by Sambe and Felton¹⁹  to circumvent this bottleneck. This approach gives reasonable energies but is prone to numerical instabilities in energy derivatives.²⁰  To overcome this problem the direct use of the (linear scaling) approximated density for the calculation of the exchange–correlation potential was suggested.^21–24  In the framework of deMon2k this has led to the development of auxiliary density functional theory (ADFT)²⁵  which is another first principles energy model available for deMon2k QM/MM simulations. As ADFT yields relative energies of the same quality as Kohn–Sham, however with a significant reduction in computational demand, we will focus in the following on ADFT as a first principles QM approach in QM/MM calculations. With the MPI parallelized ADFT implementation in deMon2k now available, large QM regions with 500 to 5000 atoms have become feasible. In these calculations, the computational bottleneck has shifted to linear algebra tasks such as matrix multiplications and matrix diagonalizations. In this respect, it is important to note that the variational fitting of the Coulomb potential and the use of an approximated density for the calculation of exchange–correlation contributions introduce additional linear algebra tasks in the form of solving inhomogeneous linear equation systems. To avoid corresponding computational bottlenecks we have adapted Krylov subspace methods, such as MINRES²⁶  and Eirola–Nevanlinna,²⁷  for the density fitting²⁸  and perturbation theory²⁹  calculations in deMon2k.

This chapter is organized as follows. In the next section, we derive the QM/MM energy and gradient working equations implemented in deMon2k. Section 1.3 describes algorithms for QM/MM structure optimizations and transition-state searches based on the trust-region method for local optimization. In the next section, the propagator equations for QM/MM molecular dynamics (MD) simulations in various closed-system ensembles using factorized Liouville operators are presented. Section 1.5 introduces auxiliary density perturbation theory (ADPT) for the calculation of analytic second-order ADFT energy derivatives. Their extension to QM/MM, currently under development in the deMon developers’ community, is discussed, too. The following section describes deMon2k magnetic shielding and excited-state calculations within the QM/MM framework. These are typical examples for QM/MM electronic property calculations employing perturbation theory. Section 1.7 presents the new QM/MM input builder (QIB) for deMon2k. A résumé of the current state-of-the-art ADFT QM/MM implementation in deMon2k is given in the last section of this chapter.

The QM LCGTO Kohn–Sham energy expression is given by:

In our notation Greek letters, here μ, ν, σ and τ denote (contracted) atom-centered Gaussian-type orbitals. For clarity of presentation, we restrict ourselves to closed-shell systems. The closed-shell density matrix elements, P_μν, are calculated as:

The upper limit in eqn (1.2), occ, denotes the number of doubly occupied molecular orbitals (MOs) in the QM system, and c_μi are the corresponding MO coefficients. For an isolated QM system, the first term in eqn (1.1) represents the mono-electronic core energy. The QM core matrix elements are given in atomic units, which we use throughout this chapter, by:

They collect the kinetic energy and nuclear attraction of the electrons. The second term in eqn (1.1) describes the Coulomb repulsion between the electrons. It includes the four-center two-electron integrals responsible for the N$bas4$ scaling of the method:

In this short-hand ERI notation ‖ represents the two-electron Coulomb operator 1/|r₁ − r₂|. It also separates functions from electron 1, on the left, from those of electron 2, on the right. In the following, we will use the same notation for two- and three-center ERIs. The last term in eqn (1.1) denotes the exchange–correlation energy functional evaluated with the Kohn–Sham density and its derivatives. In LCGTO the Kohn–Sham density can be expressed as:

As already mentioned, the two-electron Coulomb repulsion in eqn (1.1) introduces a formal N$bas4$ scaling due to the four-center ERIs. To circumvent this computational bottleneck, we employ the variational fitting of the Coulomb potential based on the following second-order energy error functional:^5,6 

The approximated electronic density, (r), also named auxiliary density, is calculated as a linear combination of auxiliary functions which we denote by Latin letters:

In deMon2k, atom-centered primitive Hermite Gaussian functions, indicated by a bar, are used as auxiliary functions.³⁰  Inserting this expansion of (r) and the LCGTO expansion of the Kohn–Sham density, eqn (1.5), into eqn (1.6) yields:

It is straightforward to show that ε$2H$ is semi-positive definite.³¹  Thus, we obtain from eqn (1.8):

Inserting this inequality into the Kohn–Sham energy expression, eqn (1.1), yields the QM DF-DFT energy expression:

The Coulomb fitting coefficients, x_k̄, in eqn (1.10) are calculated by minimizing the second-order error functional of eqn (1.6), which is accompanied by an energy maximization:³² 

Collecting these equations for all auxiliary functions yields the inhomogeneous linear equation system,

which must be solved in each self-consistent field (SCF) iteration. In eqn (1.12)G and J are the Coulomb matrix and vector, respectively, given by

To obtain the ADFT energy expression, only the functional argument of E_xc must be changed to the auxiliary density. Thus, the QM ADFT energy is given by:²² 

To also demonstrate that the ADFT QM energy is variational, we now calculate ADFT Kohn–Sham matrix elements as derivatives of the energy with respect to density matrix elements:

For the derivative of the ADFT exchange–correlation energy, E_xc[], assuming the local density approximation (LDA) for clarity of discussion, follows:

To proceed, we now define the ADFT exchange–correlation potential, v_xc[], as the functional derivative of the ADFT exchange–correlation energy with respect to the auxiliary density:

The differentiation of the auxiliary density in eqn (1.15) can be expressed with the formal solution of eqn (1.12),

Inserting eqn (1.16) and (1.17) into eqn (1.15) yields the derivative of the ADFT exchange–correlation energy:

Thus, the ADFT Kohn–Sham matrix elements can be straightforwardly calculated by the variation of the corresponding energy. This demonstrates the variational nature of the ADFT energy expression in eqn (1.13), albeit with its own variational bound. As for the DF-DFT energy, relative ADFT energies are in good agreement with their Kohn–Sham counterparts. Therefore, ADFT yields relative energies that are almost indistinguishable from their Kohn–Sham counterparts (≲ 0.1 kcal mol⁻¹) with significant reduction in the computational scaling and demand. For the explicit formulation of the ADFT Kohn–Sham matrix elements it is convenient to introduce the exchange–correlation coefficients as:

Eqn (1.19) can be reformulated as an inhomogeneous linear equation system analogous to the Coulomb fitting equation system in eqn (1.12):

In eqn (1.20)L denotes the exchange–correlation vector with elements L_k̄ = 〈k̄|ν_xc[]〉. Note that z is spin polarized. In the most recent deMon2k versions the linear equation systems of eqn (1.12) and (1.20) are solved by the Krylov subspace method MINRES.²⁸  As a result, the computational overhead associated with the calculation of Coulomb and exchange–correlation coefficients is negligible even for systems with hundreds of thousands of auxiliary functions. Thus, the use of extended auxiliary function sets such as GEN-A2*³³  presents no computational challenge. With the Coulomb and exchange–correlation coefficients at hand, the ADFT Kohn–Sham matrix can be expressed as:

Note that the ADFT Kohn–Sham matrix is free of density matrix elements, which greatly simplifies data handling in parallel applications. This is one reason for the excellent parallel scaling of ADFT.

Even though it is rather straightforward to extend the here outlined LDA ADFT discussion to functionals of the generalized gradient approximation (GGA), the computationally efficient implementation of hybrid functionals into ADFT (and DF-DFT) is less obvious. To facilitate our discussion, we use the following generic ADFT hybrid functional QM energy expression:

Here α represents the mixing parameter for the Fock exchange energy, E$xF$, of the hybrid functional. The corresponding counterpart (1 − α) for the DFT exchange energy is absorbed into the ADFT exchange–correlation energy functional and will be of no further concern to the following discussion. For a more detailed discussion of the structure of global and range-separated hybrid functionals in ADFT, we refer the interested reader to ref. 15 and 34–36. The computational drawback arises from the Fock exchange energy,

For the variational Fock potential fitting we now introduce MO transition densities, ρ_ij(r), and corresponding approximated transition densities, _ij(r):

With these MO transition densities the variational fitting of the Fock potential can be based on the following second-order error functional:^13,15 

Expanding the transition densities in eqn (1.26) yields:

The Fock exchange fitting coefficients are obtained by the maximization of the negative fitting error,

which is accompanied by the minimization of the approximated Fock exchange energy. This suggests that the simultaneous variational fitting of Coulomb and Fock potentials will benefit from systematic error compensation if the same auxiliary function set is used for both fits. In fact, if the variational nature of both potential fits is preserved, the resulting four-center ERI-free total Hartree–Fock energies for molecules with thousands of basis functions are nearly indistinguishable (< 1 kcal mol⁻¹) from their four-center ERI counterparts! Corresponding binding energy differences between four- and three-center Hartree–Fock calculations are below 10⁻² kcal mol⁻¹. From eqn (1.28) the following linear equation systems,

with

are obtained. Therefore, the Fock exchange-fitting coefficients are given by:

Inserting these Fock exchange-fitting coefficients into eqn (1.27) and the resulting approximation for the Fock exchange energy into eqn (1.22) yields the following generic ADFT hybrid functional energy expression free of four-center ERIs:

Direct implementation of eqn (1.32) leads to an algorithm that scales as (N_aux × N$bas2$ × N_occ) with N_aux and N_occ being the number of auxiliary functions and occupied MOs, respectively. Such an algorithm is useful only when (N_bas) ≫ (N_occ).^10,11  Note, however, that eqn (1.32) is invariant under orthogonal transformation of the MOs, because the density matrix is invariant under such transformations. Thus, any set of MOs obtained by an orthogonal transformation of the canonical MOs (CMOs) can be used in eqn (1.32). In particular, the CMOs (see Figure 1.1a) can be transformed into spatially localized MOs (LMOs) by minimizing or maximizing, an appropriate functional (Figure 1.1b). This facilitates efficient screening in the integral transformation from the atomic to the molecular representation.

The recently developed variational fitting of the Fock potential¹³  based on Foster–Boys localization explores this possibility by utilizing LMOs in eqn (1.32). As a result, it shows an improved computational performance for the calculation of variationally fitted Fock exchange energies. In this implementation, a two-step localization is used. First, an incomplete Cholesky decomposition³⁷  is performed followed by a tighter Foster–Boys localization.^38,39  The MO localization permits a variational fitting of Fock exchange by defining the fitting domains around each LMO. In particular, we define the fitting domains around each localized MO in terms of atomic centers. To this end, we calculate the atomic Löwdin population for each atom A in a given LMO, ψ′_i, according to

Note that in QM/MM calculations, the core Hamilton matrix elements in eqn (1.33) are augmented by the electrostatic embedding from the MM region. They are given by:

In eqn (1.34)H$μνQM$ are the QM core matrix elements given by eqn (1.3), and Q_D are the MM atomic charges. For generality of discussion, we introduce the nuclear attraction type operator defined as:⁴¹ 

This general definition allows inclusion of higher point moments on MM atoms, e.g. for polarizable force-fields.⁴²  Work in this direction is currently under development within the deMon developers’ community. However, for clarity of presentation we restrict ourselves here to MM point charges, i.e. k_x = k_y = k_z = 0 in eqn (1.35). Although the core Hamilton matrix elements are only calculated once for each geometry at the beginning of the SCF, they can become computationally demanding for larger numbers of MM atoms if standard nuclear attraction integral recurrence relations⁴³  are used. Due to the slow decay of Coulomb interactions an efficient screening of these integrals is not possible. However, it is possible to separate the QM/MM electrostatic interaction integrals in eqn (1.34) into near- and far-field integrals.⁴⁴  The far-field QM/MM electrostatic interaction integrals are obtained by the following asymptotic expansion:⁴⁵ 

with

In eqn (1.36) the sum limits Near and Far refer to MM atoms in the near- and far-field of the QM region. There are two major advantages associated with the above asymptotic expansion of the electrostatic MM embedding. First, nuclear attraction type integrals are substituted by overlap type integrals, 〈α + m|β〉, which avoid the calculation of the incomplete gamma function. Second, the far-field loop over MM atoms is factored out and, therefore, scales only with the number of MM atoms independent of the electronic expansion of the QM region. As a result, the MM embedding in deMon2k QM/MM SCF calculations produces only a small computational overhead in the range of a few percent compared to the corresponding QM calculation.⁴⁴  The QM/MM modified core matrix elements enter the corresponding ADFT Kohn–Sham matrix according to eqn (1.21) where H$μνQM$ is substituted by H_μν from eqn (1.34). In this way, the MM embedding directly alters the SCF.

In practical QM/MM applications the user has to choose the density functional approximation (DFA), the basis set, the auxiliary functions set and other specifications in the methodological frameworks of DF-DFT and ADFT. The DFAs available for QM/MM calculations in deMon2k are LDA, GGA and meta-GGA (only within DF-DFT) as well as global and range-separated hybrids. By default, the LDA is used in the form of Dirac exchange⁴⁷  and VWN correlation.⁴⁸  Because, in ADFT, LDA and GGA performance is similar, we recommend GGA functionals for QM/MM calculations. In particular, PBE⁴⁹  in combination with the GGA optimized DZVP-GGA basis and GEN-A2* auxiliary function set³³  yields reliable optimized structure parameters and fair qualitative binding energies throughout the periodic table. Note that for heavier elements scalar relativistic effects must be taken into account either by effective core potentials (ECPs) or by model core potentials (MCPs).⁵⁰  The ECP approach can also be used for the introduction of capping potentials⁵¹  in QM/MM calculations. We also recommend the PBE functional for QM/MM calculations because analytic ADFT kernel formulas⁵²  are implemented in deMon2k. This greatly facilitates QM/MM response calculations with ADPT. For most other GGA functionals the ADFT kernels are calculated by finite differences.⁵³  If van der Waals interactions are important, e.g. in QM/MM structure optimizations, the DF-DFT and ADFT LDA or GGA energy expressions must be augmented with empirical dispersion terms^54–56  because these interactions are not captured by those functionals. Although ADFT hybrid calculations show similar serial performance as corresponding GGA runs, we do not recommend them for initial QM/MM structure optimizations due to their less-than-ideal parallel scaling.⁴⁰  Instead, we advocate composite approaches⁵⁷  consisting of GGA structure optimization and single-point hybrid energy calculations, e.g. with PBE0.^58,59  At this point, it is important to note that basis set convergence for hybrid functionals is much slower than for GGA functionals. To this end, we have had good experience with augmented Dunning-type basis sets.⁶⁰  The handling of the three-center ERIs and the grid accuracy for the numerical integration of the exchange–correlation energy and potential are important for the computational performance of QM/MM calculations, too. By default, the direct SCF setting is used in deMon2k. In this approach the three-center ERIs are repeatedly calculated (twice) in each SCF iteration. If sufficient random access memory (RAM) is available, the use of the so-called mixed SCF approach^61,62  is very beneficial for the computational performance. Tightening of the adaptive grid accuracy^63,64  beyond the default value of 10⁻⁵ a.u. will increase computational demand. Thus, care is advised. For linear algebra tasks, we recommend the use of MINRES for the density fitting and the ScaLAPACK extension of deMon2k for matrix diagonalization.

With QM energy models for QM/MM calculations at hand, we now express the QM/MM energy as:

In eqn (1.38) E^MM, E^QMMM and E^NN denote the MM energy, the QMMM interaction energy and the nuclear repulsion energy for the QM/MM system, respectively. All three terms are classical and, therefore, have no dependency on the wavefunction or the electronic density. The current default QM/MM methodology implemented in deMon2k relies on first generation, non-polarizable, force-fields according to which the MM energy expression reads:

The first five terms correspond to bonded contributions, whereas the last two are non-bonded terms. The extension to a polarizable QM/MMpol methodology that further includes electrostatic induction within the framework of the point-charge dipole model will be discussed in Chapter 4. For the simplicity of our discussion here, we will use generic formulas for the individual terms, i.e. for bond stretching, Urey–Bradley interaction, angle bending and so on. The analytic forms of the generic bond stretching, Urey–Bradley interaction, angle bending and proper torsion are given by:

In the first equation, the distance is calculated as d = |A − B|. The equilibrium distance d₀ and the force constant k_b are force-field parameters for the specific bond. In the Urey term (see Figure 1.2), l represents the distance between two atoms, with the difference that these atoms are at the ends of an angle. The equilibrium distance l₀ and force constant k_l are given by the force-field, too.

The angle bending in eqn (1.42) is modeled by a harmonic potential as well. The corresponding force constant, k_ϑ, and equilibrium angle, ϑ₀, are force-field parameters for the bending. Eqn (1.43) is used for the calculation of both proper and improper torsions as depicted in Figure 1.2. The difference between proper and improper torsions is not in their mathematical formulation, but rather in how they are defined: proper torsions are defined over successively connected atoms, such that whenever going from one atom to the following, there is always a bond between them. On the other hand, for improper torsions this connectivity is not mandatory. This serves various purposes like restricting three atoms into a plane, e.g. to describe an sp² carbon binding. Figure 1.2 shows, besides the Urey interactions at the top, typical torsion definitions at the bottom. Similar to the rest of the bonded terms the parameters needed for the calculations of the torsion contributions are given by the force-field dataset (FFDS) file. The number of expansion coefficients, V_n, and phase factors, φ_n, in eqn (1.43) is four for proper torsions and three for improper torsions. All bonded MM interactions require molecular connectivity information, which is either provided by the user in the deMon2k input file⁶⁵  (see also Section 1.7) or automatically generated based on the distances between MM atoms. The latter is particularly useful for molecular assemblies, e.g. water clusters.

For the non-bonded interactions, the energy expressions used are the well-known van der Waals (aka Lennard-Jones) and point charge Coulomb potentials:

As for the bonded terms, the atomic parameters for atoms A and B in eqn (1.44) and (1.45) are defined by the force-field. To obtain the diatomic parameters ϵ_AB and σ_AB, a combination rule, which can be either an arithmetic or a geometric average, is used as it is shown below for the epsilon parameter:

The last two terms of eqn (1.38) are the mechanical QM/MM interaction and the nuclear repulsion. Both have simple analytic formulas in terms of QM and MM atomic positions. The QM/MM interaction energy, E^QMMM, is modeled in deMon2k by a Lennard-Jones potential:

This is the same analytic expression as in eqn (1.44) for the van der Waals interactions between MM atoms. The difference is that in eqn (1.46) one of the atoms is a QM atom and the other is a MM atom. The diatomic parameters ϵ_AD and σ_AD are again calculated by the above-described combination rules from atomic parameters. To this end, MM atomic ϵ and σ parameters are assigned to the QM atoms by the QM/MM keyword⁶⁵  in deMon2k. The last term in eqn (1.38) is the electrostatic repulsion between the QM point-charge nuclei with themselves and with the MM atomic point charges. It has the following form:

The MM atomic point charges, Q_D, are again taken from the FFDS file. Figure 1.3 (top) depicts the energy modules of deMon2k that can contribute to QM/MM calculations. Once the QM/MM energy is calculated, which implies a converged SCF solution for the QM part, the gradient components are calculated by differentiation of eqn (1.38) with respect to the atomic coordinates A_λ, λ = x, y, z:

As eqn (1.48) shows, all energy terms of a QM/MM calculation contribute to the corresponding gradient. From the gradient on, QM and MM atoms are treated on an equal footing in deMon2k. Thus, the following utility modules for structure optimization and chemical reaction modeling (see Section 1.3) as well as for MD simulations (see Section 1.4) are the same for QM, MM and QM/MM calculations. On the other hand, the analytic higher-order SCF energy derivatives calculated in the TD-ADFT and ADPT modules of Figure 1.3 are directly linked to the ADFT module and, therefore, are only available for the QM part of a QM/MM calculation (see Section 1.6). The QM/MM frequency analysis is linked to both the ADPT and MM module as outlined in Section 1.5.

The first term in eqn (1.48) is the generic QM energy derivative. For ADFT^15,22,35,36  QM/MM calculations this term is obtained by differentiation of eqn (1.33) with respect to a QM atomic coordinate:

The superscripts denote differentiation of molecular integrals or functions with respect to the QM atomic coordinate. Eqn (1.49) consists of three contributions. The first one is the Pulay force⁶⁸  given by the first term of eqn (1.49). It collects all derivatives of the MO coefficients. Thus, no response calculation is needed for the gradient evaluation as expected from Wigner's 2n + 1 rule.^69,70  The second contribution arises from the QM/MM core (second term in eqn (1.49)), Coulomb (third and fourth terms in eqn (1.49)) and exchange–correlation (fifth term in eqn (1.49)) derivatives. Together with the Pulay forces, these contributions define the QM ADFT energy derivatives for LDA and GGA. The last two terms in eqn (1.49) define the third contribution given by the derivative of the three-center ERI Fock exchange energy. The closed-shell energy-weighted density-matrix elements appearing in the Pulay term are given by:

Here ε_i refers to the canonical MO energies as obtained from the converged SCF. The Γ_k̄l̄ and x_k̄μi coefficients are calculated as:

The QM energy inside the QM/MM energy expression, eqn (1.38), depends, through the core matrix elements defined in eqn (1.34), on the atomic coordinates of the MM atoms, too. Therefore, the corresponding gradient component for a MM atom coordinate D_λ is given by:

We note that for these derivative integrals (electric field integrals), asymptotic expansions³¹  analogous to the core-matrix elements are available in deMon2k, too. The DF-DFT gradients for QM/MM calculations differ from the ADFT gradients in the term for the exchange–correlation in eqn (1.49).

The second term in eqn (1.48) is the MM energy derivative. It only depends on the atomic coordinates of the MM atoms. As for the MM energy, the MM derivatives can be partitioned into bonded and non-bonded derivatives. For the derivatives of the generic bond stretching and Urey–Bradley terms with respect to the MM atomic coordinate A_λ follows:

From the corresponding differentiation of the angle bending energy, we find:

Here A_λ and C_λ refer to atomic coordinates of the end atoms of the A − B − C angle. The angle derivative for the central B atom is calculated from the translation invariance as:

Due to their definitions, both proper and improper torsions share the same differentiation formulas, and thus we write the equations referring to the proper torsions only. We start by writing the gradient of eqn (1.43) in terms of the derivative of the dihedral angle φ:

As for the case of the bending energy gradients, the derivatives of the dihedral angle depend on the relative position of the atom in the dihedral. The four different cases are:

The vectors P and Q are perpendicular to the planes (ABC and BCD, respectively) that define the dihedral angle φ. They are given by the following cross products:

Finally, the energy derivatives of the non-bonded MM interactions are comprised by the gradients of the van der Waals and Coulomb potentials. These are given by:

For the QM/MM mechanical interaction, we find similar results to eqn (1.65) as the same model is used for both the E^QMMM and the E^vdW energy terms. Nonetheless, for E^QMMM we have derivatives with respect to both QM and MM atomic coordinates:

In these gradients, the QM and MM atomic coordinates are given by A_λ and D_λ, respectively. The same convention is used for the gradients of the repulsion between the QM point-charge nuclei themselves and with the MM atomic point charges, eqn (1.47). These gradients are given by:

With these final equations, all terms needed for the calculation of the QM/MM energy gradients are at hand for the different methodologies implemented in deMon2k, as depicted in Figure 1.3.

The potential energy surface (PES) of a molecular system is a multidimensional function arising from the Born–Oppenheimer approximation⁷¹  that possesses different kinds of stationary points. At these points, the first energy derivatives with respect to the nuclear coordinates, i.e. the energy gradients, are zero. For a molecule with M atoms the dimension of the PES is given by its internal degrees of freedom, which is either 3M − 5 or 3M − 6 depending on the molecule being linear or not, respectively. The different types of stationary points are indicated on the PES model illustrated in Figure 1.4.

As shown in Figure 1.4, different stationary points such as minima, maxima and saddle points can exist on a PES. The number of stationary points of the PES increases exponentially as the number of atoms increases. Fortunately, from a chemical point of view only minima and first order saddle points are of interest because they correspond to stable minimum structures and transition states, respectively. Nevertheless, the exploration of the full PES remains a formidable task for many molecular systems.

On the other hand, the PES can be efficiently explored in local regions employing numerical optimization techniques. For quantum chemical structure minimization, quasi-Newton methods^72–84  are commonly applied. They are particularly efficient for QM methods because only information on energy and gradients are needed. Because in deMon2k QM/MM calculations the QM and MM gradients are joined together (see Figure 1.3), a quasi-Newton trust region method (TRM)^72,73  is employed for local QM/MM structure optimizations by default. To this end, the PES is expanded around a reference molecular structure by a second-order Taylor series. This quadratic PES expansion is given by:

In eqn (1.71)g is the QM/MM gradient vector and B the corresponding approximation of the Hessian matrix. Both quantities are evaluated at the reference structure with energy E₀. In order to find a minimum, the quadratic expansion is optimized with respect to the step direction p. To ensure that the model function q(p) is a good approximation of the PES, we restrict the step size within a trust region defined by a hypersurface with radius h around the reference structure. Thus, the following constrained-minimization problem is obtained:

For the solution of eqn (1.72) the following Lagrange function is introduced:⁸⁵ 

From the stationary condition of this Lagrange function

the step direction p is defined as:

For the determination of the Lagrange multiplier λ that multiplies the identity matrix in eqn (1.75), standard techniques from the literature^73,85–87  are used. Similarly, to update the (approximated) Hessian matrix the BFGS update formula, named after its inventors Broyden,⁸⁸  Fletcher,⁸⁹  Goldfarb⁹⁰  and Shanno,⁹¹  is used:

In eqn (1.76)γ = g^{I + 1} − g^I and δ = x^{I + 1} − x^I hold. The superscripts I and I + 1 refer to optimization steps. Thus, QM/MM structure optimizations in deMon2k follow the same path as pure QM structure optimizations. The advantage is the robustness of the underlying algorithms. The disadvantage of the TRM in particular is the need for the diagonalization of the (approximate) Hessian matrix. This limits its applicability to QM/MM system sizes with up to around 10 000 atoms. For larger systems, alternative step methods free of second-order energy derivatives such as steepest descent (see keyword STEPTYPE⁶⁵ ) are available in deMon2k, too.

Besides the excellent performance of the TRM for molecular structure optimizations, it is also well suited for the search of transition states on the PES. In general, a transition state search is more complicated than the structure minimization just discussed due to the (first-order) saddle point nature (see Figure 1.4) of the transition states depicted in Figure 1.5.

However, combining the TRM with eigenvector following in the so-called uphill trust region method^92–94  yields a robust local transition state search algorithm. The corresponding step formula is given by:

In eqn (1.77)M is the so-called mode matrix defined as:

The position of the negative entry in the mode matrix M defines the Hessian matrix mode, which is followed uphill. From a computational point of view, it is desirable that the initial Hessian matrix for the transition state search has the right eigenvalue spectrum, i.e. one negative eigenvalue and all others positive. This can be achieved by an appropriate start structure and the explicit calculation of the initial Hessian matrix. To this end, the hierarchical transition state search algorithm⁹⁴  is currently being extended for QM/MM calculations in one of our laboratories. In transition state searches the unbiased Powell-symmetric-Broyden⁹⁵  update formula is used in deMon2k by default:

Whereas the QM/MM molecular structure optimizations in deMon2k utilize algorithms from QM calculations, the QM/MM MD simulations use MM-inspired algorithms. To this end, we assume classical equations of motion for the QM nuclei and MM atoms and introduce the concepts of ensembles for closed systems. Based on this approximation, microcanonical (NVE) and canonical (NVT) QM/MM MD simulations can be performed with deMon2k. The implementation of isothermal–isobaric (NPT) QM/MM MD simulations is currently under development in the deMon developers’ community.

For the conservation of temperature and pressure during MD simulations, it is necessary to include thermostats and barostats that can regulate these properties. For this purpose, the reformulation of the equations of motion in terms of the Liouville operator⁹⁶  simplifies the application and implementation of the propagation equations. We start our discussion by writing the Newton equations for a system of M particles consisting of QM nuclei and MM atoms:

In these equations, the index A addresses QM nuclei and MM atoms with masses m_A. Their positions and momenta are labeled by the symbols R_A and P_A. Dots indicate the corresponding time derivatives. We now rewrite eqn (1.80) and (1.81) as:

The Liouville operator used in these expressions is given by:

A more compact and general form of this operator is found if we collect the positions and momenta of all particles in the system into a single variable x = (R^M, P^M), known as the phase-space vector. With this vector, eqn (1.84) can be written as iL = ẋ · ∇_x, where ∇_x comprises the differentiation operators from all variables in x. Therefore, the evolution in time of a system is given by:

The exponential term in this equation is known as the classical propagator. Its form depends on the equations governing the progress of the system and, as such, the analytic solution of eqn (1.85) becomes very difficult as the system grows in complexity. Therefore, a numerical scheme is used for its implementation, namely the Trotter factorization⁹⁷  of the classical propagator:

Here, the operator of eqn (1.84) has been arbitrarily separated into two terms, iL₁ and iL₂, which propagate, individually, positions and momenta. Their forms are:

As a result, the propagator for NVE simulations reads:

The application of this operator to the phase-space vector of a system results in the well-known velocity Verlet algorithm.⁹⁸  The advantage of using this formulation resides in the direct translation between the resulting working equations and the programming statements.

For the generation of the NVT ensemble, a variety of different thermostats⁹⁹  are available in deMon2k. All of them work by controlling the kinetic energy fluctuations of the system by changing the particle velocities. The extended phase-space approach of the Nosé–Hoover chain (NHC)^100,101  thermostat, by which additional variables are included into the simulation, provides a methodology to simulate the canonical ensemble using the following equations of motion:

In these equations the k index loops over the N thermostats, each having a unitless position η_k, a mass Q_k with units of energy times time squared, and momentum P_ηk in units of energy times time. Only the first link of the thermostat chain couples to the system as indicated by eqn (1.91). The k ≥ 2 thermostats regulate the variables of the previous link in the chain. The forces driving the evolution of the thermostats’ dynamical variables are given by:

Here, k_B is the Boltzmann constant and T the desired macroscopic temperature. The resulting propagator builds upon that of eqn (1.89), with the difference that the new factorization will include a third term iL_T accounting for the presence of the thermostat variables η_k and P_ηk. This third term has the form:

The corresponding propagator in its factorized form is then given by:

This is essentially the same propagator as for the NVE case, the only difference being the update of the thermostat variables before and after the application of the velocity Verlet algorithm. The NHC algorithm within the formalism of the Liouville operator is implemented in the current development version of the deMon2k software.

Within the NVE and NVT ensembles, the dynamics of solute–solvent systems can be effectively performed using QM/MM, capped by a continuum solvent model (PCM). The PCM model is based on the Onsager approximation, implemented in the self-consistent reaction field DFT method.^67,102–105  Within the QM/MM scheme, the reaction field energy reads:

In eqn (1.99)g is the response function of the dielectric solvent medium with respect to the potential of the QM/MM system. For a spherical cavity with radius a, and the solvent medium with a dielectric constant ε, the response function is:

The reaction field energy is added to the QM/MM energy in eqn (1.38), and the gradient of the reaction field with respect to the atomic coordinates is added to the QM/MM energy gradient in eqn (1.48).

For QM/MM simulations at constant temperature and pressure, the Martyna, Tobias, Klein (MTK)^106,107  barostat is currently extended to QM/MM applications in deMon2k. In this algorithm, the volume of the system is introduced as a dynamical variable together with the corresponding momentum. When used in combination with the NHC thermostat, the isobaric–isothermic NPT ensemble can be sampled. The equations of motion for the particle system and barostat are:

The terms labeled with η belong to the thermostat coupled to the particles, whereas those labeled with ξ correspond to the thermostat coupled to the barostat, whose variables are labeled with ϵ. Eqn (1.92) to (1.94) describe the dynamics of both thermostats, although the ξ thermostat is subject to the forces:

The force driving the barostat, G_ϵ, is given by:

The external pressure, P_ext, is the target value for the simulation. The corresponding propagator resembles that of eqn (1.98) with the addition of the terms from the barostat and the second thermostat. However, the presence of the barostat changes how the positions and velocities of the particles are updated within the (modified) velocity Verlet algorithm. It now reads:

Here α = 1 + 1/M. The changes in the cell containing the system are considered isotropic and retrieved from the instantaneous volume.

From this discussion, it is obvious that a volume definition is needed for the simulation of systems under constant pressure. To this end, we introduce periodic boundary conditions in the form of the cyclic cluster model (CCM).^108–110  In this model, the molecules are placed inside a (large) unit cell, which defines the size of the system during the calculation, e.g. the number of atoms. This unit cell is then translated along its axes in all directions without increasing the system size for the calculation. Such a scheme is shown in Figure 1.6 for a two-dimensional unit cell containing the atoms A0, B0 and C0. The dashed objects are created by translation of the central unit cell.

In the CCM framework, the interaction between atoms A and B, for example, may not be calculated between A0 and B0. In the case of the scheme shown in Figure 1.6, the calculated interaction will be between A0 and B4. The atom B4 is chosen because it has, of all the B atoms, the smallest distance to atom A0. On the other hand, the calculated interaction for atoms A and C is between A0 and C0 because, of all the C atoms, C0 has the smallest distance to atom A0. If there are several combinations that have the same distance, all are calculated and an average is used. As in a finite system, the interactions are calculated in real space. What is different is the relative orientation of the atoms to each other. For a structure like a regular lattice of argon atoms, this procedure creates for each atom an identical and symmetric environment. Therefore, the CCM scheme is an alternative realization of periodic boundary conditions. Figure 1.7 compares the convergence of the cohesive energy of solid Lennard-Jones argon with increasing cluster size between the free cluster model (FCM) and CCM. As this figure shows, CCM clusters reach the converged cohesive energy of −8.42 kJ mol⁻¹ at much smaller cluster sizes, i.e. at smaller average coordination numbers, than FCM clusters.

Another characteristic of CCM is that the atoms inside the system cannot leave it due to the equality between atoms and their translations. Thus, if during a Born–Oppenheimer molecular dynamics (BOMD) simulation an atom leaves the unit cell, a translation will bring it back inside. Because the projection of the distance vector between two atoms onto an edge of the unit cell is not allowed to exceed half of the length of that edge, this effectively creates an interaction region of this size around each atom. Only atoms within this region generate interactions with the atom in the center. The shape of this region is equal to a Wigner–Seitz cell (WSC) constructed from the unit cell. With this, one can describe the CCM as exchanging summations over all atoms in the (large) unit cell by summations over all atoms in a WSC. A further consequence of the use of WSCs is that there are no direct interactions of an atom with its equivalent in a translated unit cell, meaning that direct defect–defect interactions do not exist. This is particularly useful in a QM/MM solvation model for the treatment of infinitely dilute systems in which the solute is seen as a defect in the periodicity of the solvent environment.

In order to extend CCM to polar systems, long-range Coulomb interactions must be added. To this end, the lattice sum algorithm of Ewald in its classical¹¹¹  formulation is available in deMon2k. The target quantity for this method is the electrostatic potential generated by the classical point charges of the MM system given by:

For the sake of simplicity, it has been assumed that the particles are inside a cubic box of length L. The vector m of integer components works, together with L, as a translation operator of the reference cell. Its norm is represented by |m|. Care must be taken in eqn (1.110) for m = (0, 0, 0) to avoid self-interaction. In the classical Ewald method, the calculation of eqn (1.110) for neutral systems is carried out by means of two rapidly and absolutely convergent series evaluated in real and reciprocal space.¹¹²  The starting point of this method is the total potential generated by the point charges of the particles inside the original cell and its images, namely:

The main idea behind the Ewald method is to use Gaussian charge distributions of the same magnitude but opposite sign of the point charges to shield them from their environment, as shown in Figure 1.8 (left). These distributions are of the form:

The parameter α is the width of the distributions. It also has the effect of biasing the total potential energy toward one of the two main contributions into which it is divided. Of course, this energy is invariant with respect to the value of α.

Due to the presence of the Gaussian charge distributions, the point-charge potential becomes short-range and can be calculated directly in real space. On the other hand, a second set of Gaussian distributions with the same sign and magnitude as the original point-charge distributions is added to cancel the first set of distributions. Thus, the potential from the point charges is split into two parts, and the energy expression of eqn (1.110) becomes:

The first term in this equation is a sum in the manner of the Coulomb potential calculated over the original and image cells (m vectors), weighted by the complementary error function . This function approaches zero as its argument grows, improving the convergence of the sum. The second term in eqn (1.113) is a sum arising from the compensating distributions and is calculated over the reciprocal vectors . The last term in eqn (1.113) accounts for self-interaction.

In practice, finite limits for the m and k sums are needed. These limits are chosen as a compromise of the desired accuracy and the computational demand for the calculation. Usually, the sum over m is carried out in the original cell only (|m| = 0), whereas the k sum is performed over about 300 vectors.^112,113  This is achieved by tuning the α parameter toward large values. For ,¹¹⁴  both sums converge at the same rate. A value of α = 6/L is used in deMon2k. To illustrate the performance of this implementation we compare MM CCM water calculations augmented by long-range Coulomb lattice sums according to eqn (1.113) against deMon2k embedding and periodic Tinker¹¹⁵  calculations. To this end, a large unit cell (CCM cluster) of four water molecules as depicted in Figure 1.9 (left) is used. For the deMon2k embedding calculations performed with eqn (1.110), the system was replicated around the original cell following both a cubic and a spherical pattern (Figure 1.9, right). The convergence of the energy of these systems toward the value calculated with eqn (1.113) of −22.40 kcal mol⁻¹ is shown in Figure 1.10 (left). The norms of the corresponding forces per atom for this system are shown in the right part of Figure 1.10.

The Ewald method as described above is implemented in deMon2k for MM systems under CCM periodic boundary conditions. Nonetheless, in the case of QM/MM systems, interactions between the QM electron density and the periodic point charges need to be considered. This same problem has been recently investigated in the literature,¹¹⁶  however, within a different working framework to the one discussed here. In QM/MM systems the potential energy has the following form:

The potential ϕ(r) is given by eqn (1.111). As it was commented before, this potential is split into two contributions by the Ewald method. Therefore, eqn (1.114) is partitioned into two integrals. Their forms, when considering that the real part of the potential is calculated in the original cell only, are:

Notice that there is no need for a self-interaction term, which only appears in eqn (1.113) due to the point-like nature of the MM charges. An efficient methodology for the evaluation of these integrals is currently under development in one of our research groups. The availability of these integrals and their gradients will make it possible to use the MTK barostat inside QM/MM calculations and thus, fully atomistic simulations under ambient conditions will become possible.

The calculation of second-order energy derivatives with respect to nuclear displacements in the form of Hessian matrix elements is essential for the characterization of optimized molecular structures. It is also of great importance for transition state searches (see Section 1.3) and structure optimizations with soft internal degrees of freedom that are common in QM/MM systems. The dimensionality of the Hessian matrix is 9M², with M being the number of atoms in the system. Although the first-order electronic energy derivatives are relatively easy to calculate, second-order derivatives are less straightforward. The simplest way to calculate second-order derivatives is to follow a finite difference method where typically two first-order derivatives are required for each atomic coordinate. This method is available for QM/MM calculations in deMon2k but requires significant computational resources and is prone to numerical instabilities. Thus, an analytic second-order energy derivative formulation for QM/MM is very desirable.

Similar to the calculation of the QM/MM gradient components, after a converged SCF solution for the QM part, Hessian matrix components of a molecular system are obtained by differentiation of eqn (1.48) once again with respect to atomic coordinates B_η, with η = x, y, z:

The first term in eqn (1.117) corresponds to the generic QM second-order derivative. For an ADFT QM-only calculation, this term is calculated by differentiation of eqn (1.13), with respect to a QM atom coordinate A_λ and a QM atom coordinate B_η, and the second, third and fourth term of eqn (1.117) will vanish. According to Wigner's 2n + 1 rule,^69,70  the resulting expression contains the first-order response density matrix as well as the response Coulomb and exchange–correlation fitting coefficients. Based on this, the expression of the second-order ADFT energy derivatives can be partitioned into two contributions as:¹¹⁷ 

The first contribution of eqn (1.118) is completely independent of the perturbed density matrix and fitting coefficients and has the form:

Here, f_xc[] denotes the ADFT exchange–correlation kernel, defined as the functional derivative of the ADFT exchange–correlation potential. In the case of pure density functionals, the exchange–correlation kernel is defined as:^52,118 

The analytic molecular integral derivatives in eqn (1.119) are calculated over the corresponding integral recurrence relations,^30,119  whereas numerical integration of the exchange–correlation potential and kernel are evaluated with the same adaptive grid that is used for the SCF.^63,64  The exchange–correlation kernel is calculated either from analytic expressions for each functional⁵²  or by finite differences if only analytic exchange–correlation potential expressions are available.⁵³  The density, P_μν, and energy-weighted density, W_μν, matrix elements are calculated from the converged canonical MO coefficients and energies. The Coulomb fitting coefficients, x_k̄, and exchange–correlation fitting coefficients, z_k̄, are also taken directly from the previously converged single-point energy calculation. The second contribution to the QM second-order ADFT energy derivative in eqn (1.118) contains only terms that depend on first-order response matrix or vector elements. It is given by:

Again, the molecular integral derivatives can be straightforwardly calculated by integral recurrence relations. For the calculation of the perturbed matrix and vector elements, we employ auxiliary density perturbation theory (ADPT). This methodology is an adaptation of McWeeny's self-consistent perturbation (SCP) theory^120–124  to ADFT. It allows the calculation of the first-order response density matrix elements with respect to QM-atom displacements by the following expression:

In eqn (1.122) and denote the first-order response Kohn–Sham and overlap matrix elements in the MO representation, which are defined by:

and

In conventional Kohn–Sham DFT, due to the fact that the perturbed Kohn–Sham matrix depends on the perturbed density matrix, eqn (1.122) has to be solved iteratively to achieve self-consistency.¹²⁰  This iterative evaluation can become a serious computational bottleneck if convergence is slow. In contrast, in the ADFT framework, the corresponding Kohn–Sham matrix only depends on the fitted density.¹²⁵  As a result, McWeeny's SCP method can be reformulated as a non-iterative approach in terms of the response of the fitting coefficients. The resulting ADPT equations yield coupled-perturbed Kohn–Sham (CPKS) analog equation systems, of significantly reduced dimensions that read:^29,125,126 

Here G denotes the Coulomb matrix, eqn (1.12). The Coulomb response matrix, A, and the kernel matrix, F, are given by:

and

The linear equation system of eqn (1.125) keeps the same form irrespective of the nature of the perturbation parameter (electric, magnetic, QM or MM nuclear displacement). The essential difference arises from the perturbation vector elements, , which are characteristic of the specific perturbation. So far, ADPT has been successfully applied for the calculation of static and dynamic polarizabilities,^29,125,126  first hyperpolarizabilities,¹²⁷  Fukui functions^128,129  and second-order nuclear displacements.¹¹⁷  For the latter, the perturbation vector becomes much more complicated due to the perturbation dependency of basis and auxiliary functions. In fact, in this case the calculation of these vector elements can become a computational bottleneck because they are needed for each individual perturbation.

To solve the ADPT equation system, eqn (1.125), we can calculate the inverse of the response matrix, and then solve the equation system analytically. This is possible because the dimension of the ADPT response matrix (N_aux × N_aux) is much smaller than the (N_occ × N_uno) × (N_occ × N_uno) CPKS dimension.¹²⁵  An important consideration is that R is perturbation-independent. This is clearly an advantage, because regardless of the number of perturbations (or perturbation types), R remains unaltered. Thus, a single inversion is needed, and R⁻¹ can be stored on disk and read when required. However, as the system size grows, the response matrix R tends to become ill-conditioned and therefore, prone to numerical instabilities during the matrix inversion. Therefore, a robust and efficient algorithm for solving large, indefinite non-symmetric systems of linear equations is needed.¹³⁰  In deMon2k we use the Eirola–Nevanlinna (EN) algorithm.^27,29,131  As the evaluation of analytic QM second energy derivatives requires the solution of three times the number of QM atoms in the ADPT equation systems, pre-calculation of the Coulomb response matrix, A, and the kernel matrix, F, is convenient. However, the calculation of the Coulomb response matrix introduces an overall (N$aux4$) computational complexity.¹²⁵  Despite its unfavorable scaling, the explicit evaluation of A and F is usually computationally advantageous for second-order energy derivatives because these matrices have to be calculated only once. Analogous to the conventional SCF approach, this implementation for the solution of the ADPT equation systems is called the conventional EN (Con-EN) method.¹¹⁷  A drawback of this implementation, besides its high order formal scaling, is its random-access memory (RAM) demand. Thus, it is only suitable for small systems with up to a few hundred atoms. To overcome this limitation, a direct implementation of the EN algorithm (Dir-EN) is available in deMon2k, too. The advantage of this method is that the explicit computation of the Coulomb response matrix and kernel matrix is avoided. Instead, the action of these matrices on a trial vector p is calculated in the iterative Dir-EN procedure. This reduces the overall (N$aux4$) complexity to (N$aux3$)²⁹  by breaking down the matrix vector multiplication into a series of sub-steps. This allows efficient ERI screening and double asymptotic expansions in the action calculations with the Coulomb response matrix. Similarly, numerical integration screening techniques can be used for the action calculation with the kernel matrix. Besides the reduction in computational complexity, the Dir-EN implementation also has a reduced RAM demand. It requires about 2 N$bas2$ matrices and a few N_aux vectors, compared to the N_bas × N_occ × N_uno RAM demand for the explicit Coulomb response matrix calculation in the Con-EN implementation. This permits ADPT calculations on systems with a thousand or more atoms.²⁹  As for the direct SCF algorithm implemented in deMon2k, the Dir-EN implementation requires the repetitive calculation of three-center ERIs in each iteration step of the EN algorithm. For many perturbations, as is needed for the calculation of analytic second-order energy-derivatives, the repetitive ERI calculation is usually the most demanding computational task in the Dir-EN method. Fortunately, three-center ERI calculation is well developed in deMon2k^30,46,119  and also very well parallelized.^132,133 

Once the ADPT equation systems are solved and the perturbed Coulomb fitting coefficients are obtained, the calculation of the remaining response matrix and vector elements is straightforward. First, the perturbed exchange–correlation fitting coefficients, , are obtained according to:¹¹⁷ 

After that, the perturbed Kohn–Sham matrix elements, which are obtained by direct differentiation of eqn (1.21) with respect to a QM nuclear coordinate, are calculated with the perturbed Coulomb and exchange–correlation fitting coefficients:

Afterward, the perturbed Kohn–Sham matrix elements are transformed into the MO representation using eqn (1.123), and the perturbed density matrix is calculated according to eqn (1.122). Once the perturbed density and Kohn–Sham matrices are available, the perturbed energy weighted density matrix W^(B_η) is calculated as:

With all perturbed matrix and vector elements at hand, the corresponding contributions of the perturbation-dependent second derivatives, eqn (1.121), can be calculated. This concludes the calculation of analytic second-order ADFT energy derivatives with respect to nuclear displacements. For the simulation of infrared (IR) and Raman spectra within the double harmonic approximation,^134,135  second- and third-order mixed energy derivatives with respect to the QM atomic coordinates and components of an external electric field are needed, too. They are implemented as analytical and semi-numerical derivatives in deMon2k.¹³⁶  Finally, we note that the analytic second-order ADFT energy derivative calculation is currently extended to effective and model core potentials employing half-numerical integral evaluation.⁵⁰ 

Although ADPT within ADFT can provide first principles response properties for systems of about a thousand atoms,¹¹⁷  a frequency analysis of a protein with many thousand atoms in water is still computationally unfeasible with a first principles QM-only methodology. For such systems with many thousands of atoms, the QM/MM methodology in deMon2k can be used.^42,44  Analogous to the QM/MM energy expression, eqn (1.38), second-order QM/MM energy derivatives can be expressed as:

The last three terms of the right-hand side of eqn (1.131) are analytic derivatives of the corresponding gradient terms discussed in Section 1.2. Therefore, they are straightforward to calculate and will not be further discussed here. The complete Hessian matrix, schematically depicted in Figure 1.11, can be separated into four regions, labeled as QM–QM, QM–MM, MM–QM and MM–MM in Figure 1.11. By symmetry, the QM–MM and MM–QM blocks are identical. In the QM–QM block, second ADFT energy derivatives with respect to QM atomic coordinates are collected. They are calculated as outlined above. Note, however, that they include, according to eqn (1.119), second derivatives of the QM/MM core Hamilton matrix elements given by eqn (1.36). The elements of the QM–QM block also contain analytic second derivatives of E^QMMM and E^NN with respect to QM atomic coordinates, which are straightforward to calculate. The QM–MM block elements contain contributions from the same three terms, namely E^QM, E^QMMM and E^NN. Again, we discuss only the terms arising from the (mixed) second derivatives of E^QM. Their calculation is simpler because the basis and auxiliary functions of the QM atoms are independent of the coordinates of the MM atoms, and the MM atoms themselves carry none of these functions. Taking this into account, the evaluation of eqn (1.118) for a QM atomic coordinate A_λ and a MM atomic coordinate D_η yields:

Note that the perturbed density matrix elements in eqn (1.132) are the same as those needed for the elements in the QM–QM block. The perturbed core–matrix elements in eqn (1.132) are given by (without employing asymptotic expansions)

and

Finally, the MM–MM block elements contain contributions from all four analytic second derivatives on the right-hand side of eqn (1.131). Again, the E^MM, E^QMMM and E^NN second derivatives with respect to MM atomic coordinates are straightforward to calculate. The E^QM analytic second derivatives with respect to MM atomic coordinates D_λ and D′_η are as follows:

with

Eqn (1.135) implies the calculation of perturbed density matrix elements, , for each corresponding MM degree of freedom. This is expected to be the principal bottleneck for the analytical Hessian calculation of a QM/MM system because there are usually many more MM atoms than QM atoms. Therefore, the ADPT equations must be solved for all MM degrees of freedom, too. Fortunately, the response equation system will have important simplifications as basis and auxiliary functions are independent of MM atom coordinates. These simplifications have been explored already in the framework of conventional DFT/MM calculations using CPKS response equations.^137,138  Despite these improvements, such calculations are still prohibitive for large systems, including enzymes. Thus, a derivation of second-order ADFT energy derivatives within the QM/MM framework is expected to be a useful tool to evaluate frequencies and infrared intensities of complex systems with many thousands of atoms. Work in this direction is currently under development within the deMon developers’ community.

The calculation of magnetic properties such as magnetic shielding and magnetizabilities of very large systems containing thousands of atoms, like those including explicit solvent molecules, can be efficiently handled by QM/MM methods. In the particular case of ADFT QM/MM magnetic property calculations,^139–143  QM regions with around 1000 atoms and 15 000 thousand basis functions can be used. The ADFT magnetic shielding tensor elements, σ_λη, are defined as the second-order ADFT energy derivative with respect to a component of the magnetic field, ℋ_λ, and a component of the nuclear magnetic moment, μ_Cη:

In a QM/MM calculation, only the QM atoms possess nuclear magnetic moments. Therefore, eqn (1.137) is only defined for the QM atoms. Furthermore, the other terms of the QM/MM energy expression, eqn (1.38), do not contribute to this mixed second energy derivative. Thus, expanding eqn (1.137) yields for the shielding tensor elements:

In this equation a and b index gauge-including atomic orbitals (GIAOs)^139,144  centered at the QM nuclei A and B. The GIAOs, ϕ_a(r), are obtained by multiplying (contracted) GTOs with the so-called London factor:

Here a(r) is the field-independent atomic orbital in the form of a (contracted) GTO, and A and G are the position vectors of the QM nucleus A and the gauge, respectively. The use of GIAOs effectively resolves the gauge origin problem in magnetic shielding calculations. The perturbed core matrix elements in eqn (1.138) remain unaltered in QM/MM shielding calculations compared to that of their QM-only counterparts. Their efficient calculation is described in ref. 139. Thus, the essential difference between QM/MM and QM shielding calculations arises from the calculation of the perturbed density matrix elements with respect to the external magnetic field components ℋ_λ. In ADFT GIAO, using LDA or GGA,¹⁴⁵  they are given according to McWeeny's SCP theory^120–124  by:

Here, o and u denote occupied and unoccupied MOs, respectively. The minus sign in eqn (1.140) compared to eqn (1.122) arises from the imaginary perturbation dependency of the GIAO basis according to eqn (1.139). The Kohn–Sham and overlap matrix elements in the MO representations, $ou(λ)$ and $ou(λ)$, are calculated according to eqn (1.123) and (1.124). The underlying perturbed matrix elements are given by:¹⁴⁵ 

and

Because in QM/MM calculations the one electron core Hamilton operator, Ĥ, is augmented by the electrostatic embedding potential from the MM environment, eqn (1.34), we can write eqn (1.142) as:

Thus, the last term in eqn (1.143) must be additionally computed for λ = x, y, z in QM/MM magnetic shielding calculations. Efficient recurrence relations are available in deMon2k for these modified nuclear attraction integrals. Therefore, the computational efficiency of the ADFT-GIAO methodology is preserved in QM/MM calculations.

Another example for QM linear response in the framework of ADFT QM/MM is the calculation of excited states. To this end, we compute excited-state properties through TD-ADFT.¹⁴⁶  In this formalism, the QM/MM excited-state energy, E*, is calculated as a single-electron excitation in the QM region from the QM/MM ground-state, eqn (1.38):

For a given number of excited states, these excitation energies, ω, are obtained by solving an equation system that resembles that of the random phase approximation (RPA):¹⁴⁷ 

Both matrices, A and B, contain the Coulomb and kernel response of the Kohn–Sham matrix, to a change in the electronic density. These matrices are free of any explicit contribution of the MM environment. For pure functionals, they are given by

and

As can be seen from eqn (1.146) the A matrix has, in the diagonal, the energy differences between unoccupied, ϵ_a, and occupied, ϵ_i, canonical MO energies. This term, which is also known as single particle transition energy, is a few orders of magnitude bigger than the other contributions and therefore dominant in the calculation of the excitation energies. In this diagonal term the use of MM embedding contributes implicitly by shifting these MO energies, resulting in calculated excitation energies that are in good agreement with experiments, even in polar solutions.^148,149 

The equation system in eqn (1.145) can be reduced to the following eigenvalue problem

with

and

The calculation of (A − B)^−1/2 is trivial if no Fock exchange is requested, which is the case of the TD-ADFT implementation in deMon2k.^150,151  Another simplification of the TD-ADFT equation system is the so-called Tamm–Dancoff approximation (TDA)¹⁵²  which is obtained by setting B = 0 as:

Although it is an approximation, it has been shown that using the TDA yields the same, or in some cases, even better accuracy than the full equation system.¹⁵² 

Regardless of the type of TD-ADFT, namely eqn (1.148) or (1.151), the TD-ADFT matrix, A or Ω, is of dimension (occ × uno)², with occ and uno being the number of occupied and unoccupied MOs, respectively. The diagonalization, or just storage, of this matrix can become computationally demanding even for very small systems (∼ 20 atoms). Therefore, to avoid storage and explicit diagonalization of the full TD-ADFT matrix, the Davidson method,^153,154  which is a direct diagonalization procedure, is implemented in deMon2k. With this method, it is possible to solve the TD-ADFT equation system for a given number of lowest excitation energies. As an example, the calculated TD-ADFT visible spectrum of Reichardt's dye at the PBE/DZVP-GGA/GEN-A2* level is shown in Figure 1.12. For these calculations the molecular structures were optimized in vacuum, water, methanol and ethanol. All solvent molecules were calculated at the MM level, accounting for a total of around 1500 atoms. The OPLS-AA force-field was used in all QM/MM calculations. As Figure 1.12 shows, the experimentally observed solvent shifts in the excitation spectra are reproduced qualitatively and correctly by our QM/MM calculations. Besides the linear TD-ADFT response just discussed, deMon2k also implements the explicit numerical propagation of TD-ADFT equations in real time (RT-TD-ADFT). This implementation is the subject of Chapter 4.

A major practical obstacle to running QM/MM simulations for complex molecular systems, such as those encountered in biology or in materials science, is to prepare input files. For most systems, deMon2k requires the specification of, for each MM atom, its atom type, as defined in the FFDS file, and importantly its connectivity, i.e. the list of other MM atoms to which it is bonded. Setting up a QM/MM simulation by hand is thus a tedious task that turns into an impossible one as soon as the system of interest contains more than a few tens of atoms. Actually, in practice, one has generally prepared and investigated the system at the MM level before switching to QM/MM, to investigate, for example, a chemical reaction energy profile or to evaluate electronic properties at a QM level of theory. Highly optimized packages for classical molecular simulations like NAMD, LAMMPS and GROMACS to name a few are readily available to generate the topology of the systems and to set up classical MM simulations from geometries in standard formats (e.g. the Protein Data Bank “PDB” format). A powerful strategy is therefore to translate the information contained in the topology files provided by one's favorite MM package to a deMon2k-compatible input file. To this end, we have devised a dedicated program, called QIB (for QM/MM input builder). QIB is currently capable of converting topology files from the AMBER package, CHARMM and GROMACS to the deMon2k format. QIB can also create a QM/MM input file from scratch, given a reference PDB file. Another essential capability of QIB is to take care of a series of technical but important details appearing in QM/MM simulations.

To illustrate the capabilities of QIB, we take the example of a QM/MM simulation of a DNA–protein complex shown in Figure 1.13. The latter has been solvated in a box of 142 Å edge of TIP3P water molecules and contains 258 550 atoms. QIB extracts important information from the AMBER Prmtop file such as the atom names, the residue pointers and the connectivity. In a second step, it identifies for each atom the corresponding atom type within the deMon2k FFDS file. Together with a coordinate file that can be provided under various formats, this is all the information needed to prepare the geometry block in deMon2k. As deMon2k does not currently implement periodic boundary conditions but can account for remote electrostatic effects with the Onsager polarizable model (see Section 4), it is advisable to prune the system to fit within spheres, typically with a 30 to 40 Å radius. Using the information on the residues to which each atom belongs, the pruning can be carried out, keeping the molecular integrity of each residue. In the present case, this leads to a molecular system encompassing around 15 600 atoms. For QM/MM simulations, the user can provide an atom selection that defines a primary QM region. Eventually, the latter is enlarged by a distance criterion, as illustrated in Figure 1.13. This option is useful to include solvation water molecules close to the primary QM region, or to simply generate a QM/MM set-up without the need to identify all atom indexes defining the QM region. A geometry file for the QM region is produced for checking. Note that QIB automatically handles the case of covalent bonds appearing across the QM/MM border. It further suggests a list of atom types for QM atoms in order to calculate Lennard-Jones interactions with MM atoms given by eqn (1.46). Other options are available in QIB, and examples of QM/MM simulations with the molecular systems depicted in Figure 1.13 are reported in Chapter 4 of this book.

In this chapter, we gave an overview of the current status of the QM/MM implementation in deMon2k and the corresponding input support program QIB. As this chapter shows, the QM/MM implementation in deMon2k is a community effort and as such has various perspectives. Although most basic computational tasks such as QM/MM structure optimizations, MD simulations and property calculations are now readily available within the QM/MM module of deMon2k, several larger implementation efforts are still under way. For instance, Chapter 4 introduces a polarizable QM/MMpol upgrade that allows electrostatic induction to be incorporated. We gave a snapshot of this situation in this chapter. The ultimate goal of the efforts described here is to provide a state-of-the-art QM/MM approach that is easy to use and versatile on the one hand and, on the other, computationally powerful enough to address challenging problems in chemistry, materials science and biology. For the latter, a major focus is on the QM methodology that is usually the computational bottleneck. Because our developments are in the framework of deMon2k, they are generally accessible and free of charge, to the academic community.

ADFT

Auxiliary density functional theory

ADPT

Auxiliary density perturbation theory

CCM

Cyclic cluster model

CMO

Canonical molecular orbital

CPKS

Coupled-perturbed Kohn–Sham

deMon2k

density of Montréal 2000

DF-DFT

Density fitted density functional theory

DFA

Density functional approximation

DFT

Density functional theory

ECP

Effective core potential

EN

Eirola–Nevanlinna

ERI

Electron repulsion integrals

FCM

Free cluster model

FFDS

Force-field dataset

GGA

Generalized gradient approximation

GIAO

Gauge-including atomic orbital

GTO

Gaussian-type orbital

IR

Infrared

LCGTO

Linear combination of Gaussian-type orbitals

LDA

Local density approximation

LMO

Localized molecular orbital

MCP

Model core potential

MD

Molecular dynamics

MM

Molecular mechanics

MO

Molecular orbital

MTK

Martyna–Tobias–Klein

NHC

Nosé–Hoover chain

NPT

Isothermal–isobaric

NVE

Microcanonical

NVT

Canonical

PCM

Polarizable continuum model

PES

Potential energy surface

QIB

QM/MM input builder

QM

Quantum mechanics

QM/MM

Quantum mechanical/molecular mechanical

RAM

Random-access memory

RI

Resolution-of-the-identity

RPA

Random phase approximation

RT-TD-ADFT

Real-time time-dependent auxiliary density functional theory

SCF

Self-consistent field

SCP

Self-consistent perturbation

TDA

Tamm–Dancoff approximation

TD-ADFT

Time-dependent auxiliary density functional theory

TRM

Trust region method

TS

Transition state

WSC

Wigner–Seitz cell

JDSR (490270), LIHS (455854) and LLS (397526) gratefully acknowledge CONACyT PhD fellowships. The work reported here has been supported by the SENER–CONACyT project 280158 and the CONACyT projects CB-258647 (BAZG) and A1-S-11929 (AMK). For the QM/MM deMon2k program development, computational resources from the CONACyT infrastructure project GIC-268251 and the SEP–Cinvestav project No. 65 are used.