
==== Front
J Mol Model
J Mol Model
Journal of Molecular Modeling
1610-2940
0948-5023
Springer Berlin Heidelberg Berlin/Heidelberg

39115689
6091
10.1007/s00894-024-06091-z
Original Paper
Exchange-correlation kernel for perturbation dependent auxiliary functions in auxiliary density perturbation theory
Hernández-Segura Luis I. luisi.segura@cinvestav.mx

1
Olvera-Rubalcava Flor A. 1
Flores-Moreno Roberto 2
Calaminici Patrizia 1
Köster Andreas M. akoster@cinvestav.mx

1
1 grid.512574.0 Chemistry Department, CINVESTAV, Av. Instituto Politecnico Nacional 2508, Col. San Pedro Zacatenco, Del. Gustavo A. Madero, Mexico City C.P. 07360 Mexico
2 https://ror.org/043xj7k26 grid.412890.6 0000 0001 2158 0196 Departamento de Química, Universidad de Guadalajara, Blvd. Gral. Marcelino García Barragán 1421, Guadalajara, Jal. C.P. 44430 Mexico
8 8 2024
8 8 2024
2024
30 9 30221 5 2024
23 7 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Context

Analytic exchange-correlation kernel formulations are of the outermost importance for density functional theory (DFT) perturbation calculations. In this paper, the working equation for the exchange-correlation kernel of the generalized gradient approximation (GGA) for perturbation dependent auxiliary functions is derived and discussed in the framework of auxiliary density functional theory (ADFT). The presented new formulation is extended to the unrestricted approach, too. A comprehensive discussion of the implementation of the GGA ADFT kernel, using either the native exchange-correlation functional implementations in deMon2k or the ones from the LibXC library, is given. Calculations with analytic exchange-correlation kernels are compared to their finite difference counterparts. The obtained results are in quantitative agreement. Nevertheless, analytic GGA ADFT kernel implementations show substantial improvement in the computational performance. Similar results are reported for analytic second derivatives of effective core potential (ECP) and model core potential (MCP) matrix elements when compared to their finite difference counterparts in molecular frequency analyses.

Method

All calculations are performed in the framework of ADFT as implemented in deMon2k. In the ADFT analytic frequency calculations, auxiliary density perturbation theory was used. The underlying two-center exchange-correlation kernel matrix elements are calculated by numerical integration either with analytic or finite difference kernel expressions. Validation calculations are performed with the VWN and PBE functionals employing DFT-optimized DZVP basis sets in conjunction with automatically generated GEN-A2 auxiliary density function sets. In the (Pt3Cu)n cluster benchmark calculations, the RPBE functional was used. For Pt atoms, the quasi-relativistic LANL2DZ effective core potential with the corresponding valence basis set was employed, whereas for Cu atoms, the all-electron DFT-optimized TZVP basis was applied. The auxiliary density was expanded by the automatically generated GEN-A2* auxiliary function set. We run all benchmark calculations in parallel on 24 cores.

Keywords

ADFT
Exchange-correlation functional derivatives
deMon2k
LibXC
Analytic frequency calculation
Finite difference exchange-correlation kernel
CONAHCYTA1-S-11929 A1-S-11929 CONAHCYT PhD fellowship798159 ECOS321168 321168 321168 issue-copyright-statement© Springer-Verlag GmbH Germany, part of Springer Nature 2024
==== Body
pmcIntroduction

Over the last three decades, Kohn-Sham density functional theory (DFT) [1, 2] methods have become the workhorse in computational chemistry and materials science [3]. Their excellent performance-to-accuracy ratio permits reliable energy and structure predictions for systems with several hundreds of atoms. Although there exists still an accuracy gap to the targeted chemical accuracy of 1 kcal/mol, it is foreseeable that this gap will close with the further development of Kohn-Sham DFT. To this end, it is important to improve density functional approximations (DFAs) such that the overall computational performance is not jeopardized. A possible framework for such developments is auxiliary density functional theory (ADFT) [4] which is based on the variational density fitting for Coulomb [5, 6] and Fock [7, 8] energies and uses the resulting auxiliary density for the evaluation of the exchange-correlation contributions [9, 10]. As a result, ADFT local density approximation (LDA) and generalized gradient approximation (GGA) as well as hybrid calculations have a formal cubic scaling with respect to the number of basis functions. No four-center electron repulsion integral (ERI) evaluations nor basis function product evaluations at grid points are needed in ADFT. Despite this enormous computational simplification, ADFT provides for a given DFA the same accuracy in terms of structure parameters, relative energies, and harmonic frequencies, to name a few, as corresponding Kohn-Sham calculations. Therefore, ADFT is particularly well suited for the reliable optimization of complex nanometric systems like large transition metal clusters [11, 12] or Born-Oppenheimer molecular dynamics simulations [13, 14].

Because the ADFT energy expression is variational, analytic energy derivatives can be straightforwardly calculated, e.g., energy gradients by the 2n+1 theorem of perturbation theory [15]. This permits an alternative perturbation theory formulation to the commonly employed coupled-perturbed Kohn-Sham (CPKS) approach [16]. We have named this alternative approach auxiliary density perturbation theory (ADPT) [17, 18] to emphasize its close connection to ADFT. Since its initial formulation, ADPT has been used to calculate polarizabilities [19–23], Fukui functions [24, 25], electron binding energies [26], alchemical derivatives [27], and various magnetic properties [28–34]. In all cases, CPKS equivalent results were obtained with a much reduced computational demand. In fact, for LDA and GGA perturbation calculations, timings comparable to the ADFT self-consistent field (SCF) approach are observed. More recently, this has been extended to time-dependent ADFT excited state calculations obtaining similar computational performance [35–37]. One reason for the improved computational performance of ADPT arises from the numerical exchange-correlation kernel calculations, i.e., the second functional derivatives of the exchange-correlation energy. In ADPT, these matrix elements contain only auxiliary functions from two centers. As a result, the corresponding kernel matrix has only the dimension of the number of auxiliary functions. Key ingredients for the efficient exchange-correlation kernel calculations are analytic formulas for the various LDA and GGA functionals. Whereas the LDA ADFT kernel formulas are rather straightforward to derive, the corresponding GGA formulas are more involved and deserve some attention. For perturbation-independent auxiliary functions, the generic formulas for closed-shell ADFT kernel calculations have been published a few years ago [38]. However, the recent development of ADFT second analytic energy derivatives [39] employing ADPT requires corresponding kernel calculations for perturbation dependent auxiliary functions. Again, computationally efficient formulations are needed because we aim for frequency analyses of systems with many hundreds to several thousands of atoms. Such calculations are not only mandatory for the characterization of the optimized structures but also for initial Hessian matrix calculations, either for structure minimizations or transition state optimizations. Again, the focus is on GGA kernel formulas which we will derive in this work. To keep the discussion most general, we present the GGA kernel working equations in the framework of the unrestricted formulation. Furthermore, we also validate the recently implemented LibXC interface [40] of deMon2k in the context of these kernel calculations. Because we found in some transition metal cluster frequency analyses a computational bottleneck arising from the finite difference calculations of integrals for effective core potential (ECP) and model core potential (MCP) second derivatives, we also report the corresponding analytic second derivative formulas that were newly implemented into deMon2k [41] in the framework of this study.

The paper is organized as follows. The subsequent theory section contains four subsections. After a brief introduction of second analytic ADFT energy derivatives, the generic working formulas for the GGA potential and kernel calculations for perturbation dependent auxiliary functions are presented. The two following subsections discuss the finite difference GGA ADFT kernel formulations that can be employed if explicit kernel expressions for a GGA functional are missing, and the analytical second derivative formulas for ECP and MCP matrix elements. Section 3 provides the computational details for the validation and benchmark calculations that we present and discuss in the following section. All calculations are performed with the deMon2k code using either the native or LibXC exchange-correlation functional implementations. Finally, in Section 5, the conclusions are summarized.

Theory

In the linear combination of Gaussian-type orbitals (LCGTO) approximation, neglecting for the sake of simplicity explicit spin dependency, the second analytic ADFT energy derivatives with respect to atomic coordinates λ and η are given by the following [39]:1 E(λη)≡∂2E∂λ∂η=∑μ,νPμν(η)Hμν(λ)+∑k¯⟨μν||k¯⟩(λ)(xk¯+zk¯)+∑μ,νPμνHμν(λη)+∑k¯⟨μν||k¯⟩(λη)(xk¯+zk¯)+∑μ,ν∑k¯Pμν⟨μν||k¯⟩(λ)xk¯(η)+zk¯(η)+∑k¯xk¯(η)⟨k¯(λ)|vxc[ρ~]⟩+∑k¯xk¯⟨k¯(λη)|vxc[ρ~]⟩+∑k¯xk¯⟨k¯(λ)|vxc(η)[ρ~]⟩-∑μ,νWμν(η)Sμν(λ)-∑μ,νWμνSμν(λη)-∑k¯,l¯Gk¯l¯(λη)xl¯12xk¯+zk¯-∑k¯,l¯Gk¯l¯(λ)(xk¯+zk¯)xl¯(η)-∑k¯,l¯Gk¯l¯(λ)zk¯(η)xl¯

In Eq. (1), superscripts in parentheses denote derivatives with respect to the corresponding atomic coordinates λ and η. Thus, E(λη) corresponds to a Hessian matrix element. The Greek letters μ and ν denote (contracted) atomic GTOs, whereas the Latin letters with a bar, k¯ and l¯, represent primitive atom-centered Hermite Gaussian auxiliary functions. Furthermore, Pμν is an element of the (closed-shell) density matrix, and Hμν is an element of the core matrix, incorporating kinetic, nuclear attraction, and any other external potential contributions, e.g., external electric or magnetic fields. The symbol || in the three-center ERI shorthand notation denotes the Coulomb operator 1/|r1-r2|. It also separates functions from electron 1, in the bra, from those of electron 2, in the ket. The Coulomb and exchange-correlation fitting coefficients are denoted by xk¯ and zk¯, respectively. The matrix elements Sμν, Wμν and Gk¯l¯ belong to the overlap, energy-weighted density, and Coulomb matrices. Superscripts on matrix elements or fitting coefficients indicate derivatives with respect to corresponding atomic coordinates.

Analytic GGA exchange-correlation potential

Of particular interest to our further discussion are the matrix elements in Eq. (1) that contain vxc and vxc(η), which are the auxiliary density exchange-correlation potential and its derivatives that result in the corresponding kernel expressions. Note that in ADFT, all these expressions are evaluated with the auxiliary density obtained from the variational fitting of Coulomb and Fock energies. In the case of a GGA functional, the exchange-correlation energy takes the following form:2 Exc[ρ~,γ~]=∫exc(ρ~,γ~)dr

In Eq. (2), exc(ρ~,γ~) denotes the density-weighted exchange-correlation energy density depending on the auxiliary density, ρ~(r), and its gradient square, γ~(r). For simplicity of notation, we omit the explicit position dependencies from the exchange-correlation energy density, potential, and kernel throughout the discussion. In the unrestricted approach, the auxiliary density is expanded as follows:3 ρ~(r)=ρ~α(r)+ρ~β(r)=∑k¯xk¯αk¯(r)+∑k¯xk¯βk¯(r)

The xk¯α and xk¯β are the spin-polarized Coulomb fitting coefficients, obtained from separate spin-dependent fitting equations. The density gradient corrections are included by the scalar γ~στ(r) field, with σ and τ being labels for the spin, either α or β:4 γ~στ(r)=∇ρ~σ(r)·∇ρ~τ(r)

Therefore, the explicit form of the unrestricted exchange-correlation energy in Eq. (2) is given by the following:5 Exc[ρ~,γ~]=Exc[ρ~α,ρ~β,γ~αα,γ~αβ,γ~βα,γ~ββ]=∫excρ~α,ρ~β,γ~αα,γ~αβ,γ~βα,γ~ββdr

The first derivative of the unrestricted exchange-correlation energy, Eq. (5), with respect to the atomic coordinate λ yields the following:6 ∂Exc[ρ~,γ~]∂λ=∑σα,β∫∂exc(ρ~,γ~)∂ρ~σ(r)∂ρ~σ(r)∂λdr+∑σα,β∑ux,y,z∫∂exc(ρ~,γ~)∂ρ~uσ(r)∂ρ~uσ(r)∂λdr

For the sake of clarity of presentation, we use explicit summation in Eq. (6). The shorthand notation ρ~uσ(r) stands for ∂ρ~σ/∂u with u being x, y, and z. Thus, ρ~uσ(r) denotes the components of the auxiliary density gradient, ∇ρ~σ(r), with respect to the electronic coordinates. Expanding the auxiliary density and density gradient derivatives as,7 ∂ρ~σ(r)∂λ=∑k¯xk¯σ(λ)k¯(r)+∑k¯xk¯σk¯(λ)(r)

8 and∂ρ~uσ(r)∂λ=∑k¯xk¯σ(λ)∂k¯(r)∂u+∑k¯xk¯σ∂k¯(λ)(r)∂u,

yields the following:9 ∂Exc[ρ~,γ~]∂λ=∑σα,β∑k¯xk¯σ(λ)∫∂exc(ρ~,γ~)∂ρ~σ(r)+∑ux,y,z∂exc(ρ~,γ~)∂ρ~uσ(r)∂∂uk¯(r)dr+∑σα,β∑k¯xk¯σ∫∂exc(ρ~,γ~)∂ρ~σ(r)+∑ux,y,z∂exc(ρ~,γ~)∂ρ~uσ(r)∂∂uk¯(λ)(r)dr.

Defining the unrestricted GGA ADFT exchange-correlation potential operator [38],10 vxcσ[ρ~,γ~]≡∂exc(ρ~,γ~)∂ρ~σ(r)+∑ux,y,z∂exc(ρ~,γ~)∂ρ~uσ(r)∂∂u,

allows the following shorthand notation for the exchange-correlation energy derivative of Eq. (6):11 ∂Exc[ρ~,γ~]∂λ=∑σα,β∑k¯xk¯σ(λ)k¯|vxcσ[ρ~,γ~]+∑σα,β∑k¯xk¯σk¯(λ)|vxcσ[ρ~,γ~]

Note that the first term of Eq. (11) will be absorbed in the Pulay term of the ADFT energy gradients [42] and, therefore, will not be explicitly calculated. Terms analogous to the second sum of Eq. (11) appear in the Hessian matrix elements given in Eq. (1) in the form of the k¯(λ)|vxc[ρ~] and k¯(λη)|vxc[ρ~] matrix elements. Whereas the first terms are explicitly evaluated in the response ADPT calculations, the second term’s contributions to the skeleton Hessian matrix are indirectly evaluated by translational invariance. This permits the use of moderate grids for the numerical integration in the Hessian matrix calculation even without the explicit inclusion of atomic weight function derivatives [43].

For the implementation in deMon2k, it is more convenient to formulate the exchange-correlation potential of Eq. (10) in terms of ρ~σ(r) and γ~στ(r). To this end, we apply the chain rule to Eq. (10) obtaining the following:12 vxcσ[ρ~,γ~]≡∂exc(ρ~,γ~)∂ρ~σ(r)+∑σ′,τ′α,β∑ux,y,z∂exc(ρ~,γ~)∂γ~σ′τ′(r)∂γ~σ′τ′(r)∂ρ~uσ(r)∂∂u

For the derivative of γ~σ′τ′(r) with respect to ρ~uσ(r) holds:13 ∂γ~σ′τ′(r)∂ρ~uσ(r)=ρ~u′(r)δστ′+ρ~u′(r)δσσ′

Inserting Eq. (13) into Eq. (12) and the result into Eq. (11) yields the explicit form of the unrestricted exchange-correlation energy derivative as implemented in deMon2k:14 ∂Exc[ρ~,γ~]∂λ=∑σα,β∑k¯xk¯σ(λ)∫[∂exc(ρ~,γ~)∂ρ~σ(r)+2∑σ′α,β∑ux,y,z∂exc(ρ~,γ~)∂γ~σσ′(r)ρ~u′(r)∂∂u]k¯(r)dr+∑σα,β∑k¯xk¯σ∫[∂exc(ρ~,γ~)∂ρ~σ(r)+2∑σ′α,β∑ux,y,z∂exc(ρ~,γ~)∂γ~σσ′(r)ρ~u′(r)∂∂u]k¯(λ)(r)dr

The corresponding exchange-correlation potential, vxcσ[ρ~,γ~], is given by the following:15 vxcσ[ρ~,γ~]≡∂exc(ρ~,γ~)∂ρ~σ(r)+2∑σ′α,β∑ux,y,z∂exc(ρ~,γ~)∂γ~σσ′(r)ρ~u′(r)∂∂u

Analytic GGA exchange-correlation kernel

We now turn to the analytic GGA ADFT exchange-correlation kernel calculation in ADPT. To this end, we derive Eq. (11) with respect to a second atomic coordinate keeping in mind that the first term of Eq. (11) is absorbed by the Pulay relation. Thus, only the second sum of Eq. (11) needs to be explicitly considered:16 ∂2Exc[ρ~,γ~]∂η∂λ|explicit=∂∂η∑σα,β∑k¯xk¯σk¯(λ)|vxcσ[ρ~,γ~]=∑σα,β∑k¯xk¯σ(η)k¯(λ)|vxcσ[ρ~,γ~]+∑σα,β∑k¯xk¯σk¯(λη)|vxcσ[ρ~,γ~]+∑σα,β∑k¯xk¯σk¯(λ)|∂vxcσ[ρ~,γ~]∂η

The first two terms of Eq. (16) are already discussed in Section 2.1 and can be calculated by the exchange-correlation potential expressions given above. Thus, we focus in the following only on the third term of Eq. (16).

To proceed, we derive the exchange-correlation potential with respect to the atomic coordinate η:17 ∂vxcσ[ρ~,γ~]∂η=∑τα,β∫∂vxcσ(ρ~,γ~)∂ρ~τ(r′)∂ρ~τ(r′)∂ηdr′+∑τα,β∑vx,y,z∫∂vxcσ(ρ~,γ~)∂ρ~vτ(r′)∂ρ~vτ(r′)∂ηdr′

Expanding the auxiliary density and density gradient derivatives according to Eqs. (7) and (8) results into the following:18 ∂vxcσ[ρ~,γ~]∂η=∑τα,β∑l¯xl¯τ(η)[∫∂vxcσ(ρ~,γ~)∂ρ~τ(r′)l¯(r′)dr′+∑vx,y,z∫∂vxcσ(ρ~,γ~)∂ρ~vτ(r′)∂l¯(r′)∂v′dr′]+∑τα,β∑l¯xl¯τ[∫∂vxcσ(ρ~,γ~)∂ρ~τ(r′)l¯(η)(r′)dr′+∑vx,y,z∫∂vxcσ(ρ~,γ~)∂ρ~vτ(r′)∂l¯(η)(r′)∂v′dr′]

Inserting the explicit expression of the GGA exchange-correlation potential from Eq. (10) into Eq. (18) yields the following:19 ∂vxcσ[ρ~,γ~]∂η=∑τα,β∑l¯xl¯τ(η)[∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~σ(r)l¯(r′)dr′+∑ux,y,z∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~uσ(r)l¯(r′)dr′∂∂u]+∑τα,β∑l¯xl¯τ(η)[∑vx,y,z∫∂2exc(ρ~,γ~)∂ρ~vτ(r′)∂ρ~σ(r)∂l¯(r′)∂v′dr′+∑u,vx,y,z∫∂2exc(ρ~,γ~)∂ρ~vτ(r′)∂ρ~uσ(r)∂l¯(r′)∂v′dr′∂∂u]+∑τα,β∑l¯xl¯τ[∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~σ(r)l¯(η)(r′)dr′+∑ux,y,z∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~uσ(r)l¯(η)(r′)dr′∂∂u]+∑τα,β∑l¯xl¯τ[∑vx,y,z∫∂2exc(ρ~,γ~)∂ρ~vτ(r′)∂ρ~σ(r)∂l¯(η)(r′)∂v′dr′+∑u,vx,y,z∫∂2exc(ρ~,γ~)∂ρ~vτ(r′)∂ρ~uσ(r)∂l¯(η)(r′)∂v′dr′∂∂u]

At this point, it is convenient to define the GGA ADFT exchange-correlation kernel according to Eq. (19). Therefore, we obtain as GGA kernel expression in ADFT:20 fxcστ(ρ~,γ~)≡∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~σ(r)+∑ux,y,z∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~uσ(r)∂∂u+∑vx,y,z∂2exc(ρ~,γ~)∂ρ~vτ(r′)∂ρ~σ(r)∂∂v′+∑u,vx,y,z∂2exc(ρ~,γ~)∂ρ~vτ(r′)∂ρ~uσ(r)∂∂v′∂∂u

As Eq. (20) shows, the exchange-correlation kernel has in general four spin contributions, namely αα,αβ,βα, and ββ. As developed in Eq. (20), the kernel is also non-local, i.e., depends from r and r′. However, for LDA and GGA, exc(ρ~,γ~) is an ordinary function. As a result, the variables r and r′ collapse in the integration over the exchange-correlation kernel as outlined in [38] and [44].

Back-substitution of the kernel definition, Eq. (20), into the exchange-correlation potential derivative, Eq. (19), yields for the matrix element in the last term of Eq. (16):21 k¯(λ)|∂vxcσ[ρ~,γ~]∂η=∑τα,β∑l¯xl¯τ(η)∬fxcστ(ρ~,γ~)l¯(r′)k¯(λ)(r)δ(r-r′)drdr′+∑τα,β∑l¯xl¯τ∬fxcστ(ρ~,γ~)l¯(η)(r′)k¯(λ)(r)δ(r-r′)drdr′=∑τα,β∑l¯xl¯τ(η)k¯(λ)|fxcστ[ρ~,γ~]|l¯+∑τα,β∑l¯xl¯τk¯(λ)|fxcστ[ρ~,γ~]|l¯(η)

Therefore, we find explicit contributions to the ADFT Hessian matrix elements from the second derivatives of the exchange-correlation energy functional with respect to the atomic coordinates η and λ:22 ∂2Exc[ρ~,γ~]∂η∂λ|explicit=∑σα,β∑k¯xk¯σ(η)k¯(λ)|vxcσ[ρ~,γ~]+∑σα,β∑k¯xk¯σk¯(ηλ)|vxcσ[ρ~,γ~]+∑σ,τα,β∑k¯,l¯xk¯σk¯(λ)|fxcστ[ρ~,γ~]|l¯xl¯τ(η)+∑σ,τα,β∑k¯,l¯xk¯σk¯(λ)|fxcστ[ρ~,γ~]|l¯(η)xl¯τ

For the kernel implementation in deMon2k, it is more convenient to formulate the exchange-correlation kernel of Eq. (20) in terms of ρ~σ(r) and γ~στ(r). Applying the chain rule to the terms involving the σ and τ density gradients in Eq. (19) and using Eq. (13), we obtain the following:23 ∂vxcσ[ρ~,γ~]∂η=∑τα,β∑l¯xl¯τ(η)∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~σ(r)l¯(r′)dr′+2∑σ′,τα,β∑l¯∑ux,y,zxl¯τ(η)∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂γ~σσ′(r)l¯(r′)dr′ρ~u′(r)∂∂u+2∑τ,τ′α,β∑l¯∑vx,y,zxl¯τ(η)∫∂2exc(ρ~,γ~)∂γ~ττ′(r′)∂ρ~σ(r)ρ~v′(r′)∂l¯(r′)∂v′dr′+4∑σ′,τ,τ′α,β∑l¯∑u,vx,y,zxl¯τ(η)∫∂2exc(ρ~,γ~)∂γ~ττ′(r′)∂γ~σσ′(r)ρ~v′(r′)∂l¯(r′)∂v′dr′ρ~u′(r)∂∂u+2∑τα,β∑l¯∑vx,y,zxl¯τ(η)∫∂exc(ρ~,γ~)∂γ~στ(r)∂l¯(r′)∂v′dr′∂∂v+∑τα,β∑l¯xl¯τ∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~σ(r)l¯(η)(r′)dr′+2∑σ′,τα,β∑l¯∑ux,y,zxl¯τ∫∂2exc(ρ~,γ~)∂ρ~τ(r′)∂γ~σσ′(r)l¯(η)(r′)dr′ρ~u′(r)∂∂u+2∑τ,τ′α,β∑l¯∑vx,y,zxl¯τ∫∂2exc(ρ~,γ~)∂γ~ττ′(r′)∂ρ~σ(r)ρ~v′(r′)∂l¯(η)(r′)∂v′dr′+4∑σ′,τ,τ′α,β∑l¯∑u,vx,y,zxl¯τ∫∂2exc(ρ~,γ~)∂γ~ττ′(r′)∂γ~σσ′(r)ρ~v′(r′)∂l¯(η)(r′)∂v′dr′ρ~u′(r)∂∂u+2∑τα,β∑l¯∑vx,y,zxl¯τ∫∂exc(ρ~,γ~)∂γ~στ(r)∂l¯(η)(r′)∂v′dr′∂∂v

Hence, the implemented GGA kernel formulation in ADFT is given by the following:24 fxcστ(ρ~,γ~)≡∂2exc(ρ~,γ~)∂ρ~τ(r′)∂ρ~σ(r)+2∑σ′α,β∑ux,y,z∂2exc(ρ~,γ~)∂ρ~τ(r′)∂γ~σσ′(r)ρ~u′(r)∂∂u+2∑τ′α,β∑vx,y,z∂2exc(ρ~,γ~)∂γ~ττ′(r′)∂ρ~σ(r)ρ~vτ(r′)∂∂v′+4∑σ′,τ′α,β∑u,vx,y,z∂2exc(ρ~,γ~)∂γ~ττ′(r′)∂γ~σσ′(r)ρ~v′(r′)ρ~u′(r)×∂∂v′∂∂u+2∑vx,y,z∂exc(ρ~,γ~)∂γ~στ(r)∂∂v′∂∂v

Finite difference GGA exchange-correlation kernel

In the last section, the working formula for the unrestricted GGA ADFT kernel, Eq. (24), in deMon2k was presented. This implementation requires second derivatives of the density-weighted exchange-correlation energy density. In cases where these derivatives are not available, the matrix elements of the GGA ADFT kernel can be calculated by a finite difference formulation. This approach requires only the implementation of the corresponding GGA exchange-correlation potential formulas that are mandatory for SCF calculations. To this end, the exchange-correlation potential is evaluated at symmetric arbitrarily small changes in the auxiliary density at a given grid point. Following the definition of derivatives, the unrestricted finite difference formulation of the GGA ADFT kernel can be expressed as follows:25 k¯(λ)|fxcστ|=∫δvxcσ[ρ~τ(r)]δρ~τ(r′)k¯(λ)(r′)dr′≈vxcσ[ρ~τ(r)+ϵk¯(λ)(r)]-vxcσ[ρ~τ(r)-ϵk¯(λ)(r)]2ϵ

Equation (25) represents the arbitrary changes in the auxiliary density by ±ϵk¯(λ)(r), where ϵ is the step size, and k¯(λ)(r) is the derivative of the auxiliary function k¯ with respect to the atomic coordinate λ. This choice is made with the intention of constructing the kernel integrals appearing in Eq. (22).

Take as an example the kernel integral in the third term of Eq. (22). It can be expressed by the following finite difference formulation:26 k¯(λ)|fxcτσ|l¯≈∫vxcσ[ρ~τ(r)+ϵk¯(λ)(r)]-vxcσ[ρ~τ(r)-ϵk¯(λ)(r)]2ϵl¯(r)dr

Equation (26) can be used to calculate the exchange-correlation kernel integrals appearing in Eq. (22) for any GGA functional for which exchange-correlation potential formulas are implemented. The integral in Eq. (26) is evaluated numerically with the same grid used in the SCF. In deMon2k, the step size ϵ is set to 10-9 a.u., which yields consistent and robust results [20].

Analytic second derivatives of ECPs and MCPs

In the original frequency analysis module of deMon2k, the second derivatives of ECPs and MCPs are calculated by finite differences from the corresponding analytic gradients. For systems with many ECP or MCP centers, these finite difference calculations can become a computational bottleneck. Therefore, we implemented within this work analytic second derivatives for ECPs and MCPs. The relevant contributions to the Hessian matrix from second derivatives of ECPs or MCPs, U, centered on atom C have the general form:27 HCλη=∑μνPμν∂2∂λ∂η⟨μ|U|ν⟩

Taking rC=r-C for ECPs, the pseudopotential operator is given by the following:28 UECP(rC)=UL(rC)+∑l=0L-1∑m=-ll|Slm(r^C)⟩Ul(rC)⟨Slm(r^C)|.

The first term on the right-hand side is a purely local central force potential located on atom C. The second term is a double sum of semi-local operators with a fully local radial part and non-local angular part, here expressed by normalized real spherical harmonics, Slm(r^C). In this notation, r^C is the unit vector in the direction of rC. For MCPs, a similar operator form is used:29 UMCP(rC)=UL(rC)+∑l=0L-1∑m=-ll|ϕlm(rC)⟩Bl⟨ϕlm(rC)|.

In Eq. (29), Bl is a constant and ϕlm(rC) is a core orbital. Since core orbitals are normalized, the terms inside the sum are true projectors, fully non-local. Analytical calculation of second derivatives was implemented taking advantage of translational invariance. Differentiation of the operators, especially the semi-local ECP operator, is not advisable. Therefore, explicit differentiation of the pseudopotential operators was avoided, only basis functions are differentiated.30 HCλη=∑μνPμν∂∂η[(δλA-δλC)⟨μ(λ)|U|ν⟩+(δλB-δλC)⟨μ|U|ν(λ)⟩]=2∑μνPμν[(δλA-δλC)(δηA-δηC)⟨μ(λη)|U|ν⟩+(δλA-δλC)(δηB-δηC)⟨μ(λ)|U|ν(η)⟩]

The Kronecker delta in Eq. (30) tests only if λ or η correspond to the same atom where the differentiated atomic basis function is located. The derivative of a Cartesian Gaussian function is a combination of Gaussians with shifted angularity: ± 1 for first derivatives and ± 2 and 0 for second derivatives. Avoiding to differentiate the pseudopotential operator allows us to apply directly the semi-numerical integration algorithm developed for ECPs [45]. Evaluation of semi-local ECPs and local ECPs and MCPs’ second derivatives has a computational cost slightly larger than the evaluation of the corresponding energy integrals. The extra cost arises from the increased angularity of the required angular integrals which are, however, independent of the number of atoms. Non-local projectors in MCPs are easily decomposed into core-valence overlap integrals, and their derivatives can be evaluated very efficiently since they are completely analog to the overlap matrix derivatives.Table 1 Analytic unrestricted ADFT harmonic frequencies (cm–1) of H2O calculated with the LDA VWN and GGA PBE functionals employing their native deMon2k and LibXC implementations. In the analytic frequency analyses, analytic kernel (ANK) and finite difference kernel (FDK) calculations are used. The H2O structure for these frequency analyses is depicted on the right. See text for further details

		

Computational details

The open-shell analytical GGA exchange-correlation kernel formula, as detailed in the preceding section, Eq. (24), has been incorporated into the LCGTO-ADFT code deMon2k [46]. For the kernel calculations, either the native or LibXC partial derivatives of exc(ρ~,γ~) are used. We denote them by native or LibXC in the following. The test calculations are performed with the Vosko, Wilk, and Nusair (VWN) LDA functional [47] in conjunction with Dirac exchange [48] or the GGA Perdew-Burke-Ernzerhof (PBE) functional [49] as implemented in deMon2k or LibXC. The Kohn-Sham orbitals are expanded with the double-zeta valence polarized (DZVP) basis set [50], while the auxiliary density is computed using the automatically generated GEN-A2 auxiliary function set [51]. Thus, the variational fitting of the Coulomb potential [5] is used, and the exchange-correlation energies, potentials, and kernels are calculated with the resulting auxiliary density. This ADFT approach is also used for the analytic frequency calculations. In all calculations, the self-consistent field (SCF) energy and auxiliary density convergence criteria are set to 10-5 a.u. and 5×10-4 a.u., respectively. For the numerical integration of the exchange-correlation energy and its derivatives, an adaptive grid with a grid tolerance of 10-5 a.u. was used.

As benchmark calculations, we performed harmonic frequency analyses of (Pt3Cu)n clusters with n being 1 to 11. The calculations are executed with the RPBE [52] GGA functional, again using either the native or LibXC implementation of it. In these calculations, the extended GEN-A2* auxiliary function sets [53] are used. For the platinum atoms, a double-zeta basis set from the Los Alamos National Laboratory was employed in conjunction with an 18-electron quasi-relativistic effective core potential (QECP| LANL2DZ) [54]. For the copper atoms, an all-electron triple-zeta valence polarized (TZVP) basis set optimized for GGA functionals [51] was used. Initial structures of the (Pt3Cu)n clusters were taken from [11, 12] and tightly optimized using an SCF energy convergence criteria and an auxiliary density convergence criteria of 5×10-8 a.u. and 10-5 a.u., respectively. The structure optimization tolerance was set to 3×10-4 a.u. for the root-mean-square forces. The frequencies are calculated using the optimized cluster structures by analytic second ADFT energy derivatives employing ADPT. The required exchange-correlation kernels were calculated with the analytical and finite difference methods. In these calculations, a fixed fine grid was used for the numerical integration of the exchange-correlation energy, potential, and kernel. All calculations were performed using a 24-processor parallel architecture.Table 2 Analytic unrestricted ADFT harmonic frequencies (cm–1) of small open-shell systems calculated with the LDA VWN and GGA PBE functionals employing their native deMon2k and LibXC implementations. In the analytic frequency analyses, analytic kernel (ANK) and finite difference kernel (FDK) calculations are used. The molecular structures are taken from [55]. See text for further details

Molecule	VWN	PBE	
	Native	LibXC	Native	LibXC	
	ANK	FDK	ANK	FDK	ANK	FDK	ANK	FDK	
2BeH	1666.5	1666.0	1671.2	1672.9	1658.6	1659.0	1665.4	1664.8	
2CH3 ν1	444.8	444.8	444.9	444.9	287.2	287.1	287.1	287.1	
2CH3 ν2	1321.4	1321.3	1321.3	1321.1	1373.2	1373.0	1373.3	1373.0	
2CH3 ν3	1321.4	1321.3	1321.3	1321.1	1373.2	1373.0	1373.3	1373.0	
2CH3 ν4	3061.5	3061.5	3061.5	3061.5	3033.2	3033.2	3033.2	3033.2	
2CH3 ν5	3237.1	3237.1	3237.2	3237.1	3201.3	3201.2	3201.3	3201.2	
2CH3 ν6	3237.1	3237.1	3237.2	3237.1	3201.3	3201.2	3201.3	3201.2	
3Li2	332.8	332.8	332.9	333.0	337.9	337.9	338.2	337.9	
3Na2	162.1	162.1	162.1	162.2	166.3	166.2	166.2	166.2	
2NH2 ν1	1452.0	1452.0	1451.9	1451.9	1493.4	1493.4	1493.4	1493.4	
2NH2 ν2	3304.1	3304.1	3304.0	3304.0	3263.7	3263.7	3263.7	3263.7	
2NH2 ν3	3407.5	3407.5	3407.3	3407.3	3369.0	3369.0	3368.9	3368.9	
3NH	3190.0	3190.0	3189.7	3189.7	3149.2	3149.2	3149.2	3149.2	
2NO	1990.2	1990.0	1990.2	1989.9	1997.5	1997.5	1997.5	1997.5	
3O2	1524.4	1524.3	1524.4	1524.2	1536.5	1536.5	1536.5	1537.0	
2PH2 ν1	1033.9	1029.0	1033.9	1029.6	1052.8	1052.5	1052.8	1052.6	
2PH2 ν2	2321.9	2320.1	2321.7	2320.8	2346.5	2346.2	2346.4	2346.3	
2PH2 ν3	2337.7	2334.3	2338.0	2335.6	2358.0	2357.7	2357.9	2357.9	
2SiF	770.0	769.7	769.8	769.6	782.2	782.2	782.2	782.2	
MAD	Ref.	0.6	0.4	0.8	Ref.	0.1	0.4	0.4	

Validation calculations

To validate the GGA ADFT kernel working equation, Eq. (24), we performed frequency analyses of the triplet H2O molecule at a not optimized geometry. To do so, we use the singlet VWN/DZVP/GEN-A2 optimized H2O structure, as depicted in Table 1, for our triplet frequency analyses. The obtained harmonic frequencies are listed in Table 1, too. Here, we compare analytic ADFT frequencies obtained with the LDA VWN functional and the GGA PBE functional. For both functionals, the native deMon2k functional implementation is compared with the corresponding LibXC implementation. For both implementations analytic (ANK) and finite difference (FDK), kernel results are listed. Each column displays the three harmonic frequencies of the triplet H2O at the optimized singlet geometry. As can be seen from Table 1, the first frequency is always imaginary (indicated by a minus sign). For the LDA VWN frequencies, excellent agreement between all methods is observed. Deviations are in the range of 1 cm–1 or below. The situation changes for the GGA PBE frequencies. Whereas the agreement of harmonic frequencies calculated with the analytic and finite difference kernels remains excellent, significant deviations between the native and LibXC kernel calculations are found. These differences are largest, up to 20 cm–1, for the imaginary frequency and reduce with increasing frequencies as Table 1 shows. We attribute this to different screenings in the native and LibXC partial derivatives of the density-weighted PBE exchange-correlation energy density. Similar deviations are found for other GGA functionals, too.

To further investigate the differences between the native and LibXC implementations of exchange-correlation functionals, we performed frequency analyses on small open-shell systems. In these unrestricted calculations, the VWN/DZVP/GEN-A2 and the PBE/DZVP/GEN-A2 level of theories were used. To be unbiased with respect to optimized structure parameters, we used structure parameters from the National Institute of Standards and Technology’s Computational Chemistry Comparison and Benchmark Database (NIST) [55] for the frequency analyses.

Table 2 compares the calculated harmonic frequencies of ten selected open-shell systems with each other. The table is structured in the same way as Table 1 adding mean absolute deviations (MADs) with respect to the native deMon2k ANK results at the bottom of the table. As Table 2 shows, a general good to excellent agreement between all frequencies is observed. This is confirmed by the MADs that are always below 1 cm–1. In particular, the difference between harmonic frequencies from the native and LibXC functional implementations is reduced well below 10 cm–1. Largest deviations are observed for 2BeH, both for the LDA VWN (∼ 5 cm–1) and GGA PBE (∼ 7 cm–1) functionals. Furthermore, Table 2 shows that the differences between native and LibXC implementations are now in the same size range as for analytic (ANK) and finite difference (FDK) kernel calculations (∼ 5 cm–1 for 2PH2 ν1 with the LDA VWN functional). This supports our previous assumption that these differences arise from different screenings in the functional implementations. We attribute the improvement of the consistency of the frequencies in Table 2 with respect to Table 1 to the molecular geometries. Whereas the molecular structure parameters for the systems in Table 2 are close to optimized structure data, the triplet H2O structure used in Table 1 is far from the corresponding minimum. This is also supported by the fact that the largest deviation between native and LibXC functional implementations is observed for the imaginary frequency of the triplet H2O example.Table 3 Computational timings (s) for finite difference (FD) and analytic (AN) ECP second derivative calculations of (Pt3Cu)n clusters, with n=1 to 11. The speed-up factor refers to frequency analyses performed on 24 cores

Cluster	FD	AN	Speed-up	
2Pt3Cu	2.2	0.2	11.5	
3Pt6Cu2	13.8	0.7	20.0	
4Pt9Cu3	43.4	1.5	28.8	
3Pt12Cu4	102.8	2.8	36.3	
4Pt15Cu5	182.8	4.1	44.4	
3Pt18Cu6	295.2	5.4	54.8	
2Pt21Cu7	504.7	8.2	61.8	
5Pt24Cu8	756.5	10.8	69.9	
2Pt27Cu9	1064.3	13.6	78.5	
3Pt30Cu10	1594.2	18.3	87.1	
6Pt33Cu11	2202.2	23.5	93.5	

Table 4 Analytic unrestricted ADFT harmonic frequencies (cm–1) of (Pt3Cu)n clusters, with n=1 to 11, calculated with the GGA RPBE functional using its native deMon2k and LibXC implementations. In the analytic frequency analyses, analytic kernel (ANK) and finite difference kernel (FDK) calculations are used. The molecular structures are re-optimizations of the structures reported in [11, 12]. See text for further details

Native	
Cluster	ANK	FDK	
	ν1	ν2	ν3	ν4	ν5	ν1	ν2	ν3	ν4	ν5	
2Pt3Cu	89.4	89.5	136.3	137.5	177.5	89.4	89.5	136.3	137.5	177.5	
3Pt6Cu2	14.2	36.6	41.6	61.6	75.3	14.2	36.6	41.6	61.6	75.3	
4Pt9Cu3	41.1	46.1	57.1	61.9	63.4	41.1	46.0	57.0	62.0	63.4	
3Pt12Cu4	16.1	37.5	37.9	47.1	47.8	16.2	37.5	37.9	47.1	47.8	
4Pt15Cu5	34.9	39.9	44.8	49.3	51.3	34.9	39.7	44.5	49.1	51.0	
3Pt18Cu6	23.7	36.0	40.4	42.7	47.8	23.3	36.2	40.3	42.8	47.9	
2Pt21Cu7	28.8	31.8	34.6	35.9	43.0	28.0	31.0	34.3	35.3	42.7	
5Pt24Cu8	23.9	29.9	34.5	35.9	39.4	24.4	30.6	34.2	36.4	39.8	
2Pt27Cu9	6.8	15.1	16.6	25.6	29.0	8.7	15.9	19.3	25.5	29.2	
3Pt30Cu10	28.0	28.8	34.6	35.9	37.9	28.1	29.3	34.3	36.2	37.1	
6Pt33Cu11	28.3	32.5	33.8	34.7	37.8	29.3	32.8	34.0	34.4	37.8	
MAD	Ref.	Ref.	Ref.	Ref.	Ref.	0.4	0.3	0.4	0.2	0.2	
LibXC	
Cluster	ANK	FDK	
	ν1	ν2	ν3	ν4	ν5	ν1	ν2	ν3	ν4	ν5	
2Pt3Cu	88.3	89.7	136.3	137.5	177.8	88.2	89.7	136.2	137.5	177.7	
3Pt6Cu2	18.1	36.2	44.4	61.9	75.1	19.1	36.4	44.5	61.4	74.9	
4Pt9Cu3	46.4	46.9	57.8	62.6	65.9	47.5	56.1	62.2	64.0	68.8	
3Pt12Cu4	14.1	37.3	40.4	46.6	48.7	12.9	36.5	39.9	45.8	48.2	
4Pt15Cu5	33.8	33.9	44.4	48.8	50.2	33.9	40.7	44.6	49.4	51.8	
3Pt18Cu6	22.7	36.3	40.7	43.2	47.7	22.2	36.1	40.1	43.1	47.6	
2Pt21Cu7	29.7	33.6	35.0	36.1	44.3	29.7	33.6	35.0	36.1	42.9	
5Pt24Cu8	27.1	29.7	32.4	36.7	39.6	27.0	30.6	32.5	36.5	39.0	
2Pt27Cu9	2.6	6.9	18.1	27.3	28.7	6.0	7.2	14.3	26.3	27.5	
3Pt30Cu10	27.9	28.6	34.2	36.2	38.2	27.9	28.6	34.2	36.2	38.2	
6Pt33Cu11	31.9	33.5	35.1	38.1	41.4	31.7	34.4	34.4	38.4	40.2	
MAD	2.4	1.8	1.1	0.8	1.0	2.4	2.3	1.5	0.9	1.3	

Benchmark calculations

To benchmark our new GGA ADFT kernel implementation, we revisited the (Pt3Cu)n clusters previously reported in [11, 12]. To this end, we re-optimized the cluster structures and performed frequency analyses with the optimized structures employing the native and LibXC implementations of the RPBE [52] functional. We used this functional because of its use in the original report. Due to the computational demand of these calculations, we profiled some of them. To our surprise, we found that a significant portion of the computational time was spent on the finite difference calculation of the ECPs. Table 3 lists these timings in column 2.

To overcome this computational bottleneck analytic second ECP (and MCP), derivatives were implemented in deMon2k as described in Section 2.4. The timings for the analytic ECP second derivative calculations are reported in column 3 of Table 3 and graphically compared with their finite difference counterparts in Fig. 1. Speed-ups of up to almost 100 are observed for the studied (Pt3Cu)n clusters! As Fig. 1 shows, the newly analytic implementation of ECP second derivatives exhibits a nearly linear scaling with respect to the number of ECP centers. Note that the semi-local integrator scaling is only due to the numerical radial integrator [45] because analytic angular integration is independent of the number of pseudopotential centers in the system.Fig. 1 Comparison of the computational timings for finite difference and analytic ECP second derivative calculation of (Pt3Cu)n clusters, with n=1 to 11. All calculations are performed on 24 cores

The five lowest harmonic frequencies of the optimized (Pt3Cu)n are listed in Table 4 for the native and LibXC functional implementations. The MADs with respect to the native deMon2k ANK frequencies are given at the end of the table. Again, the analytic (ANK) and finite difference (FDK) kernel calculations are compared for both implementations. Notably, there is an outstanding agreement between the analytic and finite difference kernel calculations as shown by the corresponding MADs that are all below 0.5 cm–1. Deviations rarely exceed three wavenumbers, underscoring the robustness of the deMon2k implementation. In contrast, the results obtained with the LibXC interface exhibit more deviations. While deviations are generally less than two wavenumbers (maximum MAD is 2.4 cm–1), some frequencies display discrepancies of up to 10 wavenumbers. Additionally, using the LibXC interface with the above-mentioned methodology, two systems revealed imaginary frequencies in the case of the FDK method, namely Pt6Cu2 and Pt27Cu9. To overcome this problem, we re-optimized these two clusters with a tighter structure optimization tolerance of 10-5 a.u. and an analytically calculated start Hessian matrix. In the subsequent frequency analyses, only positive frequencies are observed. For consistency, both ANK and FDK frequencies for Pt6Cu2 and Pt27Cu9 listed in Table 4 were obtained with this methodology. These results show that the here-described ADFT kernel implementations are numerically stable for delicate electronic structures typical of transition metal clusters.

Conclusions

The unrestricted GGA ADFT exchange-correlation kernel working equation for perturbation dependent auxiliary functions is derived and discussed. We validated and benchmarked this kernel implementation in the framework of frequency analyses using the native and LibXC exchange-correlation functional implementations in deMon2k. Overall, there is a very satisfying consistency between the native and LibXC frequencies. In both cases, deviations of less than a wavenumber are typically observed when comparing analytic and finite difference kernel calculations. Between the implementations, larger deviations of 5 cm–1 or more are occasionally observed. The largest deviation (∼ 20 cm–1) between native and LibXC frequencies is found for the imaginary mode in the non-minimum structure of 3H2O. We attribute this to the different screenings in the exchange-correlation energy density calculations. However, for (near) minimum structures, the agreement of native and LibXC frequencies is for small open-shell molecules usually excellent. In the case of (Pt3Cu)n clusters, results obtained using the LibXC interface show more significant deviations, albeit generally within 2 wavenumbers, with occasional discrepancies of up to 10 wavenumbers. In the case of Pt6Cu2 and Pt27Cu9, these deviations led to imaginary frequencies in the final optimized structures with the default convergence thresholds. By tightening these thresholds and using an analytically calculated start Hessian matrix, these problems were fixed and the clusters converged with the LibXC functional implementation to minima, too. The computational bottleneck due to the finite difference second derivative ECP (and MCP) integral calculations in the (Pt3Cu)n clusters was resolved by the implementation of corresponding analytic second derivatives. The resulting speed-ups are large and scale linearly with the number of ECP (MCP) centers. Consequently, frequency analyses of systems with hundreds of ECP (MCP) centers are now feasible with deMon2k in very reasonable times.

Acknowledgements

Flor Andrea Olvera-Rubalcava gratefully acknowledges CONAHCYT for her PhD fellowship 798159.

Author contribution

All authors contributed to the writing of the manuscript.

Funding

This work was supported by the CONAHCYT fellowship 798159 and the grants A1-S-11929 and ECOS 321168.

Data availability

No datasets were generated or analyzed during the current study.

Declarations

Conflict of interest

The authors declare no competing interests.

Dedicated to Professor Toro-Labbé on the occasion of his 70th birthday.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Hohenberg P Kohn W Inhomogeneous electron gas Phys Rev 1964 136 B864 B871 10.1103/PhysRev.136.B864
Hohenberg P, Kohn W (1964) Inhomogeneous electron gas. Phys Rev 136:B864–B87110.1103/PhysRev.136.B864
2. Kohn W Sham LJ Self-consistent equations including exchange and correlation effects Phys Rev 1965 140 A1133 A1138 10.1103/PhysRev.140.A1133
Kohn W, Sham LJ (1965) Self-consistent equations including exchange and correlation effects. Phys Rev 140:A1133–A113810.1103/PhysRev.140.A1133
3. Teale AM Helgaker T Savin A Adamo C Aradi B Arbuznikov AV Ayers P Baerends EJ Barone V Calaminici P Cances E Carter EA Chattaraj PK Chermette H Ciofini I Crawford TD De Proft F Dobson J Draxl C Frauenheim T Fromager E Fuentealba P Gagliardi L Galli G Gao J Geerlings P Gidopoulos N Gill PMW Gori-Giorgi P Görling A Gould T Grimme S Gritsenko O Jensen HJ Johnson ER Jones RO Kaupp M Köster AM Kronik L Krylov AI Kvaal S Laestadius A Levy MP Lewin M Liu SB Loos PF Maitra NT Neese F Perdew J Pernal K Pernot P Piecuch P Rebolini E Reining L Romaniello P Ruzsinszky A Salahub DR Scheffler M Schwerdtfeger P Staroverov VN Sun J Tellgren E Tozer DJ Trickey SB Ullrich CA Vela A Vignale G Wesolowski TA Xu X Yang W DFT exchange: sharing perspectives on the workhorse of quantum chemistry and materials science Phys Chem Chem Phys 2022 24 28700 28781 10.1039/D2CP02827A 36269074
Teale AM, Helgaker T, Savin A, Adamo C, Aradi B, Arbuznikov AV, Ayers P, Baerends EJ, Barone V, Calaminici P, Cances E, Carter EA, Chattaraj PK, Chermette H, Ciofini I, Crawford TD, De Proft F, Dobson J, Draxl C, Frauenheim T, Fromager E, Fuentealba P, Gagliardi L, Galli G, Gao J, Geerlings P, Gidopoulos N, Gill PMW, Gori-Giorgi P, Görling A, Gould T, Grimme S, Gritsenko O, Jensen HJ, Johnson ER, Jones RO, Kaupp M, Köster AM, Kronik L, Krylov AI, Kvaal S, Laestadius A, Levy MP, Lewin M, Liu SB, Loos PF, Maitra NT, Neese F, Perdew J, Pernal K, Pernot P, Piecuch P, Rebolini E, Reining L, Romaniello P, Ruzsinszky A, Salahub DR, Scheffler M, Schwerdtfeger P, Staroverov VN, Sun J, Tellgren E, Tozer DJ, Trickey SB, Ullrich CA, Vela A, Vignale G, Wesolowski TA, Xu X, Yang W (2022) DFT exchange: sharing perspectives on the workhorse of quantum chemistry and materials science. Phys Chem Chem Phys 24:28700–2878136269074 10.1039/D2CP02827A
4. Calaminici P Alvarez-Ibarra A Cruz-Olvera D Domínguez-Soria VD Flores-Moreno R Gamboa GU Geudtner G Goursot A Mejía-Rodríguez D Salahub DR Zuniga-Gutierrez B Köster AM Leszczynski J Kaczmarek-Kedziera A Puzyn T Papadopoulos GM Reis H Shukla KM Auxiliary density functional theory: from molecules to nanostructures Handbook of Computational Chemistry Part II: Applicaions of Computational Methods to Model Systems 2017 2 Switzerland Springer, Cham
Calaminici P, Alvarez-Ibarra A, Cruz-Olvera D, Domínguez-Soria VD, Flores-Moreno R, Gamboa GU, Geudtner G, Goursot A, Mejía-Rodríguez D, Salahub DR, Zuniga-Gutierrez B, Köster AM (2017) Auxiliary density functional theory: from molecules to nanostructures. In: Leszczynski J, Kaczmarek-Kedziera A, Puzyn T, Papadopoulos GM, Reis H, Shukla KM (eds) Handbook of Computational Chemistry Part II: Applicaions of Computational Methods to Model Systems, 2nd edn. Springer, Cham, Switzerland
5. Dunlap BI, Connolly JW, Sabin JR (1979) On first-row diatomic molecules and local density models. J Chem Phys 71:4993-4999
6. Dunlap BI Rösch N Trickey SB Variational fitting methods for electronic structure calculations Mol Phys 2010 108 3167 3180 10.1080/00268976.2010.518982
Dunlap BI, Rösch N, Trickey SB (2010) Variational fitting methods for electronic structure calculations. Mol Phys 108:3167-318010.1080/00268976.2010.518982
7. Polly R, Werner HJ, Manby FR, Knowles PJ (2004) Fast Hartree-Fock theory using local density fitting approximations. Mol Phys 102:2311–2321
8. Mejía-Rodríguez D Köster AM Robust and efficient variational fitting of Fock exchange J Chem Phys 2014 141 124114 10.1063/1.4896199 25273419
Mejía-Rodríguez D, Köster AM (2014) Robust and efficient variational fitting of Fock exchange. J Chem Phys 141:12411425273419 10.1063/1.4896199
9. Laikov DN Fast evaluation of density functional exchange-correlation terms using the expansion of the electron density in auxiliary basis sets Chem Phys Lett 1997 281 151 156 10.1016/S0009-2614(97)01206-2
Laikov DN (1997) Fast evaluation of density functional exchange-correlation terms using the expansion of the electron density in auxiliary basis sets. Chem Phys Lett 281:151–15610.1016/S0009-2614(97)01206-2
10. Köster AM Reveles JU del Campo JM Calculation of exchange-correlation potentials with auxiliary function densities J Chem Phys 2004 121 3417 3424 10.1063/1.1771638 15303904
Köster AM, Reveles JU, del Campo JM (2004) Calculation of exchange-correlation potentials with auxiliary function densities. J Chem Phys 121:3417–342415303904 10.1063/1.1771638
11. Galindo-Uribe CD, Calaminici P, Cruz-Martínez H, Cruz-Olvera D, Solorza-Feria O (2021) First-principle study of the structures, growth pattern, and properties of (PtCu), n= 1-9, clusters. J Chem Phys 154:154302
12. Galindo-Uribe CD Calaminici P Solorza-Feria O First-principle investigation of structures and energy properties of (Pt3Cu)n, n= 10?11 nanoclusters Theor Chem Acc 2023 142 23 10.1007/s00214-023-02963-4
Galindo-Uribe CD, Calaminici P, Solorza-Feria O (2023) First-principle investigation of structures and energy properties of (PtCu), n= 10?11 nanoclusters. Theor Chem Acc 142:2310.1007/s00214-023-02963-4
13. De HS Krishnamurty S Mishra D Pal S Finite temperature behavior of gas phase neutral Aun (3 ≤ n ≤ 10) clusters: a first principles investigation J Phys Chem C 2011 115 17278 17285 10.1021/jp2023605
De HS, Krishnamurty S, Mishra D, Pal S (2011) Finite temperature behavior of gas phase neutral Au (3 n 10) clusters: a first principles investigation. J Phys Chem C 115:17278–1728510.1021/jp2023605
14. Luna-Valenzuela A Pedroza-Montero JN Köster AM Calaminici P Gálvez-González LE Posada-Amarillas A Pd8 cluster: too small to melt? A BOMD study J Phys Chem A 2024 128 572 580 10.1021/acs.jpca.3c06173 38207112
Luna-Valenzuela A, Pedroza-Montero JN, Köster AM, Calaminici P, Gálvez-González LE, Posada-Amarillas A (2024) Pd cluster: too small to melt? A BOMD study. J Phys Chem A 128:572–58038207112 10.1021/acs.jpca.3c06173
15. Epstein ST (1974) The variation method in quantum chemistry. Academic Press, London
16. Bérces A, Dickson RM, Fan L, Jacobsen H, Swerhone D, Ziegler T (1997) An implementation of the coupled perturbed Kohn-Sham equations: perturbation due to nuclear displacements. Comput Phys Commun 100:247–262
17. Flores-Moreno R Köster AM Auxiliary density perturbation theory J Chem Phys 2008 128 134105 10.1063/1.2842103 18397051
Flores-Moreno R, Köster AM (2008) Auxiliary density perturbation theory. J Chem Phys 128:13410518397051 10.1063/1.2842103
18. Mejía-Rodríguez D Delgado Venegas RI Calaminici P Köster AM Robust and efficient auxiliary density perturbation theory calculations J Chem Theory Comput 2015 11 1493 1500 10.1021/ct501065g 26574360
Mejía-Rodríguez D, Delgado Venegas RI, Calaminici P, Köster AM (2015) Robust and efficient auxiliary density perturbation theory calculations. J Chem Theory Comput 11:1493–150026574360 10.1021/ct501065g
19. Carmona-Espíndola J Flores-Moreno R Köster AM Time-dependent auxiliary density perturbation theory J Chem Phys 2010 133 084102 10.1063/1.3478551 20815555
Carmona-Espíndola J, Flores-Moreno R, Köster AM (2010) Time-dependent auxiliary density perturbation theory. J Chem Phys 133:08410220815555 10.1063/1.3478551
20. Shedge SV Carmona-Espíndola J Pal S Köster AM Comparison of the auxiliary density perturbation theory and the noniterative approximation to the coupled perturbed Kohn-Sham method: case study of the polarizabilities of disubstituted azoarene molecules J Phys Chem A 2010 114 2357 2364 10.1021/jp909966f 20088563
Shedge SV, Carmona-Espíndola J, Pal S, Köster AM (2010) Comparison of the auxiliary density perturbation theory and the noniterative approximation to the coupled perturbed Kohn-Sham method: case study of the polarizabilities of disubstituted azoarene molecules. J Phys Chem A 114:2357–236420088563 10.1021/jp909966f
21. Carmona-Espíndola J, Flores-Moreno R, Köster AM (2012) Static and dynamic first hyperpolarizabilities from time-dependent auxiliary density perturbation theory. Int J Quantum Chem 112:3461–3471
22. Calaminici P, Carmona-Espíndola J, Geudtner G, Köster AM (2012) Static and dynamic polarizability of C fullerene. Int J Quantum Chem 112:3252–3255
23. Karne AS Vaval N Pal S Vásquez-Pérez JM Köster AM Calaminici P Systematic comparison of DFT and CCSD dipole moments, polarizabilities and hyperpolarizabilities Chem Phys Lett 2015 635 168 173 10.1016/j.cplett.2015.06.046
Karne AS, Vaval N, Pal S, Vásquez-Pérez JM, Köster AM, Calaminici P (2015) Systematic comparison of DFT and CCSD dipole moments, polarizabilities and hyperpolarizabilities. Chem Phys Lett 635:168–17310.1016/j.cplett.2015.06.046
24. Flores-Moreno R Melin J Ortiz JV Merino G Efficient evaluation of analytic Fukui functions J Chem Phys 2008 129 224105 10.1063/1.3036926 19071905
Flores-Moreno R, Melin J, Ortiz JV, Merino G (2008) Efficient evaluation of analytic Fukui functions. J Chem Phys 129:22410519071905 10.1063/1.3036926
25. Flores-Moreno R Symmetry conservation in Fukui functions J Chem Theory Comput 2010 6 48 54 10.1021/ct9002527 26614318
Flores-Moreno R (2010) Symmetry conservation in Fukui functions. J Chem Theory Comput 6:48–5426614318 10.1021/ct9002527
26. Flores-Ramos JA Valdez-Ruvalcaba J González-Ochoa HO Flores-Moreno R Electron binding energies from static linear response calculations Theor Chem Acc 2021 140 131 10.1007/s00214-021-02831-z
Flores-Ramos JA, Valdez-Ruvalcaba J, González-Ochoa HO, Flores-Moreno R (2021) Electron binding energies from static linear response calculations. Theor Chem Acc 140:13110.1007/s00214-021-02831-z
27. Flores-Moreno R Cortes-Llamas SA Pineda-Urbina K Medel-Juarez VM Jayaprakash GK Analytic alchemical derivatives for the analysis of differential acidity assisted by the h function J Phys Chem A 2021 125 10463 10474 10.1021/acs.jpca.1c07364 34812636
Flores-Moreno R, Cortes-Llamas SA, Pineda-Urbina K, Medel-Juarez VM, Jayaprakash GK (2021) Analytic alchemical derivatives for the analysis of differential acidity assisted by the h function. J Phys Chem A 125:10463–1047434812636 10.1021/acs.jpca.1c07364
28. Zuniga-Gutierrez B, Geudtner G, Köster AM (2011) NMR shielding tensors from auxiliary density functional theory. J Chem Phys 134:124108
29. Zuniga-Gutierrez B Geudtner G Köster AM Magnetizability tensors from auxiliary density functional theory J Chem Phys 2012 137 094113 10.1063/1.4749243 22957561
Zuniga-Gutierrez B, Geudtner G, Köster AM (2012) Magnetizability tensors from auxiliary density functional theory. J Chem Phys 137:09411322957561 10.1063/1.4749243
30. Zuniga-Gutierrez B Camacho-Gonzalez M Simon-Bastida P Bendana-Castillo A Calaminici P Köster AM Efficient calculation of nuclear spin-rotation constants from auxiliary density functional theory J Chem Phys 2015 143 104103 10.1063/1.4929999 26374014
Zuniga-Gutierrez B, Camacho-Gonzalez M, Simon-Bastida P, Bendana-Castillo A, Calaminici P, Köster AM (2015) Efficient calculation of nuclear spin-rotation constants from auxiliary density functional theory. J Chem Phys 143:10410326374014 10.1063/1.4929999
31. Zuniga-Gutierrez B Camacho-Gonzalez M Simon-Bastida P Bendana-Castillo A Calaminici P Köster AM Efficient calculation of the rotational g tensor from auxiliary density functional theory J Phys Chem A 2015 119 1469 1477 10.1021/jp505169k 24968112
Zuniga-Gutierrez B, Camacho-Gonzalez M, Simon-Bastida P, Bendana-Castillo A, Calaminici P, Köster AM (2015) Efficient calculation of the rotational g tensor from auxiliary density functional theory. J Phys Chem A 119:1469–147724968112 10.1021/jp505169k
32. Zuniga-Gutierrez B Medel-Juarez V Varona A González Ramirez HN Flores-Moreno R Calculation of the EPR g-tensor from auxiliary density functional theory J Chem Phys 2020 152 014105 10.1063/1.5130174 31914741
Zuniga-Gutierrez B, Medel-Juarez V, Varona A, González Ramirez HN, Flores-Moreno R (2020) Calculation of the EPR g-tensor from auxiliary density functional theory. J Chem Phys 152:01410531914741 10.1063/1.5130174
33. López-Estrada O Selenius E Zuniga-Gutierrez B Malola S Häkkinen H Cubic aromaticity in ligand-stabilized doped Au superatoms J Chem Phys 2021 154 204303 10.1063/5.0050127 34241155
López-Estrada O, Selenius E, Zuniga-Gutierrez B, Malola S, Häkkinen H (2021) Cubic aromaticity in ligand-stabilized doped Au superatoms. J Chem Phys 154:20430334241155 10.1063/5.0050127
34. López-Estrada O Zuniga-Gutierrez B Selenius E Malola S Häkkinen H Magnetically induced currents and aromaticity in ligand-stabilized Au and AuPt superatoms Nat Commun 2021 12 2477 10.1038/s41467-021-22715-x 33931646
López-Estrada O, Zuniga-Gutierrez B, Selenius E, Malola S, Häkkinen H (2021) Magnetically induced currents and aromaticity in ligand-stabilized Au and AuPt superatoms. Nat Commun 12:247733931646 10.1038/s41467-021-22715-x
35. Ipatov A Fouqueau A Pérez del Valle C Cordova F Casida ME Köster AM Vela A Jamorski CJ Excitation energies from an auxiliary-function formulation of time-dependent density-functional response theory with charge conservation constraint J Mol Structure: THEOCHEM 2006 762 179 191 10.1016/j.theochem.2005.07.034
Ipatov A, Fouqueau A, Pérez del Valle C, Cordova F, Casida ME, Köster AM, Vela A, Jamorski CJ (2006) Excitation energies from an auxiliary-function formulation of time-dependent density-functional response theory with charge conservation constraint. J Mol Structure: THEOCHEM 762:179–19110.1016/j.theochem.2005.07.034
36. Carmona-Espíndola J Köster AM Photoabsorption spectra from time-dependent auxiliary density functional theory Can J Chem 2013 91 795 803 10.1139/cjc-2012-0501
Carmona-Espíndola J, Köster AM (2013) Photoabsorption spectra from time-dependent auxiliary density functional theory. Can J Chem 91:795–80310.1139/cjc-2012-0501
37. Hernández-Segura LI Köster AM Efficient implementation of time-dependent auxiliary density functional theory J Chem Phys 2023 158 024108 10.1063/5.0135263 36641386
Hernández-Segura LI, Köster AM (2023) Efficient implementation of time-dependent auxiliary density functional theory. J Chem Phys 158:02410836641386 10.1063/5.0135263
38. Zuniga-Gutierrez B, Köster AM (2016) Analytical GGA exchange-correlation kernel calculation in auxiliary density functional theory. Mol Phys 114:1026–1035
39. Delgado-Venegas RI Mejía-Rodríguez D Flores-Moreno R Calaminici P Köster AM Analytic second derivatives from auxiliary density perturbation theory J Chem Phys 2016 145 224103 10.1063/1.4971292 27984884
Delgado-Venegas RI, Mejía-Rodríguez D, Flores-Moreno R, Calaminici P, Köster AM (2016) Analytic second derivatives from auxiliary density perturbation theory. J Chem Phys 145:22410327984884 10.1063/1.4971292
40. Lehtola S, Steigemann C, Oliveira MJT, Marques MAL (2018) Recent developments in LIBXC: A comprehensive library of functionals for density functional theory. SoftwareX 7:1–5
41. Geudtner G, Calaminici P, Carmona-Espíndola J, del Campo JM, Domínguez-Soria VD, Flores-Moreno R, Gamboa GU, Goursot A, Köster AM, Reveles JU, Mineva T, Vásquez-Pérez JM, Vela A, Zuniga-Gutierrez B, Salahub DR (2012) deMon2k. WIREs Comput Mol Sci 2:548–555
42. Pulay P Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules: I. Theory Mol Phys 1969 17 197 204 10.1080/00268976900100941
Pulay P (1969) Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules: I. Theory. Mol Phys 17:197–20410.1080/00268976900100941
43. Baker J Andzelm J Scheiner A Delley B The effect of grid quality and weight derivatives in density functional calculations J Chem Phys 1994 101 8894 8902 10.1063/1.468081
Baker J, Andzelm J, Scheiner A, Delley B (1994) The effect of grid quality and weight derivatives in density functional calculations. J Chem Phys 101:8894–890210.1063/1.468081
44. Parr RG, Yang W (1989) Density-functional theory of atoms and molecules. Oxford University Press, New York
45. Flores-Moreno R Alvarez-Mendez RJ Vela A Köster AM Half-numerical evaluation of pseudopotential integrals J Comput Chem 2006 27 1009 1019 10.1002/jcc.20410 16628539
Flores-Moreno R, Alvarez-Mendez RJ, Vela A, Köster AM (2006) Half-numerical evaluation of pseudopotential integrals. J Comput Chem 27:1009–101916628539 10.1002/jcc.20410
46. Köster AM Geudtner G Alvarez-Ibarra A Calaminici P Casida ME Carmona-Espindola J Dominguez VD Flores-Moreno R Gamboa GU Goursot A Heine T Ipatov A de la Lande A Janetzko F del Campo JM Mejia-Rodriguez D Reveles JU Vasquez-Perez J Vela A Zuniga-Gutierrez B Salahub DR deMon2k, Version 6 2018 Cinvestav, Mexico City The deMon developers
Köster AM, Geudtner G, Alvarez-Ibarra A, Calaminici P, Casida ME, Carmona-Espindola J, Dominguez VD, Flores-Moreno R, Gamboa GU, Goursot A, Heine T, Ipatov A, de la Lande A, Janetzko F, del Campo JM, Mejia-Rodriguez D, Reveles JU, Vasquez-Perez J, Vela A, Zuniga-Gutierrez B, Salahub DR (2018) deMon2k, Version 6. The deMon developers, Cinvestav, Mexico City
47. Vosko SH Wilk L Nusair M Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis Can J Phys 1980 58 1200 1211 10.1139/p80-159
Vosko SH, Wilk L, Nusair M (1980) Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis. Can J Phys 58:1200–121110.1139/p80-159
48. Dirac PAM Note on exchange phenomena in the Thomas atom Math Proc Cambridge 1930 26 376 385 10.1017/S0305004100016108
Dirac PAM (1930) Note on exchange phenomena in the Thomas atom. Math Proc Cambridge 26:376–38510.1017/S0305004100016108
49. Perdew JP Burke K Ernzerhof M Generalized gradient approximation made simple Phys Rev Lett 1996 77 3865 3868 10.1103/PhysRevLett.77.3865 10062328
Perdew JP, Burke K, Ernzerhof M (1996) Generalized gradient approximation made simple. Phys Rev Lett 77:3865–386810062328 10.1103/PhysRevLett.77.3865
50. Godbout N Salahub DR Andzelm J Wimmer E Optimization of Gaussian-type basis sets for local spin density functional calculations. Part I. Boron through neon, optimization technique and validation Can J Phys 1992 70 560 571
Godbout N, Salahub DR, Andzelm J, Wimmer E (1992) Optimization of Gaussian-type basis sets for local spin density functional calculations. Part I. Boron through neon, optimization technique and validation. Can J Phys 70:560–571
51. Calaminici P Janetzko F Köster AM Mejia-Olvera R Zuniga-Gutierrez B Density functional theory optimized basis sets for gradient corrected functionals: 3d transition metal systems J Chem Phys 2007 126 044108 10.1063/1.2431643 17286463
Calaminici P, Janetzko F, Köster AM, Mejia-Olvera R, Zuniga-Gutierrez B (2007) Density functional theory optimized basis sets for gradient corrected functionals: transition metal systems. J Chem Phys 126:04410817286463 10.1063/1.2431643
52. Zhang Y, Yang W (1999) Comment on “Generalized gradient approximation made simple”. Phys Rev Lett 80:890
53. Calaminici P Flores-Moreno R Köster AM A density functional study of structures and vibrations of Ta3O and Ta3O- Comput Lett 2005 1 164 171 10.1163/157404005776611420
Calaminici P, Flores-Moreno R, Köster AM (2005) A density functional study of structures and vibrations of TaO and TaO. Comput Lett 1:164–17110.1163/157404005776611420
54. Hay PJ Wadt WR Ab initio effective core potentials for molecular calculations. Potentials for K to Au including the Outermost Core Orbitals J Chem Phys 1985 82 270 283 10.1063/1.448799
Hay PJ, Wadt WR (1985) Ab initio effective core potentials for molecular calculations. Potentials for K to Au including the Outermost Core Orbitals. J Chem Phys 82:270–28310.1063/1.448799
55. National Institute of Standards and Technology, Computational Chemistry Comparison and Benchmark DataBase. https://cccbdb.nist.gov/introx.asp. Accessed September 2023
