Skip to Main Content
Skip Nav Destination

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, transition-state 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 second-order ADFT energy derivatives, nuclear magnetic resonance chemical shift calculations and excited state calculations using time-dependent ADFT. The chapter closes with the description of a transformation program used to generate deMon2k QM/MM inputs.

Over the last two decades, many quantum mechanical/molecular mechanical (QM/MM) implementations have emerged based on the pioneering works of Warshel and Karplus1  and Warshel and Levitt2  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 Nbas4 scaling, with Nbas 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-workers5,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) approaches7  because RI is a direct result of the variational fitting of the Coulomb potential and, therefore, possesses the same variational bound.8  To avoid four-center ERI calculations for hybrid functionals, the variational fitting of the Fock potentials9–15  was introduced into deMon2k,16  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-KS17  or DGAUSS18  used the semi-numerical fit of the exchange–correlation potential with auxiliary functions as proposed by Sambe and Felton19  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 MINRES26  and Eirola–Nevanlinna,27  for the density fitting28  and perturbation theory29  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:

Equation 1.1

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:

Equation 1.2

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:

Equation 1.3

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 Nbas4 scaling of the method:

Equation 1.4

In this short-hand ERI notation ‖ represents the two-electron Coulomb operator 1/|r1r2|. 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:

Equation 1.5

As already mentioned, the two-electron Coulomb repulsion in eqn (1.1) introduces a formal Nbas4 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 

Equation 1.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:

Equation 1.7

In deMon2k, atom-centered 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:

Equation 1.8

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

Equation 1.9

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

Equation 1.10

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

Equation 1.11

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

Equation 1.12

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

graphic
As for the QM Kohn–Sham energy, the QM DF-DFT 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 Exc must be changed to the auxiliary density. Thus, the QM ADFT energy is given by:22 

Equation 1.13

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:

Equation 1.14

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

Equation 1.15

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

Equation 1.16

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

graphic
as:

Equation 1.17

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

Equation 1.18

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−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:

Equation 1.19

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

Equation 1.20

In eqn (1.20)L denotes the exchange–correlation vector with elements L = 〈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 GEN-A2*33  presents no computational challenge. With the Coulomb and exchange–correlation coefficients at hand, the ADFT Kohn–Sham matrix can be expressed as:

Equation 1.21

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:

Equation 1.22

Here α represents the mixing parameter for the Fock exchange energy, ExF, 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,

graphic
which contains four-center ERIs and, therefore, will jeopardize the computational benefits of ADFT. To overcome this problem, we implemented the variational fitting of the Fock potential13,15  in deMon2k. To this end, we express the Fock exchange energy in terms of occupied MOs,
graphic
as

Equation 1.23

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

Equation 1.24
Equation 1.25

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

Equation 1.26

Expanding the transition densities in eqn (1.26) yields:

graphic
Because ε2F is (termwise) semi-negative definite35  the following inequality holds:

Equation 1.27

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

Equation 1.28

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−1) from their four-center ERI counterparts! Corresponding binding energy differences between four- and three-center Hartree–Fock calculations are below 10−2 kcal mol−1. From eqn (1.28) the following linear equation systems,

Equation 1.29

with

Equation 1.30

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

Equation 1.31

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:

Equation 1.32

Direct implementation of eqn (1.32) leads to an algorithm that scales as (Naux × Nbas2 × Nocc) with Naux and Nocc being the number of auxiliary functions and occupied MOs, respectively. Such an algorithm is useful only when (Nbas) ≫ (Nocc).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.

Figure 1.1

Canonical (a) and localized (b) molecular orbitals of n-C10H22. The atomic centers for the local auxiliary function set (c) and local basis set (d) are indicated by blue and red regions, respectively.

Figure 1.1

Canonical (a) and localized (b) molecular orbitals of n-C10H22. The atomic centers for the local auxiliary function set (c) and local basis set (d) are indicated by blue and red regions, respectively.

Close modal

The recently developed variational fitting of the Fock potential13  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 decomposition37  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

graphic
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 matrices40  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:

Equation 1.33

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:

Equation 1.34

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

Equation 1.35

This general definition allows inclusion of higher point moments on MM atoms, e.g. for polarizable force-fields.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. kx = ky = kz = 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 relations43  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.44  The far-field QM/MM electrostatic interaction integrals are obtained by the following asymptotic expansion:45 

graphic
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):

Equation 1.36

with

Equation 1.37

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.44  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 exchange47  and VWN correlation.48  Because, in ADFT, LDA and GGA performance is similar, we recommend GGA functionals for QM/MM calculations. In particular, PBE49  in combination with the GGA optimized DZVP-GGA basis and GEN-A2* auxiliary function set33  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 potentials51  in QM/MM calculations. We also recommend the PBE functional for QM/MM calculations because analytic ADFT kernel formulas52  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 DF-DFT and ADFT LDA or GGA energy expressions must be augmented with empirical dispersion terms54–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.40  Instead, we advocate composite approaches57  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.60  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 approach61,62  is very beneficial for the computational performance. Tightening of the adaptive grid accuracy63,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:

Equation 1.38

In eqn (1.38)EMM, EQMMM and ENN 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:

Equation 1.39

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:

Equation 1.40
Equation 1.41
Equation 1.42
Equation 1.43

In the first equation, the distance is calculated as d = |AB|. The equilibrium distance d0 and the force constant kb 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 l0 and force constant kl are given by the force-field, too.

Figure 1.2

Urey–Bradley interaction (top) as modeled by a spring-like system formed by the B and C atoms. Bottom: definition of a proper torsion (left) over consecutive bond connected atoms and definition of an improper torsion (right) where consecutive bond connectivity is not mandatory.

Figure 1.2

Urey–Bradley interaction (top) as modeled by a spring-like system formed by the B and C atoms. Bottom: definition of a proper torsion (left) over consecutive bond connected atoms and definition of an improper torsion (right) where consecutive bond connectivity is not mandatory.

Close modal

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 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 sp2 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, Vn, 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 file65  (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:

Equation 1.44
Equation 1.45

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:

graphic
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 OPLS-AA force-field66  by default. Alternatively, the user can select the family of AMBER force-field 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, EQMMM, is modeled in deMon2k by a Lennard-Jones potential:

Equation 1.46

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 keyword65  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:

Equation 1.47

The MM atomic point charges, QD, 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:

Equation 1.48

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.

Figure 1.3

Diagram of deMon2k energy modules for QM/MM calculations and their relation to gradient and property calculations (see also Section 1.5 and 1.6).

Figure 1.3

Diagram of deMon2k energy modules for QM/MM calculations and their relation to gradient and property calculations (see also Section 1.5 and 1.6).

Close modal

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

Equation 1.49

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 force68  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:

Equation 1.50

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

Equation 1.51
Equation 1.52

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:

Equation 1.53

We note that for these derivative integrals (electric field integrals), asymptotic expansions31  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:

Equation 1.54
Equation 1.55

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

Equation 1.56
Equation 1.57
Equation 1.58

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

Equation 1.59

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 φ:

Equation 1.60

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:

Equation 1.61
Equation 1.62
Equation 1.63
Equation 1.64

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:

Equation 1.65
Equation 1.66

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

Equation 1.67
Equation 1.68

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:

Equation 1.69
Equation 1.70

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 approximation71  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.

Figure 1.4

Stationary points on a potential energy surface (PES) model.

Figure 1.4

Stationary points on a potential energy surface (PES) model.

Close modal

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 methods72–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:

Equation 1.71

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 E0. 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:

Equation 1.72

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

Equation 1.73

From the stationary condition of this Lagrange function

Equation 1.74

the step direction p is defined as:

Equation 1.75

For the determination of the Lagrange multiplier λ that multiplies the identity matrix in eqn (1.75), standard techniques from the literature73,85–87  are used. Similarly, to update the (approximated) Hessian matrix the BFGS update formula, named after its inventors Broyden,88  Fletcher,89  Goldfarb90  and Shanno,91  is used:

Equation 1.76

In eqn (1.76)γ = gI + 1gI and δ = xI + 1xI 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 STEPTYPE65 ) 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.

Figure 1.5

Transition state (TS) nature on the model PES.

Figure 1.5

Transition state (TS) nature on the model PES.

Close modal

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

Equation 1.77

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

Equation 1.78

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 algorithm94  is currently being extended for QM/MM calculations in one of our laboratories. In transition state searches the unbiased Powell-symmetric-Broyden95  update formula is used in deMon2k by default:

Equation 1.79

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 operator96  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:

Equation 1.80
Equation 1.81

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

Equation 1.82
Equation 1.83

The Liouville operator used in these expressions is given by:

Equation 1.84

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 = (RM, PM), 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:

Equation 1.85

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 factorization97  of the classical propagator:

Equation 1.86

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

Equation 1.87
Equation 1.88

As a result, the propagator for NVE simulations reads:

Equation 1.89

The application of this operator to the phase-space vector of a system results in the well-known 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 thermostats99  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:

Equation 1.90
Equation 1.91
Equation 1.92
Equation 1.93
Equation 1.94

In these equations the k index loops over the N thermostats, each having a unitless position ηk, a mass Qk 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:

Equation 1.95
Equation 1.96

Here, kB 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 iLT accounting for the presence of the thermostat variables ηk and Pηk. This third term has the form:

Equation 1.97

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

Equation 1.98

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:

Equation 1.99

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:

Equation 1.100

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:

Equation 1.101
Equation 1.102
Equation 1.103
Equation 1.104

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:

Equation 1.105
Equation 1.106

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

Equation 1.107

The external pressure, Pext, 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:

Equation 1.108
Equation 1.109

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.

Figure 1.6

Translation scheme and interactions of atom A in a cyclic cluster model (CCM) calculation.

Figure 1.6

Translation scheme and interactions of atom A in a cyclic cluster model (CCM) calculation.

Close modal

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−1 at much smaller cluster sizes, i.e. at smaller average coordination numbers, than FCM clusters.

Figure 1.7

Cohesive energy of solid argon as a function of cluster size calculated using the free cluster model (FCM) and the cyclic cluster model (CCM).

Figure 1.7

Cohesive energy of solid argon as a function of cluster size calculated using the free cluster model (FCM) and the cyclic cluster model (CCM).

Close modal

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 classical111  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:

Equation 1.110

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.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:

Equation 1.111

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:

Equation 1.112

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 α.

Figure 1.8

One-dimensional point and Gaussian charge distributions. The original potential (right-hand side) is obtained by adding the shielded and compensating potentials.

Figure 1.8

One-dimensional point and Gaussian charge distributions. The original potential (right-hand side) is obtained by adding the shielded and compensating potentials.

Close modal

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:

Equation 1.113

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 ,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 long-range Coulomb lattice sums according to eqn (1.113) against deMon2k embedding and periodic Tinker115  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.

Figure 1.9

CCM cluster of four water molecules (a) together with their cubic (b) and spherical (c) embedding supercells of increasing size.

Figure 1.9

CCM cluster of four water molecules (a) together with their cubic (b) and spherical (c) embedding supercells of increasing size.

Close modal
Figure 1.10

Energy convergence of deMon2k embedding calculations toward the Ewald energy results (left) and comparison of the deMon2k CCM force norms against periodic Tinker calculations (right).

Figure 1.10

Energy convergence of deMon2k embedding calculations toward the Ewald energy results (left) and comparison of the deMon2k CCM force norms against periodic Tinker calculations (right).

Close modal

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:

Equation 1.114

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:

Equation 1.115
Equation 1.116

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 9M2, 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:

Equation 1.117

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:117 

Equation 1.118

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

Equation 1.119

Here, fxc[] 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 

Equation 1.120

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 functional52  or by finite differences if only analytic exchange–correlation potential expressions are available.53  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, and exchange–correlation fitting coefficients, z, 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:

Equation 1.121

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) theory120–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:

Equation 1.122

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

Equation 1.123

and

Equation 1.124

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.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 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 

Equation 1.125

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

Equation 1.126

and

Equation 1.127

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 functions128,129  and second-order 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 (Naux × Naux) is much smaller than the (Nocc × Nuno) × (Nocc × Nuno) CPKS dimension.125  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−1 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.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, 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 (Naux4) computational complexity.125  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.117  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 (Naux4) complexity to (Naux3)29  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 Nbas2 matrices and a few Naux vectors, compared to the Nbas × Nocc × Nuno 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.29  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 deMon2k30,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 

Equation 1.128

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:

Equation 1.129

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:

Equation 1.130

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.136  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.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 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:

Equation 1.131

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 EQMMM and ENN with respect to QM atomic coordinates, which are straightforward to calculate. The QM–MM block elements contain contributions from the same three terms, namely EQM, EQMMM and ENN. Again, we discuss only the terms arising from the (mixed) second derivatives of EQM. 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:

Equation 1.132

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)

Equation 1.133

and

Equation 1.134

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 EMM, EQMMM and ENN second derivatives with respect to MM atomic coordinates are straightforward to calculate. The EQM analytic second derivatives with respect to MM atomic coordinates Dλ and Dη are as follows:

Equation 1.135

with

Equation 1.136

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.

Figure 1.11

Hessian matrix structure for a QM/MM system.

Figure 1.11

Hessian matrix structure for a QM/MM system.

Close modal

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η:

Equation 1.137

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:

Equation 1.138

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:

Equation 1.139

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,145  they are given according to McWeeny's SCP theory120–124  by:

Equation 1.140

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:145 

Equation 1.141

and

Equation 1.142

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:

Equation 1.143

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.146  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):

Equation 1.144

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 

Equation 1.145

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

Equation 1.146

and

Equation 1.147

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

Equation 1.148

with

Equation 1.149

and

Equation 1.150

The calculation of (AB)−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)152  which is obtained by setting B = 0 as:

Equation 1.151

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 TD-ADFT, namely eqn (1.148) or (1.151), the TD-ADFT 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 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.

Figure 1.12

Calculated Reichardt's dye spectra in vacuum (black dashed line), water (black solid line), methanol (red line) and ethanol (violet line). Photographs of solutions containing the molecule in water, methanol and ethanol are shown in the inset, top right.

Figure 1.12

Calculated Reichardt's dye spectra in vacuum (black dashed line), water (black solid line), methanol (red line) and ethanol (violet line). Photographs of solutions containing the molecule in water, methanol and ethanol are shown in the inset, top right.

Close modal

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.

Figure 1.13

Pictogram for the automatic generation of deMon2k input files from the standard MM package topology and geometry files by QIB for QM/MM simulations. The program generates a minimal deMon2k input file and an xyz file with the geometry of the QM region encompassing link atoms.

Figure 1.13

Pictogram for the automatic generation of deMon2k input files from the standard MM package topology and geometry files by QIB for QM/MM simulations. The program generates a minimal deMon2k input file and an xyz file with the geometry of the QM region encompassing link atoms.

Close modal

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.

1

Dedicated to the memory of Rogelio Isaac Delgado-Venegas who passed away while this chapter was being completed. He is deeply missed.

1.
Warshel
A.
,
Karplus
M.
,
J. Am. Chem. Soc.
,
1972
, vol.
94
pg.
5612
2.
Warshel
A.
,
Levitt
M.
,
J. Mol. Biol.
,
1976
, vol.
103
pg.
227
3.
Hohenberg
P.
,
Kohn
W.
,
Phys. Rev
,
1964
, vol.
B864
pg.
136
4.
Kohn
W.
,
Sham
L. J.
,
Phys. Rev
,
1965
, vol.
A1133
pg.
140
5.
Dunlap
B. I.
,
Connolly
J. W. D.
,
Sabin
J. R.
,
J. Chem. Phys.
,
1979
, vol.
71
pg.
4993
6.
Mintmire
J. W.
,
Dunlap
B. I.
,
Phys. Rev. A.
,
1982
, vol.
25
pg.
88
7.
Vahtras
O.
,
Almlöf
J.
,
Feyereisen
M. W.
,
Chem. Phys. Lett.
,
1993
, vol.
213
pg.
514
8.
Dunlap
B. I.
,
Rösch
N.
,
Trickey
S. B.
,
Mol. Phys.
,
2010
, vol.
108
pg.
3167
9.
Früchtl
H.
,
Kendall
R. A.
,
Harrison
R. J.
,
Dyall
K. G.
,
Int. J. Quantum Chem.
,
1997
, vol.
64
pg.
63
10.
Hamel
S.
,
Casida
M. E.
,
Salahub
D. R.
,
J. Chem. Phys.
,
2001
, vol.
114
pg.
7342
11.
Weigend
F.
,
Phys. Chem. Chem. Phys.
,
2002
, vol.
4
pg.
4285
12.
Polly
R.
,
Werner
H.-J.
,
Manby
F. R.
,
Knowles
P. J.
,
Mol. Phys.
,
2004
, vol.
102
pg.
2311
13.
Mejía-Rodríguez
D.
,
Köster
A. M.
,
J. Chem. Phys.
,
2014
, vol.
141
pg.
124114
14.
Manzer
S. F.
,
Epifanovsky
E.
,
Head-Gordon
M.
,
J. Chem. Theory Comput.
,
2015
, vol.
11
pg.
518
15.
Delesma
F. A.
,
Geudtner
G.
,
Mejía-Rodríguez
D.
,
Calaminici
P.
,
Köster
A. M.
,
J. Chem. Theory Comput.
,
2018
, vol.
14
pg.
5608
16.
A. M.
Köster
,
G.
Geudtner
,
A.
Alvarez-Ibarra
,
P.
Calaminici
,
M. E.
Casida
,
J.
Carmona-Espíndola
,
F. A.
Delesma
,
R.
Delgado-Venegas
,
V. D.
Domínguez
,
R.
Flores-Moreno
,
G. U.
Gamboa
,
A.
Goursot
,
T.
Heine
,
A.
Ipatov
,
A.
de la Lande
,
F.
Janetzko
,
J. M.
del Campo
,
N.
Pedroza Montero
,
L. G. M.
Pettersson
,
D.
Mejía-Rodríguez
,
J.
Reveles
,
J.
Vásquez-Pérez
,
A.
Vela
,
B.
Zúñiga-Gutiérrez
and
D. R.
Salahub
, deMon2k; Version 6.1.5; The deMon developers: Cinvestav,
Mexico City
,
2020
. See also: http://www.demon-software.com
17.
M. E.
Casida
,
C.
Daul
,
A.
Goursot
,
A. M.
Köster
,
L. G. M.
Pettersson
,
E.
Proynov
,
A.
St-Amant
,
D. R.
Salahub
,
H.
Duarte
,
N.
Godbout
,
J.
Guan
,
C.
Jamorski
,
M.
Leboeuf
,
V. G.
Malkin
,
O. L.
Malkina
,
F.
Sim
and
A.
Vela
,
deMon-KS Version 3.4, deMon Software
,
Montréal
,
1996
18.
Andzelm
J.
,
Wimmer
E.
,
J. Chem. Phys.
,
1992
, vol.
96
pg.
1280
19.
Sambe
H.
,
Felton
R. H.
,
J. Chem. Phys.
,
1975
, vol.
62
pg.
1122
20.
D. R.
Salahub
,
J.
Weber
,
A.
Goursot
,
A. M.
Köster
and
A.
Vela
,
Applied Density Functional Theory and the deMon Code, 1964–2004
, in
Theory and Applications of the Computational Chemistry: The First 40 Years
, ed. C. E. Dykstra, G. Frenking, K. S. Kim and G. Scuseria,
Elsevier
,
Amsterdam
,
2005
21.
Laikov
D. N.
,
Chem. Phys. Lett.
,
1997
, vol.
281
pg.
151
22.
Köster
A. M.
,
Reveles
J. U.
,
del Campo
J. M.
,
J. Chem. Phys.
,
2004
, vol.
121
pg.
3417
23.
Birkenheuer
U.
,
Gordienko
A. B.
,
Nasluzov
V. A.
,
Fuchs-Rohr
M. K.
,
Rösch
N.
,
Int. J. Quantum Chem.
,
2005
, vol.
102
pg.
743
24.
Bienvenu
A. V.
,
Knizia
G.
,
J. Chem. Theory Comput.
,
2018
, vol.
14
pg.
1297
25.
P.
Calaminici
,
A.
Alvarez-Ibarra
,
D.
Cruz-Olvera
,
V. D.
Dominguez-Soria
,
R.
Flores-Moreno
,
G. U.
Gamboa
,
G.
Geudtner
,
A.
Goursot
,
D.
Mejía-Rodríguez
,
D. R.
Salahub
,
B.
Zuniga-Gutierrez
and
A. M.
Köster
,
Auxiliary Density Functional Theory: From Molecules to Nanostructures
, in
Handbook of Computational Chemistry
, ed. J. Leszczynski, A. Kaczmarek-Kedziera, T. Puzyn, M. G. Papadopoulos, H. Reis and M. K. Shukla,
Springer International Publishing
, 2nd edn,
2017
, ch. 18, p. 795
26.
Choi
S.-C. T.
,
Paige
C. C.
,
Saunders
M. A.
,
SIAM J. Sci. Comput.
,
2011
, vol.
33
pg.
1810
27.
Eirola
T.
,
Nevanlinna
O.
,
Linear Algebra Appl.
,
1989
, vol.
121
pg.
511
28.
Pedroza-Montero
J. N.
,
Morales
J. L.
,
Geudtner
G.
,
Álvarez-Ibarra
A.
,
Calaminici
P.
,
Köster
A. M.
,
J. Chem. Theory Comput.
,
2020
, vol.
16
pg.
2965
29.
Mejía-Rodríguez
D.
,
Delgado-Venegas
R. I.
,
Calaminici
P.
,
Köster
A. M.
,
J. Chem. Theory Comput.
,
2015
, vol.
11
pg.
1493
30.
Köster
A. M.
,
J. Chem. Phys.
,
2003
, vol.
118
pg.
9943
31.
A.
Álvarez-Ibarra
,
Asymptotic Expansion of Molecular Integrals in Self-Consistent Auxiliary Density Functional Methods
, PhD thesis, Cinvestav,
2013
32.
Köster
A. M.
,
del Campo
J. M.
,
Janetzko
F.
,
Zuniga-Gutierrez
B.
,
J. Chem. Phys.
,
2009
, vol.
130
pg.
114106
33.
Calaminici
P.
,
Janetzko
F.
,
Köster
A. M.
,
Mejia-Olvera
R.
,
Zuniga-Gutierrez
B.
,
J. Chem. Phys.
,
2007
, vol.
126
pg.
044108
34.
Mejía-Rodríguez
D.
,
Huang
X.
,
del Campo
J. M.
,
Köster
A. M.
,
Adv. Quantum Chem.
,
2015
, vol.
71
pg.
41
35.
D.
Mejía-Rodríguez
,
Low-Order Scaling Methods for Auxiliary Density Functional Theory
, PhD thesis, Cinvestav,
2015
36.
F. A.
Delesma-Díaz
,
Range-Separated Hybrid Functionals in Auxiliary Density Functional Theory
, PhD thesis, Cinvestav,
2020
37.
Aquilante
F.
,
Pedersen
T.
,
Merés
A.
,
Koch
H.
,
J. Chem. Phys.
,
2006
, vol.
125
pg.
174101
38.
Boys
S. F.
,
Rev. Mod. Phys.
,
1960
, vol.
32
pg.
296
39.
Foster
J. M.
,
Boys
S. F.
,
Rev. Mod. Phys.
,
1960
, vol.
32
pg.
300
40.
Pedroza-Montero
J. N.
,
Delesma
F. A.
,
Morales
J. L.
,
Calaminici
P.
,
Köster
A. M.
,
J. Chem. Phys.
,
2020
, vol.
153
pg.
134112
41.
Obara
S.
,
Saika
A.
,
J. Chem. Phys.
,
1986
, vol.
84
pg.
3963
42.
de la Lande
A.
,
Alvarez-Ibarra
A.
,
Hasnaoui
K.
,
Cailliez
F.
,
Wu
X.
,
Mineva
T.
,
Cuny
J.
,
Calaminici
P.
,
López-Sosa
L.
,
Geudtner
G.
,
Navizet
I.
,
Garcia Iriepa
C.
,
Salahub
D. R.
,
Köster
A. M.
,
Molecules
,
2019
, vol.
24
pg.
1653
43.
Quintanar
C.
,
Caballero
R.
,
Köster
A. M.
,
Int. J. Quantum Chem.
,
2004
, vol.
96
pg.
483
44.
Salahub
D. R.
,
Noskov
S. Y.
,
Lev
B.
,
Zhang
R.
,
Ngo
V.
,
Goursot
A.
,
Calaminici
P.
,
Köster
A. M.
,
Alvarez-Ibarra
A.
,
Mejía-Rodríguez
D.
,
Řezáč
J.
,
Cailliez
F.
,
de la Lande
A.
,
Molecules
,
2015
, vol.
20
pg.
4780
45.
Alvarez-Ibarra
A.
,
Köster
A. M.
,
Zhang
R.
,
Salahub
D. R.
,
J. Chem. Theory Comput.
,
2012
, vol.
8
pg.
4232
46.
Alvarez-Ibarra
A.
,
Köster
A. M.
,
J. Chem. Phys.
,
2013
, vol.
139
pg.
024102
47.
Dirac
P. A. M.
,
Proc. Cambridge Philos. Soc.
,
1930
, vol.
26
pg.
376
48.
Vosko
S. H.
,
Wilk
L.
,
Nusair
M.
,
Can. J. Phys.
,
1980
, vol.
58
pg.
1200
49a.
Perdew
J. P.
,
Burke
K.
,
Ernzerhof
M.
,
Phys. Rev. Lett.
,
1996
, vol.
77
pg.
3865
49b.
Perdew
J. P.
,
Burke
K.
,
Ernzerhof
M.
,
Phys. Rev. Lett.
,
1997
, vol.
78
pg.
1396(E)
50.
Flores-Moreno
R.
,
Alvarez-Mendez
R. J.
,
Vela
A.
,
Köster
A. M.
,
J. Comput. Chem.
,
2006
, vol.
27
pg.
1009
51.
Jardillier
N.
,
Goursot
A.
,
Chem. Phys. Lett.
,
2008
, vol.
454
pg.
65
52.
Zuniga-Gutierrez
B.
,
Köster
A. M.
,
Mol. Phys.
,
2016
, vol.
114
pg.
1026
53.
Shedge
S. V.
,
Carmona-Espíndola
J.
,
Pal
S.
,
Köster
A. M.
,
J. Phys. Chem. A
,
2010
, vol.
114
pg.
2357
54.
Goursot
A.
,
Mineva
T.
,
Kevorkyants
R.
,
Talbi
D.
,
J. Chem. Theory Comput.
,
2007
, vol.
3
pg.
755
55.
Jurečka
P.
,
Černý
J.
,
Hobza
P.
,
Salahub
D. R.
,
J. Comput. Chem.
,
2007
, vol.
28
pg.
555
56.
Grimme
S.
,
Anthony
J.
,
Ehrlich
S.
,
Krieg
H.
,
J. Chem. Phys.
,
2010
, vol.
132
pg.
154104
57.
Pérez-Figueroa
S. E.
,
Calaminici
P.
,
Köster
A. M.
,
J. Phys. Chem. A
,
2019
, vol.
123
pg.
4565
58.
Perdew
J. P.
,
Ernzerhof
M.
,
Burke
K.
,
J. Chem. Phys.
,
1996
, vol.
105
pg.
9982
59.
Adamo
C.
,
Barone
V.
,
J. Chem. Phys.
,
1999
, vol.
110
pg.
6158
60.
Dunning
T. H.
,
J. Chem. Phys.
,
1989
, vol.
90
pg.
1007
61.
Alvarez-Ibarra
A.
,
Köster
A. M.
,
Mol. Phys.
,
2015
, vol.
113
pg.
3128
62.
A.
Álvarez-Ibarra
,
P.
Calaminici
,
A.
Goursot
,
C. Z.
Gómez-Castro
,
R.
Grande-Aztatzi
,
T.
Mineva
,
D. R.
Salahub
,
J. M.
Vásquez-Pérez
,
A.
Vela
,
B.
Zuniga-Gutierrez
and
A. M.
Köster
,
First Principles Computational Biochemistry with deMon2k
, in
Frontiers in Computational Chemistry
, ed. Zaheer-Ul_Haq and J. D. Madura,
Bentham Science Publishers (eBook)
,
2015
,
vol. 2
, p. 281
63.
Krack
M.
,
Köster
A. M.
,
J. Chem. Phys.
,
1998
, vol.
108
pg.
3226
64.
Köster
A. M.
,
Flores-Moreno
R.
,
Reveles
J. U.
,
J. Chem. Phys.
,
2004
, vol.
121
pg.
681
65.
A. M.
Köster
,
G.
Geudtner
,
G. U.
Gamboa
,
A.
Alvarez-Ibarra
,
P.
Calaminici
,
R.
Flores-Moreno
,
A.
Goursot
,
A.
de la Lande
,
D.
Mejia-Rodriguez
,
T.
Mineva
,
L. G. M.
Pettersson
,
J. M.
Vasquez-Perez
,
B.
Zuniga-Gutierrez
,
S. B.
Trickey
and
D. R.
Salahub
,
The deMon2k Users’ Guide, Version 5.0
, Cinvestav,
Mexico-City
,
2018
, See also: http://www.demon-software.com
66.
Jorgensen
W. L.
,
Maxwell
D. S.
,
Tirado-Rives
J.
,
J. Am. Chem. Soc.
,
2015
, vol.
118
pg.
11225
67.
Hornak
V.
,
Abel
R.
,
Okur
A.
,
Strockbine
B.
,
Roitberg
A.
,
Simmerling
C.
,
Proteins
,
2006
, vol.
65
pg.
712
68.
Pulay
P.
,
Mol. Phys.
,
1969
, vol.
17
pg.
197
69.
S. T.
Epstein
,
The Variation Method in Quantum Chemistry
,
Academic Press
London
,
1974
70.
Dalgarno
A.
,
Stewart
A. L.
,
Proc. R. Soc. London, Ser. A
,
1956
, vol.
238
pg.
269
71.
Born
M.
,
Oppenheimer
J. R.
,
Ann. Phys.
,
1927
, vol.
389
pg.
457
72.
R.
Fletcher
,
Practical Methods of Optimization
,
John Wiley & Sons
,
Chichester, N.Y
.,
1980
73.
J.
Nocedal
and
S. J.
Wright
,
Numerical Optimization
,
Springer-Verlag
,
New York
,
1999
74.
Schlegel
H. B.
,
J. Comput. Chem.
,
1982
, vol.
3
pg.
214
75.
K. P.
Lawley
,
Advances in Chemical Physics: Ab Initio Methods in Quantum Chemistry Part I
,
John Wiley & Sons
,
Chichester, N.Y.
,
1987
76.
D. R.
Yarkony
,
Modern Electronic Structure Theory: Part I
,
World Scientific
,
Singapore
,
1995
77.
P. v. R.
Schleyer
,
Encyclopedia of Computational Chemistry
,
Wiley
,
New York
,
1998
78.
Farkas
Ö.
,
Schlegel
H. B.
,
J. Chem. Phys.
,
1998
, vol.
109
pg.
7100
79.
Farkas
Ö.
,
Schlegel
H. B.
,
J. Chem. Phys.
,
1999
, vol.
111
pg.
10806
80.
Farkas
Ö.
,
Schlegel
H. B.
,
Phys. Chem. Chem. Phys.
,
2002
, vol.
4
pg.
11
81.
Schlegel
H. B.
,
J. Comput. Chem.
,
2003
, vol.
24
pg.
1514
82.
F.
Jensen
,
Introduction to Computational Chemistry
,
John Wiley & Sons
,
Chichester, NY
,
2007
83.
Schlegel
H. B.
,
Wiley Interdiscip. Rev. Comput. Mol. Sci.
,
2011
, vol.
1
pg.
790
84.
E. G.
Lewars
,
Computational Chemistry: Introduction to the Theory and Applications of Molecular and Quantum Mechanics
,
Springer-Verlag
,
Dordrecht
,
2011
85.
del Campo
J. M.
,
Köster
A. M.
,
Croat. Chem. Acta
,
2009
, vol.
82
pg.
283
86.
Levenberg
K.
,
Quart. Appl. Math.
,
1944
, vol.
2
pg.
164
87.
Marquardt
D. W.
,
SIAM J.
,
1963
, vol.
11
pg.
431
88.
Broyden
C. G.
,
J. Inst. Math. Appl.
,
1970
, vol.
6
pg.
76
89.
Fletcher
R.
,
Comput. J.
,
1970
, vol.
13
pg.
317
90.
Goldfarb
D.
,
Math. Comput.
,
1970
, vol.
24
pg.
23
91.
Shanno
D. F.
,
Math. Comput.
,
1970
, vol.
24
pg.
647
92.
Culot
P.
,
Dive
G.
,
Nguyen
V. H.
,
Ghuysen
J. M.
,
Theor. Chim. Acta
,
1992
, vol.
82
pg.
189
93.
Bofill
J. M.
,
J. Comput. Chem.
,
1994
, vol.
15
pg.
1
94.
del Campo
J. M.
,
Köster
A. M.
,
J. Chem. Phys.
,
2008
, vol.
129
pg.
024107
95.
Powell
M. J. D.
,
Math. Program.
,
1971
, vol.
1
pg.
26
96.
M. E.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
,
Oxford University Press
,
Oxford
,
2010
97.
Trotter
H. F.
,
Proc. Am. Math. Soc.
,
1959
, vol.
10
pg.
545
98.
Swope
W. C.
,
Anderson
H. C.
,
Berens
P. H.
,
Wilson
K. R.
,
J. Chem. Phys.
,
1982
, vol.
76
pg.
637
99.
Gamboa
G. U.
,
Vázquez-Pérez
J. M.
,
Calaminici
P.
,
Köster
A. M.
,
Int. J. Quantum Chem.
,
2010
, vol.
110
pg.
2172
100.
Nosé
S.
,
J. Chem. Phys.
,
1984
, vol.
81
pg.
511
101.
Martyna
G. J.
,
Klein
M.
,
Tuckerman
M.
,
J. Chem. Phys.
,
1992
, vol.
97
pg.
2635
102.
Cornell
W. D.
,
Cieplak
P.
,
Bayly
C. I.
,
Gould
I. R.
,
Merz Jr
K. M.
,
Ferguson
D. M.
,
Spellmeyer
D. C.
,
Fox
T.
,
Caldwell
J. W.
,
Kollman
P. A.
,
J. Am. Chem. Soc.
,
1995
, vol.
117
pg.
5179
103.
Mineva
T.
,
Russo
N.
,
Int. J. Quantum Chem.
,
1997
, vol.
61
pg.
665
104.
Cheatham III
T. E.
,
Cieplak
P.
,
Kollman
P. A.
,
J. Biomol. Struct. Dyn.
,
1999
, vol.
16
pg.
845
105.
Wang
J.
,
Cieplak
P.
,
Kollman
P. A.
,
J. Comput. Chem.
,
2000
, vol.
21
pg.
1049
106.
Martyna
G. J.
,
Tobias
D.
,
Klein
M.
,
J. Chem. Phys.
,
1994
, vol.
101
pg.
4177
107.
Tuckerman
M.
,
Alejandre
J.
,
López-Rendón
R.
,
Jochim
A. L.
,
J. Phys. A: Math. Gen.
,
2006
, vol.
39
pg.
5629
108.
Bredow
T.
,
Geudtner
G.
,
Jug
K.
,
J. Comput. Chem.
,
2001
, vol.
22
pg.
89
109.
Janetzko
F.
,
Köster
A. M.
,
Salahub
D. R.
,
J. Chem. Phys.
,
2008
, vol.
128
pg.
024102
110.
Janetzko
F.
,
Bredow
T.
,
Geudtner
G.
,
Köster
A. M.
,
J. Comput. Chem.
,
2008
, vol.
29
pg.
2295
111.
de Leeuw
S. W.
,
Perram
J. W.
,
Smith
E. R.
,
Proc. R. Soc. Lond. A
,
1980
, vol.
373
pg.
27
112.
Toukmaji
A. Y.
,
Board Jr
J. A.
,
Comput. Phys. Commun.
,
1996
, vol.
95
pg.
73
113.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
, 2nd edn,
2017
,
Oxford University Press
,
Oxford
114.
Nijboer
B.
,
De Wette
F.
,
Physica
,
1957
, vol.
23
pg.
309
115.
Rackers
J. A.
,
Laury
M. L.
,
Lu
C.
,
Wang
Z.
,
Lagardère
L.
,
Piquemal
J.-P.
,
Ren
P.
,
Ponder
J. W.
,
J. Chem. Theory Comput.
,
2018
, vol.
14
pg.
5273
116.
Kawashima
Y.
,
Ishimura
K.
,
Shiga
M.
,
J. Chem. Phys.
,
2019
, vol.
150
pg.
124103
117.
Delgado-Venegas
R. I.
,
Mejía-Rodríguez
D.
,
Flores-Moreno
R.
,
Calaminici
P.
,
Köster
A. M.
,
J. Chem. Phys.
,
2016
, vol.
145
pg.
224103
118.
R.
Parr
and
W.
Yang
,
Density-Functional Theory of Atoms and Molecules
,
Oxford University Press
,
New York
,
1989
119.
Köster
A. M.
,
J. Chem. Phys.
,
1996
, vol.
104
pg.
4114
120.
Diercksen
G.
,
McWeeny
R.
,
J. Chem. Phys.
,
1966
, vol.
44
pg.
3554
121.
McWeeny
R.
,
Diercksen
G.
,
J. Chem. Phys.
,
1968
, vol.
49
pg.
4852
122.
Dodds
J. L.
,
McWeeny
R.
,
Raynes
W. T.
,
Riley
J. P.
,
Mol. Phys.
,
1977
, vol.
33
pg.
611
123.
Dodds
J. L.
,
McWeeny
R.
,
Sadlej
A. J.
,
Mol. Phys.
,
1977
, vol.
34
pg.
1779
124.
R.
McWeeny
,
Methods of Molecular Quantum Mechanics
, 2nd edn,
2001
,
Academic Press
,
London
125.
Flores-Moreno
R.
,
Köster
A. M.
,
J. Chem. Phys.
,
2008
, vol.
128
pg.
134105
126.
Carmona-Espíndola
J.
,
Flores-Moreno
R.
,
Köster
A. M.
,
J. Chem. Phys.
,
2010
, vol.
133
pg.
084102
127.
Carmona-Espíndola
J.
,
Flores-Moreno
R.
,
Köster
A. M.
,
Int. J. Quantum Chem.
,
2012
, vol.
112
pg.
3461
128.
Flores-Moreno
R.
,
Melin
J.
,
Ortiz
J. V.
,
Merino
G.
,
J. Chem. Phys.
,
2008
, vol.
129
pg.
224105
129.
Flores-Moreno
R.
,
J. Chem. Theory Comput.
,
2010
, vol.
6
pg.
48
130.
Neumaier
A.
,
SIAM Rev.
,
1998
, vol.
40
pg.
636
131.
Yang
U. M.
,
Gallivan
K. A.
,
Appl. Numer. Math.
,
1995
, vol.
19
pg.
287
132.
Geudtner
G.
,
Janetzko
F.
,
Köster
A. M.
,
Vela
A.
,
Calaminici
P.
,
J. Comput. Chem.
,
2006
, vol.
27
pg.
483
133.
Geudtner
G.
,
Calaminici
P.
,
Carmona-Espíndola
J.
,
del Campo
J. M.
,
Dominguez-Soria
V. D.
,
Flores-Moreno
R.
,
Gamboa
G. U.
,
Goursot
A.
,
Köster
A. M.
,
Reveles
J. U.
,
Mineva
T.
,
Vasquez-Perez
J. M.
,
Vela
A.
,
Zuniga-Gutierrez
B.
,
Salahub
D. R.
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
,
2012
, vol.
2
pg.
548
134.
E. B.
Wilson
,
J. C.
Decius
and
P. C.
Cross
,
Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra
,
Dover, New York
,
1980
135.
Neugebauer
J.
,
Reiher
M.
,
Kind
C.
,
Hess
B. A.
,
J. Comput. Chem.
,
2002
, vol.
23
pg.
895
136.
Delgado-Venegas
R. I.
,
Calaminici
P.
,
Köster
A. M.
,
Mol. Phys.
,
2019
, vol.
117
pg.
1367
137.
Cui
Q.
,
Karplus
M.
,
J. Chem. Phys.
,
2000
, vol.
112
pg.
1133
138.
Cui
Q.
,
Karplus
M.
,
J. Phys. Chem. B
,
2000
, vol.
104
pg.
3721
139.
Zuniga-Gutierrez
B.
,
Geudtner
G.
,
Köster
A. M.
,
J. Chem. Phys.
,
2011
, vol.
134
pg.
124108
140.
Zuniga-Gutierrez
B.
,
Geudtner
G.
,
Köster
A. M.
,
J. Chem. Phys.
,
2012
, vol.
137
pg.
094113
141.
Zuniga-Gutierrez
B.
,
Camacho-Gonzalez
M.
,
Simon-Bastida
P.
,
Bendana-Castillo
A.
,
Calaminici
P.
,
Köster
A. M.
,
J. Phys. Chem. A
,
2014
, vol.
119
pg.
1469
142.
Zuniga-Gutierrez
B.
,
Camacho-Gonzalez
M.
,
Simon-Bastida
P.
,
Bendana-Castillo
A.
,
Calaminici
P.
,
Köster
A. M.
,
J. Chem. Phys.
,
2015
, vol.
143
pg.
104103
143.
Zuniga-Gutierrez
B.
,
Flores-Moreno
R.
,
Medel-Juarez
V.
,
Varona
A.
,
González-Ramírez
H. N.
,
J. Chem. Phys.
,
2020
, vol.
152
pg.
014105
144.
Ditchfield
R.
,
J. Chem. Phys.
,
1972
, vol.
56
pg.
5688
145.
B. A.
Zuniga-Gutierrez
,
Magnetic Shielding, Spin–Spin Coupling Constant and Magnetizability Tensors Calculated with Auxiliary Density Perturbation Theory
, PhD thesis, Cinvestav,
2011
146.
Ipatov
A.
,
Fouqueau
A.
,
Perez del Valle
C.
,
Cordova
F.
,
Casida
M. E.
,
Köster
A. M.
,
Vela
A.
,
J. Mol. Struct.: THEOCHEM
,
2006
, vol.
762
pg.
179
147.
M. E.
Casida
,
Time-dependent density functional response theory for molecules
, in
Recent Advances in Density Functional Methods (Part I)
, ed. D. P. Chong,
World Scientific
,
Singapore
,
1995
, ch. 5, p. 155
148.
Jacquemin
D.
,
Perpète
E. A.
,
Scuseria
G. E.
,
Ciofini
I.
,
Adamo
C.
,
J. Chem. Theory. Comput.
,
2008
, vol.
4
pg.
123
149.
Shao
Y.
,
Mei
Y.
,
Sundholm
D.
,
Kaila
V. R. I.
,
J. Chem. Theory Comput.
,
2020
, vol.
16
pg.
587
150.
Jamorski
C.
,
Casida
M. E.
,
Salahub
D. R.
,
J. Chem. Phys.
,
1996
, vol.
104
pg.
5134
151.
Carmona-Espíndola
J.
,
Köster
A. M.
,
Can. J. Chem.
,
2013
, vol.
91
pg.
795
152.
Hirata
S.
,
Head-Gordon
M.
,
Chem. Phys. Lett.
,
1999
, vol.
314
pg.
291
153.
Davidson
E.
,
J. Comput. Phys.
,
1975
, vol.
17
pg.
87
154.
Cisneros
G.
,
Berrondo
M.
,
Bunge
C. F.
,
Comput. Chem.
,
1986
, vol.
10
pg.
281
Close Modal

or Create an Account

Close Modal
Close Modal