Computing Second-Order Functional Derivatives with Respect to the External Potential

Following the increasing interest in the higher-order (functional) derivatives of conceptual density functional theory, we developed and implemented a method for calculating second-order functional derivatives with respect to the external potential. Our method is theoretically exact but involves two numerical approximations: the functional derivatives are expanded in a basis set, and the values of the corresponding expansion coefficients are determined by probing the molecular environment by a finite set of external potential perturbations. Exact solutions are obtained only in the limit of a complete basis set and an infinite number of distinct perturbations. We use this method to compute the atom-condensed linear response kernel for a series of six molecules and show that the results are comparable to the ones obtained by a previously proposed, approximate approach from second-order perturbation theory. The numerical error of the current implementation is about 0.01 au. Because the present method gives exact or quasi-exact solutions, it can be used as a benchmark against which approximate approaches are assessed.