A new framework for frequency-dependent polarizable force fields

A frequency-dependent extension of the polarizable force field ``Atom-Condensed Kohn-Sham density functional theory approximated to the second-order'' (ACKS2) [J. Chem. Phys. 141, 194114 (2014)] is proposed, referred to as ACKS2$\omega$. The method enables theoretical predictions of dynamical response properties of finite systems after a partitioning of the frequency-dependent molecular response function. Parameters in this model are computed simply as expectation values of an electronic wavefunction, and the hardness matrix is entirely reused from ACKS2 as an adiabatic approximation is used. A numerical validation shows that accurate models can already be obtained with atomic monopoles and dipoles. Absorption spectra of 42 organic and inorganic molecular monomers are evaluated using ACKS2$\omega$, and our results agree well with the time-dependent DFT calculations. Also for the calculation of $C_6$ dispersion coefficients, ACKS2$\omega$ closely reproduces its TDDFT reference. When parameters for ACKS2$\omega$ are derived from a PBE/aug-cc-pVDZ ground state, it reproduces experimental values for 903 organic and inorganic intermolecular pairs with an MAPE of 3.84\%. Our results confirm that ACKS2$\omega$ offers a solid connection between the quantum-mechanical description of frequency-dependent response and computationally efficient force-field models.


I. INTRODUCTION
Molecular simulation is a powerful tool to investigate and predict the properties and dynamics of a wide range of finite and condensed-phase systems. Different levels of theory for such simulations can be roughly divided into two categories: quantum-mechanical (QM) electric structure methods and force-field (FF) methods. Kohn-Sham Density Functional Theory (DFT) offers a good trade-off between computational efficiency and accuracy, which has been successfully applied to study both ground and excited state properties of inorganic and organic molecular systems.
However, Kohn-Sham DFT (with density fitting) has an O(N 3 ) computational complexity where N is the number of basis functions, which impedes its application to large-scale systems, especially for molecular dynamics (MD) simulations. On the other hand, FF methods, only using a set of functions of molecular geometry with several empirical parameters to determine molecular energy, have been developed and applied to large-scale MD simulations. In general, FF methods have a much lower computational cost than QM methods, at the expense of a reduced, yet acceptable, computational accuracy. One significant problem in traditional FF methods is the lack of electronic polarization, i.e., the response of the electron distribution to the environment, which led the molecular modeling community to develop polarizable force fields (PFFs) [1][2][3][4][5][6][7][8][9][10][11][12] to improve the predictive accuracy of long-range interactions by adding explicit polarization effects.
Classical PFF models can be classified into four broad categories: induced dipoles or multipoles, [1][2][3][4][5] the Drude oscillator [6][7][8] (also known as "core-shell" 13 or "charge-on-spring" 14 ), Electronegativity Equalization Method (EEM) 9,10 , and the explicit electron model. 11,12 It should be noted that several models similar to the EEM exist, such as fluctuating charges (FQ or FlucQ), [15][16][17][18][19] charge equilibration (QE, QEq, 20 or CHEQ 21 ), and chemical potential equalization (CPE) 22 . Note that EEM is often not considered to be a proper PFF: because it only uses variable atomic charges, it cannot describe out-of-plane polarization of planar molecules. 22 One can represent these effects in EEM by including additional off-nuclear sites. [23][24][25] This leads to slightly more expensive models and additional model choices, such as the location of the extra sites and the associated parameters.
Here, we will treat EEM on the same footing as other PFFs because its mathematical structure is completely analogous.
The EEM was first developed by Mortier et al. from basic DFT equations according to Sanderson's principle of electronegativity equalization, 26 i.e., electrons flow until all electronegativities are equalized during the formation of a molecule. The ability of EEM (or one of its close rela-tives) to predict ground-state charge distributions has been validated extensively, e.g. for inorganic liquids, 16 inorganic solids, [27][28][29] organic molecules, 18,[30][31][32] and biomolecular systems. 19,[33][34][35] One appealing aspect of the EEM is that atomic monopolar fluctuations are included, and that one may systematically extend it with higher atomic multipoles, as in CPE. 22 However, the EEM also has two fundamental issues in simulations of extended systems or charge transfer during chemical reactions as follows: [36][37][38][39][40] (1) the EEM always predicts a cubic scaling of the dipole polarizability with system size, which is not correct for dielectric systems where a linear scaling exists in the macroscopic limit; (2) EEM obtains fractional molecular charges even when two molecules are well separated and gives a slowly decaying intermolecular charge transfer. The main reason for both problems is that long-range charge flow is always allowed in EEM, which is only realistic in metallic systems. Shortly after these shortcomings were first discussed in the literature, several ad hoc improvements were proposed, for instance, by applying artificial constraints to molecular charges 18,38 or dipoles, 41 or introducing a distance-dependent function to penalize long-range charge transfer. 42 Later, more consistent models were proposed, such as the split-charge equilibration (SQE): by adding a bond-hardness term to EEM, SQE becomes capable of modeling dielectric materials. 37,39 Recently, Verstraelen et al. developed a new polarizable force field (PFF), namely Atom-Condensed Kohn-Sham DFT approximated to second-order (ACKS2), 40 which has been validated and applied to molecular systems, and which addresses both fundamental issues of the EEM.
ACKS2 improves upon the EEM by adding an electronic kinetic energy term via the Legendre transform of the Kohn-Sham kinetic energy, enabling one to expand the kinetic energy to secondorder in the atomic populations and atomic Kohn-Sham potentials. It was shown that this new kinetic energy term is equivalent to the bond-hardness energy in the SQE. Later, a generalization of the ACKS2 model for arbitrary variational wavefunctions was proposed. 43 The ACKS2 formalism was numerically validated by a direct computation of the ACKS2 parameters from a reference Kohn-Sham DFT calculation. With these parameters, the static linear response properties of molecules obtained with an underlying theory, e.g., Kohn-Sham DFT, can be reproduced precisely in the limit of a complete basis sets for density and potential fluctuations.
To use ACKS2 for approximate polarization energy calculations, a tailored empirical parameterization was incorporated into ReaxFF, as an alternative to the EEM, which turned out to be essential in the development of eReaxFF. 44 Later, Gütlein et al. developed a new PFF using Gaussian basis sets, based on the ACKS2 model, 45 where generic parameters were proposed for elements C and H, such as widths of Gaussian basis functions. This parametrization can predict DFT data, such as response properties and interaction energies, of a series of hydrocarbons. Later, they proposed two ACKS2 variants, f-ACKS2 and scf-ACKS2, to handle condensed phases by dividing them into molecular fragments, which are coupled by Coulomb interactions. 46 This work extends ACKS2 to the time-and frequency-dependent domain to compute the dynamical response properties of finite molecules, e.g., frequency-dependent polarizabilities, C 6 dispersion coefficients, and molecular absorption spectra. The resulting ACKS2ω model, through its use of atomic multipole expansions, allows one to investigate the contribution of charge-flow effects on dynamic response properties, which is a topic of ongoing research. One of the future applications envisioned for ACKS2ω is an efficient approximation of the non-local dispersion energy, i.e. for which interatomic charge fluctuations become important, 47-54 such as carbon nanomaterials, 48,49 traditional semiconductors, 50,51 and low-dimensional materials. 48,[52][53][54] To our best knowledge, the time-dependent or frequency-dependent extension of ACKS2 has not been investigated yet.
The derivation of ACKS2ω relies on time-dependent DFT (TDDFT). TDDFT has received much attention in recent decades because of its good trade-off between numerical efficiency and accuracy compared to wavefunction-based methods. The standard formalism of TDDFT builds on the Runge-Gross theorem, 55 stating that there exists a one-to-one correspondence between densities and potentials for any fixed initial many-body state. TDDFT uses many DFT concepts, e.g. the whole interacting many-body system is replaced by a non-interacting Kohn-Sham system that shares the same density. However, it also has unique requirements, such as the initial state dependence and time-dependent exchange-correlation functional. Readers are referred to recent books 56,57 and review articles [58][59][60][61] for further details about TDDFT. Within the TDDFT regime, the linear response of the electron density can be calculated precisely, assuming one disposes of the exact exchange-correlation functional. In practice, TDDFT is a useful and relatively simple scheme to compute molecular linear-response properties. [62][63][64][65][66] To derive the ACKS2ω model, we use a quasi-energy formalism, also known as the Floquet theory, which is often used in time-dependent response model development [66][67][68][69][70][71][72] and can be seen as a special case of TDDFT. 55,73 The benefits of the quasi-energy language have been discussed elsewhere. 72 One of the advantages of quasi-energy DFT is that the initial state dependence of the Runge-Gross theorem is not necessary anymore, and instead a periodic boundary condition is applied to an external perturbation. Moreover, the quasi-energy ansatz enables a unified frame-work for dynamic linear-response property calculations based on wavefunction and DFT methods, as it is generally valid for both variational and non-variational methods. 72 This allows dynamical response properties to be derived by the virtues of the energy derivative approach using the quasi-energy formalism. However, the validity of Floquet theory has been a point of discussion elsewhere [74][75][76][77] and principle requires the following approximations: (i) a finite basis is employed; 76,77 (ii) the driving external perturbation is weak and off-resonant. [76][77][78] Point (i) ensures an adiabatic limit exists for the "Floquet ground state", on which the energy minimum principle is valid, and the latter determines the one-to-one density-potential mapping. Point (ii) suggests that Floquet theory is only a proper approximate method to treat liner-response problems. Nonetheless, for molecular force-field development, these two conditions are generally satisfied; thus, one can still use Floquet theory to develop frequency-dependent force-field methods.
The goal of this paper is to introduce ACKS2ω as a model with a structure reminiscent of conventional PFFs, yet with all parameters defined in terms of an underlying electronic structure theory. This direct connection to electronic structure theory assures a solid foundation for more pragmatic parameterizations of frequency-dependent PFFs. [79][80][81][82][83][84] The main distinction with the work of Misquitta and Stone, e.g. ISA-Pol, 85 is that ACKS2ω explicitly describes the transformation from non-interacting to interacting response functions through electron-electron interactions, whereas ISA-Pol directly models the interacting response relying on time-dependent Hartree-Fock or TDDFT results. Potential applications of frequency-dependent PFFs include modeling of optical spectra of complex systems at a force-field cost, [79][80][81] approximations of dispersion in force fields, 85 computationally efficient non-local dispersion corrections including many-body effects for DFT calculations 82,83 or frequency-dependent polarizable embedding. 84  The remainder of this paper is organized as follows: the basic theory of the ACKS2ω model is introduced in Section II, followed by computational details of the numerical validation in Section III. Results and discussion are presented in Section IV. A summary is given in Section V. Atomic units are used, unless noted otherwise.

II. THEORY
This section derives the ACKS2ω formalism from TDDFT, using the quasi-energy formalism. Hence, a brief review of the quasi-energy method is first presented in Section II A. The quasi-energy expressions have been employed in response theory because of their straightforward definition. 70,89 The theory of ACKS2ω within quasi-energy ansatz is presented in Section II B, where a set of equations is derived in analogy with ACKS2. 43 A dynamic linear-response theory is then provided in the frequency domain in Section II C. As a result, a set of linear-response equations are generated when two finite basis sets are applied (Section II D). Finally, in Section II E, we apply ACKS2ω in Kohn-Sham DFT with a semi-local exchange-correlation (xc) functional.

A. Quasi-energy formalism
The quasi-energy formalism or Floquet theory is employed in the original TDDFT derivation. 55 The definitions of time-dependent exchange-correlation potentials and kernels are much more accessible using the quasi-energy notation. In this section, we use the notation of Salek et al. 72 and Christiansen et al. 70 From the time-dependent Schrödinger equation, we have whereĤ 0 andV 1 (t) are the time-independent zeroth-order Hamiltonian and time-dependent perturbation operator, respectively. The exact solution ofĤ 0 is given by a ground-state calculation and perturbation theory can be developed through a generalization of the Rayleigh-Schrödinger theory. 72 We consider the time-dependent perturbation theory 70,78 where F(t) is a purely time-dependent function. The vector |0 represents the complete wavefunction whose phase-isolated part is denoted by |0 , where time dependence of the vectors is understood. The time-dependent Schrödinger equation is then rewritten as in which L(t) =Ḟ(t) is the time-dependent quasi-energy. Langhoff et al. showed that |0 can be normalized at all times. 78 Consequently, the time-dependent quasi-energy can be expressed as In order to use Hellmann-Feynman theorem, in analogy to the static case, the external timedependent perturbation is assumed to be periodic, i.e.,V 1 (t) =V 1 (t + T ). 72 Here, we use the symbol T , referred to as the time period, instead of T used in Ref. 72 to avoid ambiguity with kinetic energy functional T , which is well-known in DFT or TDDFT. Using curly brackets for time-averaging, the quasi-energy is defined as Techniques for time-independent problems, such as the variational principle 70 as well as the Hellman-Feynman theorem, 72 can be applied to the quasi-energy, to obtain a solution or a response to a perturbation for the time-dependent case. The quasi-energy formalism is also generally applicable to wavefunction methods and DFT 71,72 , at least for weak and off-resonant perturbation and when using finite basis sets approximations. 76,77 In the context of TDDFT, the quasi-energy can be written a functional of the time-dependent density, ρ(r,t): 64,72 where A more extensive discussion of the quasi-energy formalism in DFT can be found in Ref. 72, 76, and 77.

B. Decomposition of the quasi-energy into explicit and implicit terms
The functional in DFT can be written as the sum of a universal functional and the interaction of the electrons with external potential. The universal functional is a general term for all systems, which can be divided into two parts suggested by the ACKS2 model: the explicit functional (E exp ), which is a known functional of the density, and the implicit functional (E imp ). Similarly, in quasienergy TDDFT the corresponding explicit functional Q exp and implicit functional Q imp are given as The implicit functional uses an auxiliary N-fermion wavefunction through a constrained-search With the method of Lagrange multipliers, it can also be written out explicitly where u(r,t) is a function that specifies a Lagrange multiplier at every point in terms of both For a given external potential, v(r,t), the N-electron ground state is solved by minimizing the following Lagrangian with respect to ρ(r,t) and maximizing it with respect to u(r,t) and µ(t): where the quasi-energy is now a functional of the density and auxiliary potential Finally, the stationary point is defined by the following sets of Euler-Lagrange equations, 64

C. Linear response
In linear-response theory, it is more convenient to first apply a Fourier-transform to all timedependent functions, 56 which yields formally similar Euler-Lagrange equations: In analogy with the ACKS2 method, a static DFT ground state is taken as the reference to which fluctuations are considered, including the external potential v 0 (r), density ρ 0 (r), auxiliary potential u 0 (r), and equalized chemical potential µ 0 . The main difference with the original ACKS2 method, is that we now consider frequency-dependent fluctuations: where ∆v(r, ω) is a frequency-dependent perturbation and ∆ρ(r, ω), ∆u(r, ω), and ∆µ(ω) are corresponding frequency-dependent responses to the perturbation. After substitution in the Euler-Lagrange equations, one obtains In the limit of a small perturbation, the first two equations can be linearized The second order functional derivatives in Eqs. (26) and (27) are the hardness kernel of the explicit functional and response kernel of implicit functional, referred to as η exp (r, r , ω) and

D. Expansion in a finite basis
In a practical simulation, one must expand the density and the auxiliary potential fluctuations in a finite basis, where f m (g n ) denote density (potential) basis functions, while C m (ω) (U n (ω)) denote the expansion coefficients of the induced density (potential) changes ∆ρ (∆u). It is worth mentioning that both basis sets are time-independent, and only the corresponding coefficients are dynamic, suggesting that the basis sets used in the ACKS2 model can also be employed herein.
Substitution of the basis-set expansion in Eq. (26), followed by a multiplication with f k (r) and integration over r, leads to: O kn = f k (r)g n (r)dr (32) Similar manipulations of Eqs. (27) and (25), lead to algebraic equations as follows: with The linear system can be written in block matrix notation as follows: where a bold letter denotes a block matrix and the subscript gives its dimension. For example, gives the expansion coefficients for the induced density, i.e., C M,1 (ω), which is also known as the response vector. The shape subscript will be omitted below for the sake of visual clarity.
According to linear response theory, the induced density can be expressed through an interacting response function, χ: where V (r, ω) is the frequency-dependent external perturbation. ACKS2ω can be used to approximate the interacting response, for which two strategies can be followed. The first and most straightforward option is to invert the block matrix in Eq. (38) directly. The top-left block of the inverse transforms the input V (ω) into the output C(ω), and is therefore the ACKS2ω approximation of the interacting response.
One may also derive a useful closed expression for the interacting response matrix. This derivation starts by solving the second row of Eq. (38) inverse of the implicit response matrix is needed, because it contains at least one zero eigenvalue, due to the charge distribution not being sensitive to a constant shift in Kohn-Sham potential. Next, one uses this solution to eliminate U (ω) from the first row of Eq. (38), yielding where we introduced a shorthand A = −η + O(χ imp ) * O T to facilitate the derivation. This equation can be solved, assuming A is non-singular: When A is singular, the density response is not well-defined. While we cannot exclude this possibility of A being singular, we never encountered this issue in our numerical results. The Lagrange multiplier is found by substituting this result in the third row of Eq. (38): The final expression for the interacting response is: The second term in χ (due to normalization) is responsible for a proper zero eigenvalue in the interacting response, which also guarantees the constraint of particle conservation, i.e., D T χ = 0.
is defined with respect to the density basis f m , whereas U (ω) is expanded in terms of g n . Hence, χ is the interacting response expanded in the density basis. If desired, one may also express it in the potential basis by assuming that perturbations in the external potential have the following form: The new expansion coefficients, V n (ω), can always be transformed to those from Eq. (33): With that, one can left-multiply Eq. (43) with O T and substitute V (ω) with OV (ω): where χ is finally identified as the interacting response matrix in the potential basis.

E. ACKS2ω matrix elements
So far, the precise forms of the implicit and explicit functionals were not specified, and one may use, in principle, any definition whose sum equals (or approximates) the total quasi-energy.
In this section, we show that for Kohn-Sham DFT, a convenient choice can be made, leading to straightforward equations for all coefficients in the ACKS2ω equations.
The Kohn-Sham DFT expression of the quasi-energy is where, T s [ρ] and S s [ρ] are defined using the determinant of the non-interacting system (|0 s ): The quasi-energy xc functional is thereby defined as Within the quasi-energy formalism, we employ the adiabatic approximation by replacing Q xc [ρ] with a time-independent counterpart {E xc [ρ]} T .
With a semi-local xc functional, the explicit part of the energy functional in Kohn-Sham DFT takes the following form: where, φ i is the i-th occupied (Occ.) molecular spatial orbital.
The expressions for the explicit hardness matrix and implicit response matrix with the Lehmann representation 56,57 take the following forms in frequency space, respectively: a∈Vir.
where, n i and n a are 1 and 0 for occupied and virtual (Vir.) molecular spatial orbitals, respectively.
After simple manipulations, the matrix elements of χ imp have a simplified form where Ω ia = ε a − ε i and the Dirac notation is used. For closed-shell molecules, a factor four is applied instead of two in Eq. (55) and the sum is restricted to spin-up orbitals only.

F. Dynamic polarizability and related properties
The dynamic linear-response properties, such as dipole polarizabilities, can be derived from the frequency-dependent response matrix. Moreover, with a proper choice of density and potential basis functions, distributed polarizabilities 90 are also available.
The potential basis set is constructed as where R m (r − R a ) is a real solid harmonic with angular quantum number and magnetic quantum number m, using the position of nucleus a as origin. The functions w a (r) (0 ≤ w a (r) ≤ 1) are socalled atom-in-molecule (AIM) weights functions, determining which proportion of the molecular ground-state density is assigned to atom a. 91 The rationale is that the integral of the product of such a potential basis function and a (response) density can be interpreted as distributed multipole moment.
The density basis set can be organized similarly, using atom-centered functions with unit monopole, dipole, etc. moments. One may use Gaussian basis sets, 45 or construct a basis that is bi-orthogonal to the potential basis set, 43 which is also the approach followed in this work.
Both basis sets are truncated by discarding all solid harmonics for which > max , where max is a user-specified threshold. For instance, max = 0 and max = 1 corresponds to fluctuating charges only and charges+dipoles, respectively. The total basis size is therefore 2( max + 1) 2 times the number of atoms.
For the remainder of this section, it is convenient to change the index notation and introduce a compound index t (or u) for the components of multipole moments, instead of separate indices and m. With these, we can define a distributed polarizability as: where χ k(a,t) n(b,u) are elements of the response matrix defined in Eq. (43). When considering the case max = 1, distributed polarizability comprises several physically distinct blocks: the distributed charge-flow, α ab qq , charge-dipole, α ab qβ , dipole-charge, α ab αq , and dipole-dipole, α ab αβ , components. The subscript q is used to denote the monopole block and α(β ) refers to dipole blocks, whereby both α and β could be x, y, or z. The total molecular dipole polarizability can be recovered as follows: 91 where r a α and r b β are corresponding atomic coordinates. The isotropic dipole polarizability of molecules is formally given as: Eqs. (58)-(59) remain valid for frequency-dependent polarizabilities. When a real frequency is used, α becomes complex (when iη = 0). The dipole strength function S(ω) can be computed using the imaginary part of α with 92 The excitation energies are thus determined from the positions of the peaks, and the area under the peaks represents the oscillator strength of the corresponding transition. In addition, the isotropic C 6 dispersion coefficient can be computed with imaginary frequencies (ω = iu) using Casimir-Polder equation: 93 where, superscripts A and B denote different molecules.

G. Two-site ACKS2ω model
Before studying the absorption spectra of molecular systems, it is instructive to start with a simple model system to illustrate the workings of the ACKS2ω approach in analogy to the TDDFT method. We consider a two-site system of an electric dipole along x-axis with coordinates of -1 and 1, and only s-type density and potential functions are considered. The two-site ACKS2ω problem can be solved analytically after making the following assumptions: where non-negative numbers a and b denotes atomic hardness and the classical Coulomb interaction between the density functions, respectively. The xc contribution to the hardness is neglected entirely, i.e., the RPA is used.

The non-interacting response matrix is
where Ω 12 is the excitation energy between the ground state and the first excited state, determining the pole position of the non-interacting response function. In addition, an arbitrary amplifier d is added to create a more general non-interacting response matrix.
One may solve the interacting response of this model analytically, by working out Eq. (43). Doing so leads to an interacting response matrix of the form Therefore, the pole of the interacting response function is where the condition a > b is generally satisfied with the RPA approximation, because the Hartree kernel is positive definite. This inequality implies that the pole of the interacting response is shifted to higher frequencies compared to the non-interacting response, when using RPA.

III. COMPUTATIONAL DETAILS
This section presents the details of the numerical validation of ACKS2ω using the TS42 database and its TS27 subset. As discussed in Section II D, both potential and density basis functions are frequency-independent. Here, the basis sets are constructed in the same way as in one of our previous works on ACKS2 43 . The MBIS method 94 was used to define the AIM weight functions, w a (r). Using the same procedure as in Ref. 43, the density basis set consists of non-interacting responses to each potential basis function, when used as a perturbation. This set is augmented with one Fukui function, to have at least one basis function that is not normpreserving. Later, a linear transformation is applied to the density basis functions, such that they become bi-orthogonal to the potential basis, enabling one to label them as s-type, p-type, etc.
This definition of the density basis set is relatively expensive, and is geared towards an accurate description of linear response functions. The goal of this work is to test the validity of the ACKS2ω theory, justifying the use of carefully constructed basis functions. The bi-orthogonality of the density and potential basis functions also simplifies the elements of ACKS2ω working matrix in Eq. (38). More specifically, the overlap matrix O becomes an identity matrix, and all non-zero elements of the vector D are equal to 1, corresponding to integrals from s-type density functions. Moreover, when the adiabatic approximation is used, the hardness matrix becomes frequency-independent and takes the same form as in the ACKS2 paper. 43 More specifically, integrals including the hardness kernel of the explicit functional were evaluated numerically using a pruned Becke-Levedev integration grid. 95 Integrals over the Hartree kernel were implemented with a Becke-Poisson solver. 96 Although the analytic expressions for LDA and PBE kernels are already provided in LibXC, integrals involving the exchange-correlation kernel were evaluated with the finite difference method used to be consistent with the previous work. 43 The only frequencydependent parameters of ACKS2ω are found in the non-interacting response tensor, which can be evaluated using Eq. (54) on the same numerical grids.
The workflow of the ACKS2ω parameterization in this work is analogous to our earlier assessment of the ACKS2 theory: 43 (1) The molecular structures listed in the T42 data set are optimized using DFT at B3LYP/aug- To demonstrate the accuracy of ACKS2ω on dynamic linear-response properties, we first studied absorption spectra of 42 molecules from the TS42 data set. We only estimate the absorption spectra at the PBE/aug-cc-pVDZ level to save computational resources. 300 real frequencies rang- In the following context, we use the notation ACKS2ω@X to specify the ACKS2ω model In comparison, the bottleneck of TDDFT is the transformation of electron repulsion integrals (ERI) from atomic orbitals (AO) to molecular orbitals (MO), which cannot be applied in largescale systems due to its significantly expensive O(N occ N 4 AO ) complexity from the conventional transformation or O(N 2 occ N 2 vir N aux ) from density fitting procedures where N AO is the size of the basis set and N occ , N vir and N aux ≈ 3N AO are the number of occupied, virtual orbitals and the size of auxiliary basis sets, respectively. [107][108][109][110] Finally, for the matrix inversion to calculate the interacting response matrix, a similar linearalgebra technology can be employed for both ACKS2ω and TDDFT with a semi-local exchangecorrelation functional. However, the dimension of ACKS2ω working matrix with l max = 1 is 4N atom , in practice much less than the TDDFT counterpart, i.e., N nocc × N nvir .
Because our current implementation is merely a prototype, we only illustrate timings of the

B. Comparison of ACKS2ω and TDDFT C 6 coefficients
This section presents a comparison of C 6 coefficients obtained with ACKS2ω and TDDFT, to validate the ACKS2ω formalism. Figure 3 shows a parity plot of C 6 dispersion coefficients for 903 molecular pairs, based on the TS42 data set, computed with s-type and sp-type ACKS2ω models.
The ACKS2ω parameters are derived from LDA/aug-cc-pVDZ calculations and LrTDLDA/augcc-pVDZ C 6 coefficients are used as a reference. All s-type ACKS2ω parametrizations significantly underestimate C 6 coefficients, averaging around 93% with respect to the reference values.
The reason for this large error is that the dipole polarizability is significantly underestimated due to less complete basis sets in s-type ACKS2, 43 . This effect is amplified in the isotropic C 6 coefficient, because it scales quadratically with the dipole polarizability. In contrast, very small errors can already be obtained with the sp-type ACKS2ω model, i.e. just considering fluctuating atomic charges and dipoles.   ACKS2ω parameters were derived using different xc functionals to study the impact of the functional. Some interesting observations can be made: (1) The positive MPE implies that all computational models in this work slightly overestimate C 6 values to some extent.
This suggests LrTDLDA systematically overestimates C 6 values when extensive basis sets are used. The MAPE (and MPE) increases with the size of the basis, and it is found to be converged at around 7%.
(3) The MAPE (MPE) of ACKS2ω has a similar trend as LrTDLDA, but two significant differences should be pointed out. First, the ACKS2ω model has a larger MAPE than MPE, indicating that the overestimation is not systematic in the ACKS2 model, i.e., for some molecules, ACKS2ω underestimates C 6 coefficients. Second, ACKS2ω@LDA has a lower MAPE (MPE) than the corresponding LrTDLDA for all cases, implying that the improvement compared to LrTDDFT is somewhat systematic. One possible expansion is that TDDFT always gives an overestimated C 6 with LDA functional, and the incomplete basis in ACKS2ω compensates for the overestimation of TDDFT to some extent, leading to a lower MAPE. This type of error compensation is obviously fortuitous, and should therefore not be relied upon blindly.
(4) The comparison between ACKS2ω@LDA and ACKS2ω@PBE shows a systematic improvement when using PBE instead of LDA. One possible reason could be that the contribution of correlation functional can be described more precisely with PBE than LDA. This inspired us to also test the ACKS2ω@LDAx model, which entirely ignores the correlation contribution in the hardness kernel, compared to ACKS2ω@LDA. As expected, it shows a significant deviation from the ACKS2ω@LDA model. Interestingly, ACKS2ω@LDAx performs slightly better than ACKS2ω@PBE, except for the aug-cc-pVDZ basis, where the latter has a lower MAPE (3.84%). Besides the potentially beneficial error compensation, these results also show that ACKS2ω is quite sensitive to the choice of functional, because the TDDFT reference also exhibits this sensitivity.  overestimated and underestimated, respectively. As discussed above, LrTDLDA overestimates C 6 coefficients and ACKS2ω@LDA slightly compensates for the overestimation due to the less complete bio-orthogonal basis sets. However, the compensation is too small to offer a significant improvement compared to TDRSHLDA. Interestingly, the results are improved and even better than PBE by removing the electron correlation contribution in response calculations, which again demonstrates that the effects of xc functional on C 6 coefficients with ACKS2ω are of significant importance.

V. CONCLUSIONS AND OUTLOOK
A frequency-dependent ACKS2 model, referred to as ACKS2ω, is proposed using the quasienergy formalism. This allows us to approach ACKS2ω directly using the variational principle, in the same way as the time-independent ACKS2 derivation. For the xc contribution to the hardness, an adiabatic approximation is applied, and the hardness is thereby frequency-independent, just as in the ACKS2 model. For the frequency-dependent non-interacting response matrix, a similar Lehmann representation is employed as in the static case. Given the hardness and non-interacting A comparison of ACKS2ω, using the LDA and PBE kernels, to experimental data and results obtained with range-separated calculations, provides some additional insights. The deviation of the ACKS2ω C 6 coefficients from the experiment is quite sensitive to the choice of the functional, with PBE and LDAx (correlation omitted from the hardness kernel), being favorable choices.
The best ACKS2ω model, ACKS2ω@PBE with the aug-cc-pVDZ basis, gives C 6 coefficients of molecules with an MAPE of 3.84%, a slight improvement over the TS method (6.3%). This improvement is significant because it potentially provides a better way to describe non-local dispersion energy between molecular dimers, which is a topic of ongoing research. In addition, the ACKS2ω@LDAx model has only a slightly inferior accuracy of C 6 coefficients than TDRSHLDA response calculations, with MAPEs of 5.29% and 5.17%, respectively, for 27 homodimers from the TS27 database.
The strength of the ACKS2ω model is that, once parameterized, calculations involving response kernels have a low computational cost, comparable in complexity to conventional polarizable force fields. In this work, the parametrization is relatively expensive to obtain accurate predictions, thereby showing that ACKS2ω is capable of reproducing its LrTDDFT reference.
In future work, simpler force-field like parametrizations could be used instead, to overcome the computational bottleneck of the current parametrization. We anticipate that this approach will be beneficial, e.g. for calculations of long-range correlation energies or optical properties of extended systems, for which DFT would be infeasible.

SUPPLEMENTARY MATERIAL
Additional display items showing spectra for the remaining molecules of the TS42 set, not included in Fig. 2

DATA AVAILABILITY STATEMENT
The data that supports the findings of this study are available within the article and its supplementary material.