Efficient representation of the linear density‐density response function

We present a thorough derivation of the mathematical foundations of the representation of the molecular linear electronic density‐density response function in terms of a computationally highly efficient moment expansion. Our new representation avoids the necessities of computing and storing numerous eigenfunctions of the response kernel by means of a considerable dimensionality reduction about from 103 to 101. As the scheme is applicable to any compact, self‐adjoint, and positive definite linear operator, we present a general formulation, which can be transferred to other applications with little effort. We also present an explicit application, which illustrates the actual procedure for applying the moment expansion of the linear density‐density response function to a water molecule that is subject to a varying external perturbation potential. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.


Introduction
Almost any situation in physics and chemistry can be described by a state function and its temporal evolution. The latter is generally governed by a differential equation, for example, the Newton equation of motion, the time-dependent Schrödinger or Dirac equation, the Navier-Stokes equation, the Maxwell equations, and others. In many of these situations, a viable strategy for simplifying the numerical solution of the underlying differential equations is the separation of both equation and solution into two parts of "large" and "small" amplitudes, respectively, which are then technically easier to deal with. The formal framework of this idea is commonly called perturbation theory, and the formal solution to this type of problem is often called Green's function.
The standard perturbation theory ansatz is based on the assumption that there is one given "large" part and one given "small" contribution of the differential operator, so that two new specific differential equations result from this large/small separation. However, there are situations in which a whole variety of "small" contributions exists (but only one "large" differential equation). Examples for these situations from molecular quantum chemistry are the calculation of vibrational modes, where there is one perturbation Hamiltonian for the displacement of each atom in each Cartesian direction or molecular dynamics simulations within embedding schemes, that is, where a given molecule is subject to varying chemical environments. In such situations, the same perturbation theory equation has to be solved repeatedly for numerous different perturbations.
In the specific case of quantum chemistry, the formal solution to this problem is called the density-density response function, which means a linear operator that acts on the perturbation and yields the electronic density response to that perturbation. This function is very general, as it yields the response density for any perturbation potential, but in turn, it is not readily available explicitly in numerical form. Nevertheless, there is a path for the explicit calculation of this complex quantity by means of a particular mathematical representation (using the eigenfunctions and eigenvalues of the density-density-response function seen as a matrix operator). [1,2] The straightforward calculation has been shown to be accurate yet computationally cumbersome. Here, we provide a mathematical study of a recently developed, considerably more efficient variant for the computation and the storage of this quantity (called the moment representation of this response function [3,4] ).
Compact linear operators often appear as integral transforms and can be viewed as continuous generalizations of matrices, where the corresponding integral kernel fulfills some weak conditions, that is, the kernel must not be singular and must decay to zero fast enough at infinity. The compactness of these operators ensures that the matrix elements are well defined. For X R 3 and χ(r, r 0 ) 2 L 2 (X × X) eq. (1) defines a linear Hilbert-Schmidt operatorT : L 2 X ð Þ !L 2 X ð Þ, which is in particular compact. [5]T : f r ð Þ↦ ð χ r,r 0 ð Þf r 0 ð Þd 3 r 0 ð1Þ If the kernel is furthermore symmetric (χ(r, r 0 ) = χ(r 0 , r) 8 r, r 0 2 X) and positive ( Ð Ð χ(r, r 0 )f(r 0 )f(r)d 3 r d 3 r 0 ≥ 0 8 f(r) 2 L 2 [X]) the operator is self-adjoint and positive.
The properties of this class of operators are well-understood from a mathematical point of view. In physics or quantum chemistry, these integral transforms are often derived from the inverse HamiltonianĤ − 1 and have many applications, for example, the Lippmann-Schwinger equation for the calculation of scattering in quantum mechanics. The Lippmann-Schwinger equation relates the scattered wave function with the interaction that produces the scattering (the scattering potential). Another example is the frequency dependent density-density response function, which includes an externally defined energy shift of the molecular orbital spectrum (originating from a photon of given frequency), which in turn leads to a breakdown of the linear response regime when the photon energy matches an electronic excitation energy. While this latter situation is within our view, we are still working on how to adapt the present representation of the linear density-density response function to this energy-dependent version.
In this article, the operator defined by the static linear density-density response functionT : V pert r ð Þ↦n resp r ð Þ will be considered as an example for the illustration of our theorems for linear, self-adjoint and positive definite operators.
For any given specific perturbation potential V xy pert r ð Þ the corresponding electronic response density n resp xy r ð Þ can be calculated via standard quantum chemical approaches, for example, by means of a perturbation-theory based calculation. [6][7][8][9][10][11][12] The challenge, however, is to develop a formalism for any (arbitrary) perturbation potential to a level where no explicit quantum chemical calculation is needed any more once the perturbation potential is specified. The core application for this approach will be molecular dynamics simulations, where the potential experienced by a given molecule varies slightly from any molecular dynamics step to the next. Once the kernel χ(r, r 0 ) is calculated, the calculation of molecular response density breaks down to integration according to eq. (1). The computational costs of this procedure are negligible compared to a single quantum chemistry calculation.
Unfortunately, there is no closed expression for the calculation of the kernel χ(r, r 0 ) and hence for the calculation of density responses. The (molecular density-density) response function is a six dimensional nonlocal kernel χ(r, r 0 ) and its storage on a grid would also require excessive storage dimensions, which decrease the efficiency of any real application.
A possible attempt to overcome these problems is to apply the spectral theorem of compact and self-adjoint operators.
Thus, we can express the action of the operatorT to a vector f(r) as:T using eigenvalues λ i and eigenfunctions χ i (r), which can be stored without difficulties. [1,2,13,14] The eigenfunctions and eigenvalues satisfyTχ i r ð Þ = λ i χ i r ð Þ. hχ i , fi denotes, the inner product of functions χ i (r) and f(r) in L 2 (X) and is defined According to the spectral theorem, the sequence of the eigenvalues converges to zero (lim i!∞ λ i = 0), which allows approximation of eq. (2) via a finite sum, that is, the truncation of the eigensystem representation at a finite index to the desired numerical precision. In particular, the sequence of eigenvalues is square summable, as the operatorT belongs to the second shadow class as a Hilbert-Schmidt operator. [5] Thus the eigenvalues decay fast enough to satisfy the condi- For our specific example ofT : V pert r ð Þ↦n resp r ð Þ with χ(r, r 0 ) equal to the molecular density-density response function, there is a suitable way to compute and store the most significant subset of eigenvalues and eigenfunctions. We were able to calculate the eigensystem representation and could already show that the calculation of converged response densities would require several thousand eigenvalues and eigenstates. [1,2] We want to emphasize that we furthermore have to calculate several thousand scalar products for the evaluation of each single perturbing potential. This is a computationally demanding protocol.
In this article, we present a transformation of this subset of eigenstates, which condenses the physically relevant information into a new set of states of considerably reduced dimensionality (moment expansion). Our aim in this manuscript is to provide a thorough mathematical basis for this transformation by means of two theorems.
For the specific case of the static linear density-density response function and a specific basis of the perturbing potential, we recently published a recipe for the calculation of the eigensystem representation [1,2] as well as the moment expansion. [3,4] Here, we generalize these results to an entire class of linear operators and arbitrary basis sets.

Conditions/Restrictions for the Application of the Theorems
In this article, we focus on linear operators defined by integral transforms according to eq. (1). Thus, both domain and image consist of square integrable functions. A closer inspection of eq. (2) reveals that only domain vectors f(r)with an nonvanishing overlap with the eigenstates (9 i 2 N : hχ i , fi 6 ¼ 0) will contribute to the image vector. Unfortunately, we do not have any information about the shape of the eigenstates in general.
However, we assume the linear operatorT will have compact support, that is, For concrete evaluations of the operatorT, we have to choose a basis {P n (r)| n 2 N} for the representation of a domain vector f(r): with c n being the n-th basis expansion coefficient. A precise approximation of a function f(r) in a local environment (at a compact subspace) is often possible within a finite number of basis functions. For example, an analytic function can always be expanded within a Taylor series (using the monomials as basis of the vector space). The terms of the Taylor expansion are ordered with respect to descending importance and it is even possible to approximate the introduced error by a truncated expansion.
We can validate, if the conditions for the domain vectors are fulfilled for a specific choice ofT, for example, for the operator T : V pert r ð Þ↦n resp r ð Þ obtained by setting the kernel χ(r, r 0 ) to the molecular density-density response function. This operator links an arbitrary perturbing potential to the resulting molecular response density. Due to the locality of the quantum mechanical Hamilton operator, these response densities occur only in the vicinity of a molecule, which corresponds to the condition of compact support. Thus, we are only interested in an (approximate) representation of the perturbing potential in an area around the responding molecule, that is, in particular not around the physical origin of the perturbation.
A truncated multipole expansion (e.g., in terms of Racah normalized irregular solid harmonics [15,16] ) allows such a finite representation of the perturbing potential.
A demonstration of the truncated expression of the perturbing potential of a water molecule is presented in the last section of this article.

Non-Orthogonal Basis Functions
The Dirac notation is a standard tool in quantum chemistry in order to express manipulations of linear operators and wave functions in an abstract formalism. In particular, it allows the efficient transformations of the representation with respect to different bases. The theoretical justification of the Dirac notation is the Riesz representation theorem, which states that all separable infinite Hilbert spaces are isometric isomorph.
We defined the operatorT as linear map between two infinite separable Hilbert spaces, which justifies the use of the Dirac notation even through the physical meaning of a bra/ket vector differs from the conventional one (potentials and densities as opposed to wave functions). Here we will briefly review the Dirac notation, in particular highlighting deviations from standard notation for non-orthogonal basis sets.
The resolution of identity Id jpi can be expressed for an orthogonal basis set {| p n i | n 2 N} as: and for a non-orthogonal basis {| q n i | n 2 N} we have to write: with S jqi being the overlap matrix with the elements (S jqi ) nm = hq n , q m i. The matrix S q is known as Gramian matrix of the vectors {| q n i | n 2 N} and it is a non-singular matrix, if the {| q n i | n 2 N} are linearly independent. Thus, the existence of the inverse S −1 jqi is always ensured. A derivation of eq. (5) and an explanation of the treatment of non-orthogonal basis sets is given in the Supporting Information.
Within an orthogonal basis {| p n i | n 2 N}, the linear opera-torT can be expressed aŝ with the element of the transformation matrix: Herein, T The index k in {p k } will only be used to refer to the entity of basis vectors {| p k i| k 2 N}, but is not for a specific indexing.
Within a non-orthogonal basis {| q n i | n 2 N} of the image and an orthogonal basis {| p n i | n 2 N} of the domain, the linear operatorT can be expressed aŝ with the element of the transformation matrix: Derivation of the Moment Expansion be a compact, positivedefinite and self-adjoint linear operator with eigensystem representation from eq. (2). Further, let {P n | n 2 N} be an orthonormal basis of the domain ofT. Then: and 2) The operatorT can be expressed aŝ The restrictionT Proof: Proof of sub-item 1: We choose the non-normalized ffiffiffi ffi as basis of the image and the orthonormal {P i | i 2 N} as basis of the domain ofT. According to eq. (9) we can express the elements of the infinite transformation matrix as following: As the transformation matrix ofT represents a bounded linear operator on the separable Hilbert space ℓ 2 , there is a QR into the product of an orthogonal matrix Q and an upper triangular matrix R: [17,18] T ffiffiffi ffi Equation (16) reads element-wise: ð17Þ ffiffiffi ffi λ l p χ l ,P j we obtain: We define so that R ij = hξ i , P j i. The fact that R is an upper triangular matrix implies that: hξ i , P j i = 0 8 i > j. The last expression proves equation (10) from sub-item 1.
In the next step, we prove equation (11) from sub-item 1. Equation (19) shows that Q represents a change-of-basis matrix which we will denote by W ffiffiffi ffi We can write for the inverse basis transformation: On the one hand, applying the basis transformation to yields by definition (compare commutative diagram of standard basis transformations in the Supporting Information): and on the other hand, the result of the transformation can be obtained from eqs. (16) and (20): Combining eqs. (21) and (22) we get: T Proof of sub-item 2: We choose the non-orthonormal {ξ i | i 2 N} as basis of the image and the orthonormal {P i | i 2 N} as basis of the domain ofT. According to eq. (11) we can express the operator T as: For arbitrary elements jfi of the domain ofT, we can write: Finally, we have to show: the restrictionT P N ½ is obtained by truncating eq. (26) after N terms. For an arbitrary vector j f i 2 P N ½ with j f i = P N n = 1 j P n i P n ,f h i, we have: □ General Remarks: 1. Sub-item 2 of this theorem corresponds to the following statement: The complete image due to domain vectors from the restriction f 2 P N ½ can be calculated using eq. (13) within N terms. Remarks for the special case ofT : V pert ↦n resp : 1. For the specific choice of Racah normalized irrgelular solid harmonics as basis functions for the expression of the pertubing potential, Scherrer et al. [3] showed in 2016, the existence of the moment expansion. In this article, we put the moment expansion on a more solid mathematical fundament and generalize it to arbitrary self-adjoint, positive, and compact operators (and arbitrary bases of the domain of these operators). 2. The set {P n | n 2 N} denotes the basis functions for the domain ofT, which is different for every molecule. For practical application we aim to use standard basis sets of quantum chemistry such as solid harmonics or polynomials. The support of these standard basis sets is in general different from the compact support of the linear operatorT. Thus inner products hP i , V pert i are taken on the wrong vector space, which introduces some problems for actual calculations of the response density. The problem can be circumvented by expressing V pert within the set of {P 1 , P 2 , …,P n } and employing eq. (13), because the ξ i vanish outside the domain ofT. The overlap ξ i ,V pert 8 V pert 2 P N ½ is automatically taken at the correct subspace. A detailed discussion of this problem is given in the Supporting Information.

Direct Moment Expansion
Ahlert et al. [4] developed already an idea how to obtain the moment expanded states without calculation of thousands of eigenstates. In the next theorem we put this approach to the next level and show how to obtain the first N moment expand states within N ab inito calculations (by means of the a cholesky decompostion of a N × N matrix): Theorem 2: LetT be the linear operator, {| P n i | n 2 N} the orthonormal basis of the domain ofT and {| ξ n i | n 2 N} the moment expanded states from Theorem 1. We define and two matrices (Θ ξ ) and Θñ À Á with the elements: Then: 1Þ Proof of the first sub-item: We start with the definition from eq. (29): Equation (32) shows that Θ ξ is a change-of-basis matrix. The inverse of the change-of-basis matrix is the change-of-basis matrix of the inverse basis transformation: Θ ξ is an upper triangular matrix (compare eq. (10) and eq. (30)), which in turn ensures the full rank of the matrix and hence the existence of the inverse (Θ ξ ) −1 .
Proof of the second sub-item: The symmetry of Θñ is obtained from: Multiplying eq. (32) from the left with hP j j we obtain: This is equivalent to the following matrix equation: Equation (37) actually defines a Cholesky decomposition of the symmetric matrix (Θ ξ ) into a product of an upper triangular matrix (Θ ξ ) and its transposed (Θ ξ ) T . This upper triangular matrix is nothing but the moment matrix of the moment expanded states.
□ Remarks: 1. A Cholesky decompostion is a decomposition of a symmetric, positive-definite matrix into the product of a lower triangular matrix and its transpose. We assumedT as positive operator. Thus we can obtain the symmetry as well as the positive definiteness of Θñ from Θñ À Á ij =ñ i ,P j = P i ,TP j .
2. In general Θñ and Θ ξ are infinite matrices. In the next remark, we show that only the N × N restrictions of these matrices are necessary for the calculation of the first N moment expanded states {| ξ 1 i, | ξ 2 i, …, | ξ N i}. 3. According to the last theorem, we want to summarize the effort for the calculation of the first N moment expanded states {| ξ 1 i, | ξ 2 i, …, | ξ N i}. We have to calculate: • N evaluations ofT according to

Algorithm of the direct moment expansion
Based on the second theorem, we want to report an explicit protocol for the calculation of the moment expanded states: • Choose a basis {| P k i | k 2 N} for the domain of the linear operatorT .
• Compute the set of functions jñ i i≔T j P i i j i 2 1, …,N f g È É .
• Calculation of an N × N matrix with elements Θñ À Á ij =ñ i ,R j .
For the actual transformation two steps are necessary: First step: Cholesky decomposition (via library algorithms e.g., lapack) of Θñ = Θ ξ À Á T Θ ξ À Á . Second step: The inverse of the N × N moment expanded moment matrix (Θ ξ ) −1 transforms the set {ñ i |i ≤ N} into the moment expanded states {ξ i |i ≤ N} according to eq. (33).
A schematic representation of the direct moment expansion is given in Figure 1.

Application of Theorem 1 and 2 to the Linear Density-Density Response Function General description
The calculation of the electronic response due to a perturbation potential is common to numerous important applications from spectroscopy to intermolecular interactions ( Figure 2). One example is fragmentation methods, where a large system is partitioned into smaller fragments that are more efficient to compute. The mutual interaction between the fragments can be realized via a perturbation theory approach, where each fragment perturbs the other ones via its electrostatic potential (acting as a perturbation potential). [19][20][21][22][23][24][25][26] Following these ideas, many embedding or subsystem density functional theory (DFT) methods [27] like frozen density embedding [28][29][30][31][32][33][34] or embedded correlated wave functions methods [34][35][36][37] were developed in the literature.
Accurate electrostatic interaction energies should take polarization effects into account. For the specific choice of the molecular linear density-density response function as kernel χ(r, r 0 ) the integral transform in eq. (1) constitutes an operatorT : V pert r ð Þ↦n resp r ð Þ which maps perturbing potentials V pert (r) to response densities n resp (r). These response densities fully incorporates polarization effects.
In Theorem 1 and 2 we prove the validity of an efficient protocol for the calculation of the density response due to approximated/truncated perturbing potentials V pert r ð Þ 2 P N ½ ≔Span P 1 r ð Þ,P 2 r ð Þ,…,P N r ð Þ ð Þ using the operatorT. We present a generalized representation (i.e., independent of the specific shape of the perturbing potential) of the response density, which can be computed and stored with a reasonable computational effort by means of the moment expansion. The latter is designed to optimally reflect the shape of physically meaningful (yet arbitrary) perturbation potentials which will be generated by nearby molecules.

Efficiency gain due to moment expansion
The quantum chemical calculation of the (linear) response of a molecule to an external perturbation potential can be done (1) via a conventional self-consistent perturbation theory calculation [6][7][8][9][10][11][12] or (2) via application of the generalized densitydensity reponse function χ(r, r 0 ). For the determination of a single response density for one given perturbation potential (or at least a limited number of such perturbations), the direct calculation (1) is more efficient than our approach (2) due to the comparably high computational cost of obtaining χ(r, r 0 ) explicitly. The actual application of χ(r, r 0 ) to any specific perturbation potential in turn is quite small, because only a limited number of scalar products of the specific perturbation potential function with the basis functions for χ(r, r 0 ) is necessary (in particular, no self-consistency cycle is involved any more). Thus, the advantage of using an explicit representation of the density-density response lies in repeated applications of response calculations for numerous different perturbation potentials for the same molecule, that is, when only one initial determination of χ(r, r 0 ) is required. For the original eigensystem formulation of χ(r, r 0 ), the break-even point depends on the number of eigenstates that are needed for an accurate representation, typically around 1000-5000 (e.g., for a water molecule). The number of required self-consistent perturbation theory cycles for the initial determination of χ(r, r 0 ) roughly equals the number of desired eigenstates.
The efficiency gain of the moment expansion stems from the reduction of that number (the number of eigenstates required for an accurate representation of χ(r, r 0 ) by about two orders of magnitude. In our new representation, the focus is shifted from optimally representing the local properties of χ(r, r 0 ) toward optimally adapting the representation to the shape of the perturbation potential V pert (r). Hence, our moment expansion shifts the break-even point compared to conventional self-consistent perturbation theory calculations from 1000 to 5000 toward a number of about 10-100 conventional calculations. While this figure of course still means that for a small number of response calculations, the conventional path is more advantageous, the perspective changes when a very large number of response calculations for the same molecule is required. A case in point for this situation is molecular dynamics simulations of a molecular liquid. From the perspective of a given molecular in the liquid, each new molecular dynamics frame constitutes a new chemical environment, to which the molecular density will respond. This response is currently computed self-consistently (Born-Oppenheimer molecular dynamics) or iteratively (Car-Parrinello scheme). Our approach could be used to switch the determination of the density response to a new molecular dynamics step to a nonself-consistent way at considerably reduced computational cost without sacrificing quantum chemical accuracy.
What is the physical interpretation of this dimensionality reduction? The spatial shape of the eigenstates of the densitydensity response function resembles a kind of symmetrized set of molecular orbitals. The eigenstates often appear like atomcentered (or bond-centered) functions of localized nature and strong nodal structure, that is, strong local oscillations. Furthermore, subsequent eigenstates often cover the same spatial area but just show a higher nodal structure. From a more general perspective, an analogy to fitting a function via a polynomial interpolation comes into mind: Fitting N data points yields an interpolation polynomial of order N, which however is rarely the optimal choice for a relatively smooth function. This problem of the original eigensystem representation of the densitydensity response function is cured by the moment expansion formulation. The specific impact of the two theorems for the linear density-density response function is discussed in the following.

Relevance of Theorem 1
The eigensystem representation offers a straightforward but inefficient possibility for the calculation of the density response. The infinite expansion in eq. (2) can be truncated after M eigenstates to the M largest eigenvalues (M~10 3 − 10 4 ), [1,2] as the sequence of the eigenvalues decays to zero. In order to obtain the response density due to a perturbing potential V pert (r) 2 Span(P 1 (r), P 2 (r), …, P N (r)), the calculation of M dot products/projections between the eigenstates and the specific perturbing potential V pert (r) is necessary in eq. (38).
Utilizing the moment expanded states resulting from Theorem 1, a drastic reduction of the dimensionality of the expression of the response density is achieved. The full density response due to V pert (r) 2 Span(P 1 (r), P 2 (r), …, P N (r)) can be calculated within N moment expanded states with N < < M (subitem 2 of Theorem I): Relevance of Theorem 2 The second theorem gives rise to a much more efficient algorithm for the calculation of the moment-expanded states (referred to as direct moment expansion). For the brute force construction of the transformation Q derived in theorem 1 (which corresponds to the algorithm of the moment expansion published in 2016 [3] ) a large number of M eigenstates are calculated, followed by M × N Givens rotations to condense the physically information. Via the direct moment expansion, we are able to calculate the first N moment expanded states from the first N direct density expanded states (by a Cholesky decomposition of a N × N matrix). Each direct density expanded state can be computed via a single quantum chemical perturbation theory calculation. A concrete, step-by-step protocol for the calculation of the direct moment expansion is given in the previous section.

Numerical Example: Density Response of a Water Molecule
In this section, we will demonstrate that the linear operator defined by the linear density-density response function fulfills the conditions for the application of the moment expansion, that is, we will show that the relevant domain vectors (perturbing potentials arising from neighboring molecules) can be expressed within a finite number of basis functions.
We choose the perturbative effect of one water molecule on an adjacent one as an elementary example (Figure 2). For simplicity, we use the Hartree potential using the partial charges of one of the most common force field water models (TIP3P) as perturbation, while the perturbed water molecule is represented within standard density functional theory, complemented with the density-density response function computed using density functional perturbation theory calculations (as practial realizations of the operatorT). The Hartree potential of the TIP3P partial charges reads: [49] V frag TIP3P r ð Þ = X Nn For 10 different distances to the origin of the TIP3P water molecule, we choose six representative points (spanning an octahedron) and expanded the perturbing potential around these positions using up to 54 monomials. Subsequently, we calculated the overlap of the polynomial approximated potential and the Hartree potential on a 4 Å × 4 Å × 4 Å cube centered at the origin of the potential expansion-note, that the origin of the perturbing potential (center of mass of a TIP3P water molecule) and the origin of the potential expansion (center of the responding water molecule) do not coincide-The overlap O jf i jf i of a function jf i and jfi was taken according to: Finally, we averaged the six overlaps of the Hartree potential and the expanded potential for each distance to the origin and presented the results with respect to different numbers of basis functions in Figure 3.
Furthermore, we demonstrated that we can use the moment expansion in order to express the full density response of water . Averaged and normalized overlap of the Hartree potential and the potential obtained from the expansion in monomials. The different colors refer to different number of basis functions. Please note, that the origin of the perturbing potential (center of mass of a TIP3P water molecule) and the origin of the potential expansion do not coincide! The overlap was taken at a 4 Å × 4 Å × 4 Å cube centered at the origin of the expansion of the perturbing potential. For each distance to the spatial origin of the potential expansion, we averaged over six overlaps (corresponding to six different spatial origins of the expansion of the perturbing potential, which span an octahedron). [Color figure can be viewed at wileyonlinelibrary.com] molecule within a few moment expanded states. For each distance to the perturbing TIP3P water molecule, we reused the octahedral alignment of the origins of the expansion of the perturbing potentials from the previous paragraph as center of masses of a responding water molecule.
Due to the perturbing potential, we calculated the response density of these water molecules via the moment-expanded states and via Density Functional Perturbation Theory (DFPT) (as reference calculations).
We calculated the overlap of the resulting densities and averaged all overlaps corresponding to the same distance to the perturbing TIP3P water molecule. The obtained values are depicted in Figure 4. From Figure 3 and Figure 4, we can deduce that already within nine states we can obtain reasonable good approximations of response densities or perturbing potentials. For small intermolecular distances, we can systematically improve the approximation by addition of further moment expanded states/basis functions. Within 54 states we obtain more than 90% of the response density/ the perturbing potential, even in the regime of the first maximum of the waterwater radial distribution function (~2.5 Å). This is a remarkable reduction compared the employment of several of thousands eigenstates of the linear density-density response function.

Conclusion
In this article, we generalize the concept of the moment expansion to arbitrary compact, positive and self-adjoint linear opera-torsT. This class of linear operators includes all integral transforms with a positive, symmetric and continuous kernel χ(r, r 0 ). We show that only comparably weak conditions have to be fulfilled (such as functions of the domain has to be Taylor expandable in an local environment) in order to allow for an efficient explicit representation (i.e., calculation and memory storage) of the linear operatorT, which enables the evaluation ofT with a drastically reduced costs compared to its eigensystem representation.
We have demonstrated that the moment expansion can be considered as a QR decomposition of the transformation matrix. The straightforward calculation of the infinite orthogonal matrix Q would in particular include the calculation of thousands of eigenvalues and eigenstates. We presented an efficient version of the moment expansion (termed as direct moment expansion), which allows the calculation of the first N moment expanded states at the computational costs of about N evaluations ofT.
For the specific case of the linear density-density response functions as kernel χ(r, r 0 ), we presented a numerical example for the efficient calculations of converged molecular response densities via the moment expansion. We demonstrated that perturbing potentials from neighboring molecules could be truncated within a few (ffi 34) basis functions.
The actual scope of our work exceeds far beyond this specific example. All Green's-function based formalisms in physics and chemistry have the problem that the Green's function itself can rarely be expressed/computed explicitly because of its dimensionality. In most cases, it is merely considered as a handy nomenclature for formally expressing that a linear differential equation can be solved (inverted) for a specific right-hand-side function. Our work shows that this problem can possibly be overcome in many situations. While we do not believe that we have a comprehensive overview of all possible applications ourselves, an indicative list of quantum chemical problems would include (a) fragmentation/embedding approaches, (b) van-der-Waals interaction energies via a Lundqvist expression [50] and correlation energies within the random phase approximation, and (c) resolution of identity approaches. Possible fields of application beyond this directly related area are (d) polarizable force fields, (e) optimal control theory, and (f) signal processing (impulse response function). However, we believe that our Ansatz may also find applications in the wider area of "fundamental solutions" (in mathematical terms) where spectral theory/Fredholm integral equations are applicable.
Keywords: density-density response function Á molecular interaction Á density functional perturbation theory Á linear compact operator

Appendix Computational details
The moment expansion has been implemented in our development version of the CPMD [51] electronic structure package. The calculations have been performed using Density Functional . Averaged normalized overlap of the response density (of a water molecule) from the moment expansion and a reference DFPT calculation (due to the perturbing potential generated by the TIP3P water). The different colors refer to different number of moment-expanded states. We averaged over six density overlaps for each distance between the perturbing TIP3P water molecule and the responding water molecule. The six positions of the responding water molecules span an octahedron for each distance. Please note, that we compare the density differences (ffi response densities) due to a perturbation and not the resulting overall densities. [Color figure can be viewed at wileyonlinelibrary.com] Perturbation Theory [8,11,12,52,53] with Troullier-Martins [54] pseudo potentials in the Becke [55] and Lee-Yang-Parr [56] approximation for the exchange correlation kernel. The perturbing potential was expressed by a polynomial expansion. The basis was formed by the set of polynomials x l y m z n with l, m, n 2 N, which fulfill the condition l + m + n ≤ i max . The largest basis set consisted out of 54 polynomials (i max = 5). The calculation of the reference density response with CPMD was done for relaxed water geometry. We used a plane wave basis for all involved quantities, be it the electronic density, the perturbing potentials or the moment-expanded states. The real space representation naturally gives rise to a regular grid. For our chosen parameters, the grid increment is about 0.08 Angstroms, whereas the spread of the Gaussians in eq. (40) is chosen as σ = 0.5a 0 .