Multiscale Dynamics Simulations: Nano and Nanobio Systems in Complex Environments
 1.1 Introduction
 1.2 Energies and Gradients
 1.3 Local Structure Minimization and Transitionstate Search
 1.4 Molecular Dynamics
 1.5 Secondorder Energy Derivatives
 1.6 QM/MM Magnetic Shielding and Excitedstate Calculations
 1.7 Transformation Program for deMon2k Input Generation
 1.8 Summary
 Abbreviations
Chapter 1: QM/MM with Auxiliary DFT in deMon2k^{1}

Published:24 Sep 2021

Special Collection: 2021 ebook collection
Juan D. SamaniegoRojas, LuisIgnacio HernándezSegura, Luis LópezSosa, Rogelio I. DelgadoVenegas, Badhin Gomez, JeanChristophe Lambry, Aurelien de la Lande, Tzonka Mineva, José Alejandre, Bernardo A. ZúñigaGutiérrez, Roberto FloresMoreno, Patrizia Calaminici, Gerald Geudtner, Andreas M. Köster, 2021. "QM/MM with Auxiliary DFT in deMon2k^{1}", Multiscale Dynamics Simulations: Nano and Nanobio Systems in Complex Environments, Dennis R. Salahub, Dongqing Wei
Download citation file:
This chapter describes the theoretical background of the quantum mechanical/molecular mechanical (QM/MM) implementation in deMon2k within the framework of auxiliary density functional theory (ADFT). It aims to give the reader an overview of the current state of the art of this QM/MM implementation and perspectives for its future development. To this end, we first derive the ADFT working equations for the QM and QM/MM energy and gradient expressions. Based on the joint QM/MM gradient expression, we present algorithms for QM/MM structure optimizations, transitionstate searches and molecular dynamics simulations. The use of auxiliary density perturbation theory (ADPT) in the framework of QM/MM is discussed using illustrative implementations including analytic secondorder ADFT energy derivatives, nuclear magnetic resonance chemical shift calculations and excited state calculations using timedependent ADFT. The chapter closes with the description of a transformation program used to generate deMon2k QM/MM inputs.
1.1 Introduction
Over the last two decades, many quantum mechanical/molecular mechanical (QM/MM) implementations have emerged based on the pioneering works of Warshel and Karplus^{1 } and Warshel and Levitt^{2 } 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 forcefield 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 wellknown Hartree–Fock method, it incorporates electron correlation through the exchange–correlation functional. In the linear combination of Gaussiantype 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 twoelectron Coulomb repulsion into fourcenter electron repulsion integrals (ERIs). Already in very early molecular LCGTODFT Kohn–Sham implementations, the calculation of the fourcenter ERIs was identified as a critical computational bottleneck. To overcome this bottleneck, Dunlap and coworkers^{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 fourcenter ERIs. Note, however, that the variational bound shifts below the corresponding fourcenter ERI Kohn–Sham bound. For relative energies, this shift in the variational bound has no effect. Thus, Kohn–Sham LCGTODFT approaches incorporating the variational fitting of the Coulomb potential yield the same relative energies as their fourcenter ERI counterparts at a considerably reduced computational cost. The same holds for socalled resolutionoftheidentity (RI) approaches^{7 } because RI is a direct result of the variational fitting of the Coulomb potential and, therefore, possesses the same variational bound.^{8 } To avoid fourcenter ERI calculations for hybrid functionals, the variational fitting of the Fock potentials^{9–15 } was introduced into deMon2k,^{16 } too. The resulting fourcenter ERIfree 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 (DFDFT).
The computational bottleneck in mediumsized (50 to 500 atoms) DFDFT 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 twoelectron Coulomb repulsion, it has received much less attention. Early LCGTODFT implementations like deMonKS^{17 } or DGAUSS^{18 } used the seminumerical fit of the exchange–correlation potential with auxiliary functions as proposed by Sambe and Felton^{19 } to circumvent this bottleneck. This approach gives reasonable energies but is prone to numerical instabilities in energy derivatives.^{20 } 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)^{25 } 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^{26 } and Eirola–Nevanlinna,^{27 } for the density fitting^{28 } and perturbation theory^{29 } 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 transitionstate searches based on the trustregion method for local optimization. In the next section, the propagator equations for QM/MM molecular dynamics (MD) simulations in various closedsystem ensembles using factorized Liouville operators are presented. Section 1.5 introduces auxiliary density perturbation theory (ADPT) for the calculation of analytic secondorder 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 excitedstate 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 stateoftheart ADFT QM/MM implementation in deMon2k is given in the last section of this chapter.
1.2 Energies and Gradients
The QM LCGTO Kohn–Sham energy expression is given by:
In our notation Greek letters, here μ, ν, σ and τ denote (contracted) atomcentered Gaussiantype orbitals. For clarity of presentation, we restrict ourselves to closedshell systems. The closedshell 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 monoelectronic 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 fourcenter twoelectron integrals responsible for the N$bas4$ scaling of the method:
In this shorthand ERI notation ‖ represents the twoelectron Coulomb operator 1/r_{1} − r_{2}. 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 threecenter 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 twoelectron Coulomb repulsion in eqn (1.1) introduces a formal N$bas4$ scaling due to the fourcenter ERIs. To circumvent this computational bottleneck, we employ the variational fitting of the Coulomb potential based on the following secondorder 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, atomcentered primitive Hermite Gaussian functions, indicated by a bar, are used as auxiliary functions.^{30 } 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 semipositive definite.^{31 } Thus, we obtain from eqn (1.8):
Inserting this inequality into the Kohn–Sham energy expression, eqn (1.1), yields the QM DFDFT energy expression:
The Coulomb fitting coefficients, x_{k̄}, in eqn (1.10) are calculated by minimizing the secondorder error functional of eqn (1.6), which is accompanied by an energy maximization:^{32 }
Collecting these equations for all auxiliary functions yields the inhomogeneous linear equation system,
which must be solved in each selfconsistent field (SCF) iteration. In eqn (1.12)G and J are the Coulomb matrix and vector, respectively, given by
As for the QM Kohn–Sham energy, the QM DFDFT energy is variational, due to the variational fitting, but with a lower energy bound. For relative energies this shift of the variational bound has no effect.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:^{22 }
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),
as: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 DFDFT 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^{−1}) 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.^{28 } 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 GENA2*^{33 } 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 DFDFT) 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 rangeseparated hybrid functionals in ADFT, we refer the interested reader to ref. 15 and 34–36. The computational drawback arises from the Fock exchange energy,
which contains fourcenter ERIs and, therefore, will jeopardize the computational benefits of ADFT. To overcome this problem, we implemented the variational fitting of the Fock potential^{13,15 } in deMon2k. To this end, we express the Fock exchange energy in terms of occupied MOs,asFor 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 secondorder error functional:^{13,15 }
Expanding the transition densities in eqn (1.26) yields:
Because ε$2F$ is (termwise) seminegative definite^{35 } the following inequality holds: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 fourcenter ERIfree total Hartree–Fock energies for molecules with thousands of basis functions are nearly indistinguishable (< 1 kcal mol^{−1}) from their fourcenter ERI counterparts! Corresponding binding energy differences between four and threecenter Hartree–Fock calculations are below 10^{−2} kcal mol^{−1}. From eqn (1.28) the following linear equation systems,
with
are obtained. Therefore, the Fock exchangefitting coefficients are given by:
Inserting these Fock exchangefitting 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 fourcenter 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^{13 } 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 twostep localization is used. First, an incomplete Cholesky decomposition^{37 } 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
After ordering these atomic populations from the largest to the smallest, we sum them up until a threshold value (default is 0.9995 for final energy calculation) is reached. All atoms that contribute to this sum define the atomic domain for the given LMO. The auxiliary functions on these atoms define the local auxiliary function set given by the blue region in Figure 1.1c. The corresponding local basis set is defined by all basis functions of the domain atoms and augmented by basis functions from neighboring atoms with significant overlap (by default ≥10^{−6}) into the domain. This is the red region in Figure 1.1d. Because of this localization, each LMO has a particular Coulomb matrix. The computational cost for computing all these local G and G ^{− 1} matrices^{40 } for all occupied LMOs is, in larger systems, more than overcompensated by the reduced dimensionality of these matrices. Thus, the final QM energy expression used for ADFT QM/MM calculations has the form: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$\mu \nu 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:^{41 }
This general definition allows inclusion of higher point moments on MM atoms, e.g. for polarizable forcefields.^{42 } 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^{43 } 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 farfield integrals.^{44 } The farfield QM/MM electrostatic interaction integrals are obtained by the following asymptotic expansion:^{45 }
The superscript A indicates the center of the asymptotic expansion. There are always two expansions, one for the center where atomic orbital α is located and the other for the center where the atomic orbital β is located. The extension of this asymptotic expansion to higher moments is possible as shown in the double asymptotic expansion of ERIs.^{46 } Thus, QM/MM core Hamilton matrix elements in deMon2k are calculated as (only the asymptotic expansion with respect to center A is explicitly shown):with
In eqn (1.36) the sum limits Near and Far refer to MM atoms in the near and farfield 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 farfield 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.^{44 } The QM/MM modified core matrix elements enter the corresponding ADFT Kohn–Sham matrix according to eqn (1.21) where H$\mu \nu 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 DFDFT and ADFT. The DFAs available for QM/MM calculations in deMon2k are LDA, GGA and metaGGA (only within DFDFT) as well as global and rangeseparated hybrids. By default, the LDA is used in the form of Dirac exchange^{47 } and VWN correlation.^{48 } Because, in ADFT, LDA and GGA performance is similar, we recommend GGA functionals for QM/MM calculations. In particular, PBE^{49 } in combination with the GGA optimized DZVPGGA basis and GENA2* auxiliary function set^{33 } 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).^{50 } The ECP approach can also be used for the introduction of capping potentials^{51 } in QM/MM calculations. We also recommend the PBE functional for QM/MM calculations because analytic ADFT kernel formulas^{52 } 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.^{53 } If van der Waals interactions are important, e.g. in QM/MM structure optimizations, the DFDFT 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 lessthanideal parallel scaling.^{40 } Instead, we advocate composite approaches^{57 } consisting of GGA structure optimization and singlepoint 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 Dunningtype basis sets.^{60 } The handling of the threecenter 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 threecenter ERIs are repeatedly calculated (twice) in each SCF iteration. If sufficient random access memory (RAM) is available, the use of the socalled 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^{−5} 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, nonpolarizable, forcefields according to which the MM energy expression reads:
The first five terms correspond to bonded contributions, whereas the last two are nonbonded terms. The extension to a polarizable QM/MMpol methodology that further includes electrostatic induction within the framework of the pointcharge 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_{0} and the force constant k_{b} are forcefield 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_{0} and force constant k_{l} are given by the forcefield, too.
The angle bending in eqn (1.42) is modeled by a harmonic potential as well. The corresponding force constant, k_{ϑ}, and equilibrium angle, ϑ_{0}, are forcefield 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^{2} 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 forcefield 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^{65 } (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 nonbonded interactions, the energy expressions used are the wellknown van der Waals (aka LennardJones) 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 forcefield. 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:
By default in deMon2k, an arithmetic mean is used to obtain σ_{AB}, whereas the geometric mean is used for the calculation of ϵ_{AB}. All parameters are read from the FFDS file, which contains the OPLSAA forcefield^{66 } by default. Alternatively, the user can select the family of AMBER forcefield parameters,^{67 } which are available in a common format with corresponding FFDS files, too.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 LennardJones 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 abovedescribed combination rules from atomic parameters. To this end, MM atomic ϵ and σ parameters are assigned to the QM atoms by the QM/MM keyword^{65 } in deMon2k. The last term in eqn (1.38) is the electrostatic repulsion between the QM pointcharge 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 higherorder SCF energy derivatives calculated in the TDADFT 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^{68 } 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 threecenter ERI Fock exchange energy. The closedshell energyweighted densitymatrix 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^{31 } analogous to the corematrix elements are available in deMon2k, too. The DFDFT 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 nonbonded 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 nonbonded 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 pointcharge 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.
1.3 Local Structure Minimization and Transitionstate Search
The potential energy surface (PES) of a molecular system is a multidimensional function arising from the Born–Oppenheimer approximation^{71 } 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, quasiNewton 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 quasiNewton 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 secondorder 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_{0}. 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 constrainedminimization problem is obtained:
For the solution of eqn (1.72) the following Lagrange function is introduced:^{85 }
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,^{88 } Fletcher,^{89 } Goldfarb^{90 } and Shanno,^{91 } 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 secondorder energy derivatives such as steepest descent (see keyword STEPTYPE^{65 }) 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 (firstorder) 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 socalled 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 socalled 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^{94 } is currently being extended for QM/MM calculations in one of our laboratories. In transition state searches the unbiased PowellsymmetricBroyden^{95 } update formula is used in deMon2k by default:
1.4 Molecular Dynamics
Whereas the QM/MM molecular structure optimizations in deMon2k utilize algorithms from QM calculations, the QM/MM MD simulations use MMinspired 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^{96 } 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 phasespace 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^{97 } of the classical propagator:
Here, the operator of eqn (1.84) has been arbitrarily separated into two terms, iL_{1} and iL_{2}, 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 phasespace vector of a system results in the wellknown velocity Verlet algorithm.^{98 } 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^{99 } are available in deMon2k. All of them work by controlling the kinetic energy fluctuations of the system by changing the particle velocities. The extended phasespace 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 selfconsistent 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 twodimensional 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 LennardJones 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^{−1} 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, longrange Coulomb interactions must be added. To this end, the lattice sum algorithm of Ewald in its classical^{111 } 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 selfinteraction. 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.^{112 } 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 pointcharge potential becomes shortrange 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 pointcharge 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 selfinteraction.
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 ,^{114 } 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 longrange Coulomb lattice sums according to eqn (1.113) against deMon2k embedding and periodic Tinker^{115 } 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^{−1} 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,^{116 } 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 selfinteraction term, which only appears in eqn (1.113) due to the pointlike 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.
1.5 Secondorder Energy Derivatives
The calculation of secondorder 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^{2}, with M being the number of atoms in the system. Although the firstorder electronic energy derivatives are relatively easy to calculate, secondorder derivatives are less straightforward. The simplest way to calculate secondorder derivatives is to follow a finite difference method where typically two firstorder 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 secondorder 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 secondorder derivative. For an ADFT QMonly 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 firstorder response density matrix as well as the response Coulomb and exchange–correlation fitting coefficients. Based on this, the expression of the secondorder ADFT energy derivatives can be partitioned into two contributions as:^{117 }
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^{52 } or by finite differences if only analytic exchange–correlation potential expressions are available.^{53 } The density, P_{μν}, and energyweighted 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 singlepoint energy calculation. The second contribution to the QM secondorder ADFT energy derivative in eqn (1.118) contains only terms that depend on firstorder 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 selfconsistent perturbation (SCP) theory^{120–124 } to ADFT. It allows the calculation of the firstorder response density matrix elements with respect to QMatom displacements by the following expression:
In eqn (1.122) and denote the firstorder 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 selfconsistency.^{120 } 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.^{125 } As a result, McWeeny's SCP method can be reformulated as a noniterative approach in terms of the response of the fitting coefficients. The resulting ADPT equations yield coupledperturbed 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,^{127 } Fukui functions^{128,129 } and secondorder nuclear displacements.^{117 } 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.^{125 } An important consideration is that R is perturbationindependent. 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^{−1} can be stored on disk and read when required. However, as the system size grows, the response matrix R tends to become illconditioned and therefore, prone to numerical instabilities during the matrix inversion. Therefore, a robust and efficient algorithm for solving large, indefinite nonsymmetric systems of linear equations is needed.^{130 } 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, precalculation 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.^{125 } Despite its unfavorable scaling, the explicit evaluation of A and F is usually computationally advantageous for secondorder 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 (ConEN) method.^{117 } A drawback of this implementation, besides its high order formal scaling, is its randomaccess 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 (DirEN) 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 DirEN procedure. This reduces the overall (N$aux4$) complexity to (N$aux3$)^{29 } by breaking down the matrix vector multiplication into a series of substeps. 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 DirEN 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 ConEN implementation. This permits ADPT calculations on systems with a thousand or more atoms.^{29 } As for the direct SCF algorithm implemented in deMon2k, the DirEN implementation requires the repetitive calculation of threecenter ERIs in each iteration step of the EN algorithm. For many perturbations, as is needed for the calculation of analytic secondorder energyderivatives, the repetitive ERI calculation is usually the most demanding computational task in the DirEN method. Fortunately, threecenter 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:^{117 }
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 perturbationdependent second derivatives, eqn (1.121), can be calculated. This concludes the calculation of analytic secondorder 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 thirdorder 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 seminumerical derivatives in deMon2k.^{136 } Finally, we note that the analytic secondorder ADFT energy derivative calculation is currently extended to effective and model core potentials employing halfnumerical integral evaluation.^{50 }
Although ADPT within ADFT can provide first principles response properties for systems of about a thousand atoms,^{117 } a frequency analysis of a protein with many thousand atoms in water is still computationally unfeasible with a first principles QMonly 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), secondorder QM/MM energy derivatives can be expressed as:
The last three terms of the righthand 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 righthand 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 secondorder 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.
1.6 QM/MM Magnetic Shielding and Excitedstate Calculations
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 secondorder 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 gaugeincluding 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 socalled London factor:
Here a(r) is the fieldindependent 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 QMonly 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,^{145 } 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(\lambda )$ and $ou(\lambda )$, are calculated according to eqn (1.123) and (1.124). The underlying perturbed matrix elements are given by:^{145 }
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 ADFTGIAO 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 excitedstate properties through TDADFT.^{146 } In this formalism, the QM/MM excitedstate energy, E*, is calculated as a singleelectron excitation in the QM region from the QM/MM groundstate, 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):^{147 }
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 TDADFT implementation in deMon2k.^{150,151 } Another simplification of the TDADFT equation system is the socalled Tamm–Dancoff approximation (TDA)^{152 } 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.^{152 }
Regardless of the type of TDADFT, namely eqn (1.148) or (1.151), the TDADFT matrix, A or Ω, is of dimension (occ × uno)^{2}, 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 TDADFT 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 TDADFT equation system for a given number of lowest excitation energies. As an example, the calculated TDADFT visible spectrum of Reichardt's dye at the PBE/DZVPGGA/GENA2* 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 OPLSAA forcefield 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 TDADFT response just discussed, deMon2k also implements the explicit numerical propagation of TDADFT equations in real time (RTTDADFT). This implementation is the subject of Chapter 4.
1.7 Transformation Program for deMon2k Input Generation
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 deMon2kcompatible 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 setup 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 LennardJones 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.
1.8 Summary
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 stateoftheart 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.
Abbreviations
 ADFT
Auxiliary density functional theory
 ADPT
Auxiliary density perturbation theory
 CCM
Cyclic cluster model
 CMO
Canonical molecular orbital
 CPKS
Coupledperturbed Kohn–Sham
 deMon2k
density of Montréal 2000
 DFDFT
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
Forcefield dataset
 GGA
Generalized gradient approximation
 GIAO
Gaugeincluding atomic orbital
 GTO
Gaussiantype orbital
 IR
Infrared
 LCGTO
Linear combination of Gaussiantype 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
Randomaccess memory
 RI
Resolutionoftheidentity
 RPA
Random phase approximation
 RTTDADFT
Realtime timedependent auxiliary density functional theory
 SCF
Selfconsistent field
 SCP
Selfconsistent perturbation
 TDA
Tamm–Dancoff approximation
 TDADFT
Timedependent 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 CB258647 (BAZG) and A1S11929 (AMK). For the QM/MM deMon2k program development, computational resources from the CONACyT infrastructure project GIC268251 and the SEP–Cinvestav project No. 65 are used.
Dedicated to the memory of Rogelio Isaac DelgadoVenegas who passed away while this chapter was being completed. He is deeply missed.