Theor Chem Acc (2016) 135:6
1 3
derivative of the energy with respect to u ( r ) is ρ ( r ), the
electronic charge density at r . Therefore, χ can be written
as the functional derivative of the approximate electron
density (corresponding to the quantum chemical method
used) at a point r , ρ ( r ) with respect to u ( r ′), the potential at
another point r ′, at constant electron number N :
In this paper, we will concentrate on the static limit, but
generalization to the frequency-dependent case is straightforward at the Hartree–Fock and DFT levels, and we plan
to implement it.
Although the density response function is a fundamental quantity in theory, see, e.g., [ 6 – 12 ], there have been
few numerical studies of it, partly because it is a function
of 7 real variables (6 in the static limit), but also because
the diffi culty of its calculation. Yang [ 13 ], Kudrnovsky
and coworkers [ 14 ], Senet [ 15 ], and Ayers [ 16 ] discuss the
numerical evaluation of this quantity, but have apparently
not found a practical formulation, as none of these papers
has a single numerical result or example. The problem is
connected with the fact that evaluating χ ( r , r ′) at a single
pair of points ( r , r ′) involves a numerically awkward delta
function perturbation in the potential. Kudrnovsky and
coworkers [ 14 ] discuss a simplifi ed version of this function
in which the potential perturbation consists of varying the
nuclear charge of the atoms which can be evaluated simply using analytical derivative methods. This leads us to
naturally coarse-grained approximations to the DRF. The
foremost of these is the distributed multipole polarizability expansion of Stone and coworkers [ 3 , 17 ]. This method
partitions the DRF to local atomic multipole polarizabilities by using weight functions to defi ne the atomic regions.
Le Sueur and Stone [ 18 ] reported diffi culties in calculating physically meaningful distributed polarizabilities from
ab initio calculations, although these have been at least
partially resolved by the density fi tting method and its latest modifi cation [ 19 – 21 ]. Nevertheless, Rob and Szalewicz
[ 21 ] remark that the values obtained using the Misquitta–
Stone method [ 19 ] for distributed polarizabilities are still
not satisfactory, often containing large cancelling terms. A
further radical simplifi cation was introduced by Morita and
Kato [ 22 , 23 ]. In their technique, the charges and the electric potential are condensed to atoms, and thus polarization
effects are described as charge fl owing from one atom to
another. This method was applied very successfully by Hu
and Yang [ 24 ]. Condensing the charges to the nuclear positions greatly simplifi es both the calculations and the presentation, reducing χ from a function of 6 real variables to
an M × M matrix, where M is the number of atoms in the
system. However, its limits are clear. Consider, for instance,
(1)
χ(r, r
, 0) = χ(r, r
) =
δ 2 E
δu(r)δu(r )
N
=
δρ(r)
δu(r )
N
planar molecules like benzene and naphthalene (Table 1 ).
A polarizability model, expressed as derivatives of atomic
charges, gives obviously zero polarizability perpendicular
to the molecular plane. The calculated static dipole polarizability of benzene is 76.5 atomic units in-plane and 41.4 au
out-of-plane at the PW91/aug-cc-pVTZ level. The second
value is clearly not negligible compared to the fi rst, showing the importance of the atomic polarizability which is
neglected in the Morita–Kato–Yang scheme. The situation
is similar for naphthalene (Table 1 ). This limitation does
not apply to the distributed polarizability method.
The fi rst numerical results for the DRF we are aware of
were published by Savin et al. [ 25 ] for He and Be isoelectronic series. However, the method they used is restricted
to spherical electron densities and potentials, and thus its
applicability in chemistry is limited. A practical method,
based on standard coupled-perturbed Kohn–Sham, CPKS,
or its predecessor, the coupled-perturbed Hartree–Fock,
CPHF [ 26 , 27 ], was introduced recently by Yang et al. [ 28 ].
This paper presents a lucid rederivation of the CPHF/CPKS
equations and uses the fact that, in a fi nite basis set, the
only possible fi rst-order changes in the density are linear
combinations of φ i ( r ) φ a ( r ), where φ i and φ a are occupied
and virtual (vacant) orbitals, respectively. The results were
used extensively by Geerlings and coworkers who have
done the most signifi cant work in calculating the density
response function numerically at the single-determinant
level, and correlating the results with chemical concepts.
Their results, with a number of chemical applications, are
summarized in a recent review article [ 29 ] which also gives
a detailed exposition of such time-honored techniques
as the CPHF and CPKS methods. The method used [ 28 ]
avoids the problem with the delta functions by carrying
out the perturbations in molecular orbital (MO) basis and
transforming it to real space. In agreement with the number of possible fi rst-order density variations, this procedure
leads to nV distinct perturbations, where n is the number
of occupied and V is the number of virtual orbitals. Each
perturbation requires the solution of the coupled-perturbed
Hartree–Fock or Kohn–Sham (CPHF and CPKS) equations, giving 34 × 610 = 20740 CPHF equations in the
Table 1 Electric dipole polarizabilities of benzene and naphthalene
in the coupled and uncoupled Kohn–Sham perturbation theory
Geometry optimized using the BP86/cc-pVTZ level. Aug-cc-pVTZ
basis set
Molecule
Method
α ( xx )
α ( yy )
α ( zz )
Benzene
CPKS
76.5
76.5
41.4
Benzene
UCKS
134.9
134.9
67.5
Naphthalene
CPKS
163.3
117.9
61.2
Naphthalene
UCKS
293.4
224.3
105.5
10
Reprinted from the journal
1 3
derivative of the energy with respect to u ( r ) is ρ ( r ), the
electronic charge density at r . Therefore, χ can be written
as the functional derivative of the approximate electron
density (corresponding to the quantum chemical method
used) at a point r , ρ ( r ) with respect to u ( r ′), the potential at
another point r ′, at constant electron number N :
In this paper, we will concentrate on the static limit, but
generalization to the frequency-dependent case is straightforward at the Hartree–Fock and DFT levels, and we plan
to implement it.
Although the density response function is a fundamental quantity in theory, see, e.g., [ 6 – 12 ], there have been
few numerical studies of it, partly because it is a function
of 7 real variables (6 in the static limit), but also because
the diffi culty of its calculation. Yang [ 13 ], Kudrnovsky
and coworkers [ 14 ], Senet [ 15 ], and Ayers [ 16 ] discuss the
numerical evaluation of this quantity, but have apparently
not found a practical formulation, as none of these papers
has a single numerical result or example. The problem is
connected with the fact that evaluating χ ( r , r ′) at a single
pair of points ( r , r ′) involves a numerically awkward delta
function perturbation in the potential. Kudrnovsky and
coworkers [ 14 ] discuss a simplifi ed version of this function
in which the potential perturbation consists of varying the
nuclear charge of the atoms which can be evaluated simply using analytical derivative methods. This leads us to
naturally coarse-grained approximations to the DRF. The
foremost of these is the distributed multipole polarizability expansion of Stone and coworkers [ 3 , 17 ]. This method
partitions the DRF to local atomic multipole polarizabilities by using weight functions to defi ne the atomic regions.
Le Sueur and Stone [ 18 ] reported diffi culties in calculating physically meaningful distributed polarizabilities from
ab initio calculations, although these have been at least
partially resolved by the density fi tting method and its latest modifi cation [ 19 – 21 ]. Nevertheless, Rob and Szalewicz
[ 21 ] remark that the values obtained using the Misquitta–
Stone method [ 19 ] for distributed polarizabilities are still
not satisfactory, often containing large cancelling terms. A
further radical simplifi cation was introduced by Morita and
Kato [ 22 , 23 ]. In their technique, the charges and the electric potential are condensed to atoms, and thus polarization
effects are described as charge fl owing from one atom to
another. This method was applied very successfully by Hu
and Yang [ 24 ]. Condensing the charges to the nuclear positions greatly simplifi es both the calculations and the presentation, reducing χ from a function of 6 real variables to
an M × M matrix, where M is the number of atoms in the
system. However, its limits are clear. Consider, for instance,
(1)
χ(r, r
, 0) = χ(r, r
) =
δ 2 E
δu(r)δu(r )
N
=
δρ(r)
δu(r )
N
planar molecules like benzene and naphthalene (Table 1 ).
A polarizability model, expressed as derivatives of atomic
charges, gives obviously zero polarizability perpendicular
to the molecular plane. The calculated static dipole polarizability of benzene is 76.5 atomic units in-plane and 41.4 au
out-of-plane at the PW91/aug-cc-pVTZ level. The second
value is clearly not negligible compared to the fi rst, showing the importance of the atomic polarizability which is
neglected in the Morita–Kato–Yang scheme. The situation
is similar for naphthalene (Table 1 ). This limitation does
not apply to the distributed polarizability method.
The fi rst numerical results for the DRF we are aware of
were published by Savin et al. [ 25 ] for He and Be isoelectronic series. However, the method they used is restricted
to spherical electron densities and potentials, and thus its
applicability in chemistry is limited. A practical method,
based on standard coupled-perturbed Kohn–Sham, CPKS,
or its predecessor, the coupled-perturbed Hartree–Fock,
CPHF [ 26 , 27 ], was introduced recently by Yang et al. [ 28 ].
This paper presents a lucid rederivation of the CPHF/CPKS
equations and uses the fact that, in a fi nite basis set, the
only possible fi rst-order changes in the density are linear
combinations of φ i ( r ) φ a ( r ), where φ i and φ a are occupied
and virtual (vacant) orbitals, respectively. The results were
used extensively by Geerlings and coworkers who have
done the most signifi cant work in calculating the density
response function numerically at the single-determinant
level, and correlating the results with chemical concepts.
Their results, with a number of chemical applications, are
summarized in a recent review article [ 29 ] which also gives
a detailed exposition of such time-honored techniques
as the CPHF and CPKS methods. The method used [ 28 ]
avoids the problem with the delta functions by carrying
out the perturbations in molecular orbital (MO) basis and
transforming it to real space. In agreement with the number of possible fi rst-order density variations, this procedure
leads to nV distinct perturbations, where n is the number
of occupied and V is the number of virtual orbitals. Each
perturbation requires the solution of the coupled-perturbed
Hartree–Fock or Kohn–Sham (CPHF and CPKS) equations, giving 34 × 610 = 20740 CPHF equations in the
Table 1 Electric dipole polarizabilities of benzene and naphthalene
in the coupled and uncoupled Kohn–Sham perturbation theory
Geometry optimized using the BP86/cc-pVTZ level. Aug-cc-pVTZ
basis set
Molecule
Method
α ( xx )
α ( yy )
α ( zz )
Benzene
CPKS
76.5
76.5
41.4
Benzene
UCKS
134.9
134.9
67.5
Naphthalene
CPKS
163.3
117.9
61.2
Naphthalene
UCKS
293.4
224.3
105.5
10
Reprinted from the journal