Density Fitting in Hybrid Functional Calculations

Computing the exact exchange energy in hybrid functional calculations can be accelerated using density fitting of two-electron integrals, in combination with adaptively-compressed exchange AdaptivelyCompressedExchange (which is enabled by default).

Octopus 17 implements density fitting for molecular systems using interpolative separable density fitting (ISDF), originally described in jcp.2015.09.014. Our implementation uses weighted kmeans clustering to identify reasonably optimal interpolation points on the real-space grid using the electron density as the weight function, and from these points constructs an equal number of interpolation vectors. We outline the fundaments of the theory in the following section, then go on to describe the options for using this feature in Octopus. For more details of the procedure, one can refer to jctc.7b00807.

The method is under active development, and is subject to change. ISDF is currently enabled for ground state calculations on molecular systems.

Theory

The product of two Kohn-Sham states defined at any point in space, $r$, can approximately be written by sampling their products on a small subset of interpolation points, ${r_\mu}$, and multiplying by a set of interpolation vectors ${\zeta_\mu(r)}$: $$ \begin{align} \phi_i(r),\psi_j(r) &\approx \sum_{\mu=1}^{N_\mu} \zeta_\mu(r), \phi_i\bigl( r_\mu\bigr), \psi_j\bigl( r_\mu\bigr). \end{align} $$ In matrix form this reads: $$ \begin{align} Z \approx \Theta,C, \end{align} $$ where $Z$ defines the matrix of product states (Z_{k,ij}=\phi_i(r_k)\psi_j(r_k)), $\Theta$ contains the interpolation vectors, and $C$ defines the fitting coefficients, which in this case are the product of Kohn-Shamn states defined at interpolation points, (C_{\mu,ij}=\phi_i(r_\mu)\psi_j(r_\mu)). Since the system is overdetermined, we can impose a Galerkin condition and solve in a least‐squares sense: $$ \begin{align} \Theta = Z,C^T,(C,C^T)^{-1}. \end{align} $$ Eq 3). has a separable form, with both (Z,C^T) and (C,C^T) factoring into products of quasi‐density matrices: $$ \begin{align} \bigl[Z,C^T\bigr]{r\mu} = P\phi(r, r_\mu);P_\psi(r, r_\mu), \end{align} $$ $$ \begin{align} \bigl[C,C^T\bigr]{\mu\nu} = P\phi(r_\mu, r_\nu);P_\psi(r_\mu, r_\nu), \end{align} $$ where the density matrices defined as: $$ \begin{align} P_\phi(r,r’) = \sum_{i=1}^m \phi_i(r),\phi_i(r’), ;; P_\psi(r,r’) = \sum_{j=1}^n \psi_j(r),\psi_j(r’). \end{align} $$ Finally, substituting this compressed form into the exact‐exchange operator gives: $$ \begin{align} \bigl(V_X[{\phi_i}]\psi_j\bigr)(r) ;\approx; -,\sum_{\mu=1}^{N_\mu} P_\phi(r,r_\mu); V_{\zeta_\mu}(r); \psi_j(r_\mu), \end{align} $$ where $V_{\zeta_\mu}(r)$ defines a potential: $$ \begin{align} V_{\zeta_\mu}(r) = \int \frac{\zeta_\mu(r’)}{|r - r’|} dr’. \end{align} $$ Importantly, evaluating (V_{\zeta_\mu}(r)) for each (\mu) requires solving one Poisson‐like equation per interpolation vector. In the absence of density fitting, one would need to solve a separate Poisson solve for every orbital pair ((i,j)), i.e. (O(N_e^2)) solves, whereas the ISDF scheme reduces this to (N_\mu) solves, with (N_\mu \sim N_e), greatly reducing the overall cost of computing the action of the exchange operator on the Kohn-Sham states.

Running with ISDF

When performing a hybrid functional calculation with ACE enabled, density fitting can also be enabled by setting ACEWithISDF = yes. The user must also specify the number of interpolation points (hence the number of interpolation vectors) to perform the density fitting with. This is controlled by ISDFNpoints and currently has no default, as it is system-specific. It’s important to choose a sufficiently large number of points/vectors to converge the exact exchange energy, however choosing too many results in an ill-conditioned Gram matrix of fitting coefficients, which hinders its inversion. In practice, the user does not need to worry about choosing too large a value, as Octopus internally always using the pseudo-inverse ( Moore–Penrose inverse) for ISDF. Details on choosing an optimal value for this variable are described below.

We can get a very good indication of the optimal number of interpolation points to use in ISDF by choosing an excessive number, and computing the rank of the fitting coefficient Gram matrix, $CC^T$. If the rank is equal to the number of points, this suggests one can choose more, or has the optimal number. If the rank is less than ISDFNpoints, then the rank defines the upper bound beyond which the exact exchange energy will not increase in precision (the same argument is not true of the two-electron integrals). Computing rank of the fitting coefficient Gram matrix can be done by setting CheckISDFNpoints = yes.

The user may also set the size of the subspace in which ACE is approximately exact, with ACESize = nocc. To be consistent, this also defines the number of Kohn-Sham states used in sum over states performed in ISDF. By default, this is set to take all states used in the calculation, however for this should be set to the number of occupied states when using ISDF. Choosing less states will affect the density fitting’s ability to capture the features of the occupied manifold, whereas choosing more will cause the Choleksy decomposition of the ACE algorithm to fail.

For ground state calculations, Octopus recomputes the interpolation points every SCF iteration. Whilst this might sound expensive, only the first step takes a noticeable amount of time. Changes in the density in subsequent SCF steps are relatively small, so performing a weighted kmeans calculation every SCF is not prohibitive.