# Optimization of constrained density functional theory

David D. O’Regan<EMAIL_ADDRESS>School of Physics, CRANN and AMBER,
Trinity College Dublin, Dublin 2, Ireland  Gilberto Teobaldi
<EMAIL_ADDRESS>Stephenson Institute for Renewable Energy and
Department of Chemistry, The University of Liverpool, L69 3BX Liverpool,
United Kingdom  Beijing Computational Science Research Center, Beijing 100094,
China

###### Abstract

Constrained density functional theory (cDFT) is a versatile electronic
structure method that enables ground-state calculations to be performed
subject to physical constraints. It thereby broadens their applicability and
utility. Automated Lagrange multiplier optimisation is necessary for multiple
constraints to be applied efficiently in cDFT, for it to be used in tandem
with geometry optimization, or with molecular dynamics. In order to facilitate
this, we comprehensively develop the connection between cDFT energy
derivatives and response functions, providing a rigorous assessment of the
uniqueness and character of cDFT stationary points while accounting for
electronic interactions and screening. In particular, we provide a new, non-
perturbative proof that stable stationary points of linear density constraints
occur only at energy maxima with respect to their Lagrange multipliers. We
show that multiple solutions, hysteresis, and energy discontinuities may occur
in cDFT. Expressions are derived, in terms of convenient by-products of cDFT
optimization, for quantities such as the dielectric function and a condition
number quantifying ill-definition in multi-constraint cDFT.

###### pacs:

71.15.-m, 31.15.E-, 71.15.Qe, 71.15.Dx

##  I Background

Constrained density functional theory (cDFT) Kaduk  _et al._ (2012) is a
generalization of density functional theory (DFT) Hohenberg and Kohn (1964);
Kohn and Sham (1965) in which external constraints are applied in order to
simulate excitation processes, to calculate response properties, or to impose
a physical condition that is not met by the unconstrained approximate
exchange-correlation functional. Such constraints may be applied to
expectation values of the charge or spin density, their sums, differences, and
moments within pre-defined spatial regions Wu and Van Voorhis (2006a, 2005);
Segal  _et al._ (2007); Kowalczyk  _et al._ (2010); Dederichs  _et al._
(1984); Wu and Van Voorhis (2006b); Oberhofer and Blumberger (2009a); Sit  _et
al._ (2006); Tim Kowalczyk  _et al._ (2011); Behler  _et al._ (2008); Rudra
_et al._ (2006, 2007); Schmidt  _et al._ (2008). They are necessarily chosen
on the basis of physical intuition and experience. Constraining potentials
that are non-local or orbital-dependent may also be introduced, moving beyond
formal DFT Wu and Van Voorhis (2006a). cDFT enables individual excited states
to be studied within the well-established framework of ground-state DFT Wu and
Van Voorhis (2005, 2006a); Segal  _et al._ (2007); Difley  _et al._ (2008);
Kowalczyk  _et al._ (2010); Sena  _et al._ (2011); Turban  _et al._ (2016),
particularly those excited states which may be represented as the ground-state
for some potential. While these excitations are not guaranteed to match the
neutral excitations of the system, yielded by time-dependent DFT Runge and
Gross (1984), for example, their description may nonetheless benefit from
physical conditions, such as charge transfer, which may be absent from the
approximate functional but reintroduced using cDFT. As such, cDFT offers
important insights that are challenging to obtain otherwise Wu and Van Voorhis
(2005); Difley  _et al._ (2008); Difley and Van Voorhis (2011); Nakamura  _et
al._ (2008); Sit  _et al._ (2006).

In practice, cDFT has proven to be a very efficient approach for simulating
neutral excitations in molecular systems, particularly in cases where a clear
spatial delineation may be made between charge (or spin) donor and acceptor
regions Wu and Van Voorhis (2005, 2006a); Difley  _et al._ (2008); Segal  _et
al._ (2007); Kowalczyk  _et al._ (2010); Sena  _et al._ (2011); Roychoudhury
_et al._ (2016); Souza  _et al._ (2013); Turban  _et al._ (2016); Difley and
Van Voorhis (2011); Nakamura  _et al._ (2008). cDFT is a significant asset,
therefore, to the simulation of exciton formation, where the incorrect long-
ranged behaviour of conventional local or semi-local exchange-correlation
functionals may be partially corrected by using appropriately constructed
constraints Sau  _et al._ (2008); Souza  _et al._ (2013); Roychoudhury  _et
al._ (2016); Brooke  _et al._ (2015); Turban  _et al._ (2016). It has also
been used to calculate electron transfer Sit  _et al._ (2006); Wu and Van
Voorhis (2006c); Ding  _et al._ (2010); Oberhofer and Blumberger (2009b,
2010a); Kubas  _et al._ (2014, 2015); Oberhofer and Blumberger (2010b, 2012);
McKenna and Blumberger (2012); Blumberger and McKenna (2013), excitation
energy transfer Yeganeh and Van Voorhis (2010); Yost  _et al._ (2014), and
exchange coupling parameters Rudra  _et al._ (2006, 2007); Evans  _et al._
(2008) for use in model Hamiltonians, as well as Coulomb interaction
parameters for methods such as DFT+$U$ Anisimov  _et al._ (1991); Meider and
Springborg (1998); Danielsen (1986); Meider and Springborg (1999); Hybertsen
_et al._ (1989); Petukhov  _et al._ (2003). Moreover, cDFT has been shown to
provide an effective correction for the self-interaction error exhibited by
approximate functionals when calculating diabatic free-energy surfaces for
electron-transfer reactions Sit  _et al._ (2006). As a promising antidote to
static correlation error in approximate functionals, cDFT has been used to
generate small, efficient basis sets for configuration interaction
calculations, by enabling the most relevant charge and spin states to be
straightforwardly sampled Wu  _et al._ (2007, 2009); Nakamura  _et al._
(2008); Kaduk and Van Voorhis (2010); Kaduk  _et al._ (2014). For a recent
comprehensive review of cDFT, we refer the reader to Ref. Kaduk  _et al._ ,
2012.

cDFT may yet play unforeseen roles in future first-principles atomistic
simulation. As the field moves increasingly towards the automated construction
and interrogation of materials databases generated using high-throughput DFT
approaches Jain  _et al._ (2011); Curtarolo  _et al._ (2013), for example, it
could be used in the large-scale screening of candidate charge-transfer and
energy-transfer materials or, as we describe below, to screen for the average
local microscopic dielectric functions of complex materials and interfaces. In
order for the great utility and potential of the cDFT approach to be fully and
routinely realized in the simulation of charge-transfer excitations, and in
degenerate or strongly-interacting systems, it must be efficiently automated,
reliable, and convenient for users. For this, robust optimisation algorithms
for the Lagrange-multipliers enforcing the constraint functionals of cDFT are
desirable, and indeed necessary in cases of multiple simultaneous constraints
being applied, such as on charge and spin Kaduk  _et al._ (2012); Wu and Van
Voorhis (2006b, a); Oberhofer and Blumberger (2009a); Sena  _et al._ (2011).
Additionally, automated Lagrange multiplier updates at each ionic
configuration are indispensable when performing geometry optimisation Wu and
Van Voorhis (2005, 2006a); Segal  _et al._ (2007); Difley  _et al._ (2008);
Kowalczyk  _et al._ (2010); Sena  _et al._ (2011) or molecular dynamics Sit
_et al._ (2006); Oberhofer and Blumberger (2010a, 2009a, b); R̆ezác̆  _et al._
(2012) in tandem with cDFT.

Critical to both the theoretical underpinning and viability of cDFT
optimization is the nature of the energy landscape with respect to the
Lagrange multipliers that determine the strength of its constraining
potentials. In particular, certainty about the nature and uniqueness of any
stationary points at which the constraints are satisfied is a prerequisite to
efficiently locating them numerically. Wu and Van Voorhis (W&VV) Wu and Van
Voorhis (2005) carried out the pioneering and enabling work in this area,
analysing the relevant derivatives, and their principal results have been
subsequently synopsized in numerous works Kaduk  _et al._ (2012); Wu and Van
Voorhis (2006a); Souza  _et al._ (2013); Sena  _et al._ (2011); Ratcliff  _et
al._ (2015). It was concluded by W&VV Wu and Van Voorhis (2005) on the basis
of non-degenerate perturbation theory that a non-trivial stationary point, for
an arbitrary constraint on the electron density, arises only at a maximum of
the total-energy with respect to a cDFT Lagrange multiplier, and that this
solution is unique. This is a central result in cDFT, suggesting the
feasibility of its routine automated optimization, which has been extended to
multivariate cases in Refs. Wu and Van Voorhis, 2006a and Kaduk  _et al._ ,
2012.

##  II Introduction and Motivation

In this work, we rigorously generalize the latter result, building upon the
foundation provided by W&VV’s cDFT stationary point classification, first
showing that the analysis becomes inconclusive when electronic screening
effects are considered. Specifically, we find that while the cDFT energy
curvature 111The “curvature” is used here as a convenient shorthand for the
second derivative. We do not imply the geometric curvature, which equals the
second derivative only at stationary points. formula derived by W&VV is
appropriate for updating cDFT Lagrange multipliers during the density update
step, or inner loop, of self-consistent field DFT algorithms Wu (2016); Wu and
Yang (2003); Wu and Van Voorhis (2006a), it is not applicable to the self-
consistently relaxed total-energy relevant to the global classification of
cDFT solutions. Fig. 3 illustrates the large discrepancy between the cDFT
energy curvatures calculated using W&VV’s formula (solid circles) and those
evaluated using finite-differences (open squares), hence the necessity to
revisit the topic here.

Figure 1: (Color online) Multiple solutions and hysteresis in spin-
constrained, open-shell, stretched molecular hydrogen. Shown are (top) the
magnetic moment difference $\Delta M$ between constrained regions versus the
Lagrange multiplier $V_{\textit{c}}$, (middle) the cDFT total-energy $W$
versus $V_{\textit{c}}$ corresponding to a target moment difference $\Delta
M_{\textit{c}}=0~{}\mu_{B}$, and (bottom) the DFT energy component
$E_{\textrm{DFT}}$ versus $\Delta M$. “OS(-)” indicates data obtained using
the open-shell ground-state with negative $\Delta M$ at
$\left(V_{\textit{c}}=0\right)$, “OS(+)” the symmetry-related open-shell
solution with positive $\Delta M\left(V_{\textit{c}}=0\right)$, with data
adapted from the latter, and “CS” the meta-stable closed-shell solution.

In addition, Fig. 1 shows cDFT data that exhibits effects not hitherto
discussed in the relevant literature, to our knowledge, namely multiple
solutions, hysteresis, and energy discontinuities, and thus further motivates
the present study. Here, the ONETEP linear-scaling simulation code Skylaris
_et al._ (2005), as discussed in Section IV, was used to carry out cDFT
calculations on the hydrogen molecule stretched to an internuclear distance of
$3.2~{}\mathrm{a}_{0}$, which, using the PBE functional Perdew  _et al._
(1996) with no spin-orbit coupling, lies just beyond the Coulson-Fischer point
at which an open-shell singlet ground-state becomes favoured Coulson and
Fischer (1949); Ruzsinszky  _et al._ (2005). A constraint was placed on the
difference of spin magnetic moments, $\Delta M$, between the two hydrogen
atoms, defined on the basis of their isolated $1s$ valence pseudo-orbitals.
The constraint target was set to $\Delta M_{\textit{c}}=0~{}\mu_{B}$, and the
Lagrange multiplier $V_{\textit{c}}$ was defined such that increasing its
value increased the spin-dependent potential acting to decrease $\Delta M$. It
was found that the unpolarized, closed shell “CS” state is meta-stable,
starting from which any non-zero value of $V_{\textit{c}}$ initiates a
collapse to one of two symmetry-related, degenerate open-shell ground-states,
“OS(-)” and “OS(+)” defined in the caption of Fig. 1.

The constraint proved capable of traversing between the two energy basins
associated with these states (the “CS” state lies at the saddle point
connecting them), insofar as that it caused abrupt switching between states at
a critical constraining potential of $V_{\textit{c}}\approx\pm 81.1$ meV. The
top panel of Fig. 1 shows that multiple energy minima are possible for
particular values of the cDFT Lagrange multiplier $V_{\textit{c}}$, that
certain ranges of the constraint target (in this case $\Delta M_{\textit{c}}$)
may be inaccessible to cDFT, and that its response diverges simultaneously
with an energy discontinuity at the transition between states. The middle
panel shows that the total-energy $W$ may exhibit non-differentiable cusps at
such points, close to which cDFT optimisation is impracticable. The bottom
panel shows the DFT contribution to the total-energy, $E_{\textrm{DFT}}$, as a
function of the constrained quantity, and the two symmetry-related “OS” energy
basins in question. Geometrically, it is inevitable that at least some
solutions within the inaccessible interval connecting these two curves would,
in principle, exhibit negative curvatures. We will show that such solutions,
unless meta-stable like the “CS” state, are unstable due to anomalous response
function sign, and thus cannot be realized using cDFT.

Prompted by these results, in the following we provide generalized energy
curvature expressions which ensure that the stable stationary points of non-
trivial linear constraints in the density may occur only at maxima of the
total-energy with respect to their Lagrange multipliers, thereby cementing the
theoretical basis of automated cDFT optimization. Our approach is general up
to arbitrary orders in response and it also lifts the assumption of orbital
non-degeneracy made in W&VV’s treatment. It also allows for the maximizability
of the cDFT energy to imply the uniqueness of such maxima only when the
unconstrained system is devoid of multiple electronic minima, a specific
counterexample of which is demonstrated by Fig 1. The consistency of our
analytical approach is validated throughout this work by means of numerically
verified equalities, in some cases newly established, between integrated
linear-response functions and components of the cDFT total-energy curvature.
As such, the present work represents a timely, comprehensive treatment of
density response from the energy landscape perspective. We provide energy-
curvature relationships for each of the integrated response and inverse-
response functions, both interacting and non-interacting, which we generalize
to multiple constraints.

These response functions are by-products of cDFT Lagrange multiplier
optimization, and they can be used, for example, to calculate the average
dielectric constant in a particular region, which is of critical importance
for supercell convergence acceleration Turban  _et al._ (2016) and implicit
solvation Dziedzic  _et al._ (2011); Andreussi  _et al._ (2012) schemes. A
formula for the cDFT optimization condition number is provided, intended for
identifying systems in which extremization of the cDFT total-energy is to be
expected to present challenges, and in which Newton’s method may be apposite,
as has previously been suggested within the cDFT context in Refs. Wu and Van
Voorhis, 2006a; Oberhofer and Blumberger, 2009a; Kaduk  _et al._ , 2012.

Given the broad span of the physical and mathematical issues necessarily
considered and revised, and to favor the readability of the article, it is
organized into Sections, as follows. In Section III, we present and discuss
the previously overlooked role of dielectric screening in the energy
derivatives in self-consistent cDFT. In Section IV, we introduce the formal
connections between the Lagrange multiplier curvature of the different
contributions to the cDFT total-energy and the integrated electronic response
function, both in the interacting and non-interacting cases. These newly
derived formulae are then applied to globally characterize the stationary
points of cDFT, both in the single constraint (Section V) and multiple
constraint (Section VI) cases. In Section VII, we provide a synopsis of our
main findings and conclusions.

##  III Dielectric screening of Constrained DFT energy derivatives

We begin our analysis with the cDFT Kaduk  _et al._ (2012) constraint
functional as per the definition and notation introduced in the founding
article on cDFT automation by W&VV Wu and Van Voorhis (2005). We consider an
electronic system treated using Kohn-Sham density-functional theory (DFT)
Hohenberg and Kohn (1964); Kohn and Sham (1965), subject to an arbitrary
constraint on its electron density (see W&VV’s Eq. 1, Ref. Dederichs  _et al._
, 1984 and footnote 222A single, strictly local occupancy constraint is
considered in Ref. Wu and Van Voorhis, 2005. We retain these restrictions so
as not to obscure the fundamental aspects under consideration. These
conditions are typically lifted in practical cDFT calculations, bringing us
into multivariate optimisation of constrained Kohn-Sham spin-density
functional theory, which may also be non-local or orbital dependent.) of the
linear form

$\displaystyle
C\left[\rho\right]=N\left[\rho\right]-N_{\textit{c}};\;\;N\left[\rho\right]=\sum_{\sigma}\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\rho^{\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}.$
(1)

Here, $\rho^{\sigma}\left(\mathbf{r}\right)$ is the electronic density of spin
$\sigma$, $w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)$ is an arbitrary
local weight function describing a spatial region of particular interest in
the system, and $N_{\textit{c}}$ is the target electron number to be enforced
on that region. In order to apply this constraint, a term is added to the
conventional DFT total-energy, $E_{\textrm{DFT}}\left[\rho\right]$, to build
the functional given by (c.f. W&VV’s Eq. 4)

$\displaystyle W\left[\rho,V_{\textit{c}}\right]={}$ $\displaystyle
E_{\textit{c}}\left[\rho,V_{\textit{c}}\right]+E_{\textrm{DFT}}\left[\rho\right],\;\;\mbox{with}$
(2) $\displaystyle E_{\textit{c}}\left[\rho,V_{\textit{c}}\right]={}$
$\displaystyle V_{\textit{c}}C\left[\rho\right],$ (3)

where $V_{\textit{c}}$ is the Lagrange multiplier. Minimizing $W$ with respect
to the density via the Kohn-Sham orbitals $\phi_{i\sigma}$, for a given
$V_{\textit{c}}$, under the condition that these orbitals are orthonormalized
for each spin, gives rise to the Kohn-Sham equations Kohn and Sham (1965)
including a constraining potential
$V_{\textit{c}}w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)$. This
minimization, which is equivalent to solving the constrained Kohn-Sham
equations, does not correspond to the free extremization $\delta
W/\delta\phi^{\ast}_{i\sigma}=0$ invoked by W&VV, since the functional
derivative cannot encode the orbital orthonormality constraint. Rather, it
instead corresponds to extremizing the Lagrangian

$\displaystyle\Omega\left[\rho,V_{\textit{c}}\right]={}$ $\displaystyle
W\left[\rho,V_{\textit{c}}\right]$ (4)
$\displaystyle-\sum_{\sigma}\sum_{ij}^{N_{\sigma}}\varepsilon_{ij\sigma}\left(\int\phi_{i\sigma}^{\ast}\left(\mathbf{r}\right)\phi_{j\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}-\delta_{ij}\right).$

Here, we have assumed that the system is Kohn-Sham insulating, for simplicity,
with $N^{\sigma}$ electrons per spin $\sigma$. Following application of the
condition $\delta\Omega/\delta\phi^{\ast}_{i\sigma}=0$, we may perform a
unitary transformation among the resulting equations. This also transforms the
orbitals, yet it presents no difficulties since
$\Omega\left[\rho,V_{\textit{c}}\right]$ is invariant under such
transformations for density functionals. Diagonalising the matrix Lagrange
multiplier $\varepsilon_{ij\sigma}$, thereby, returns the Kohn-Sham cDFT
equations of W&VV’s Eq. 5, with eigenvalues $\varepsilon_{i\sigma}$. Thus, we
may succinctly write, at the physically relevant minimum of W, the expression
$\delta W/\delta\phi^{\ast}_{i\sigma}=\hat{H}_{\sigma}\phi_{i\sigma}$ or
indeed its complex conjugate $\delta
W/\delta\phi_{i\sigma}=\phi^{\ast}_{i\sigma}\hat{H}_{\sigma}$, where
$\hat{H}_{\sigma}$ is the Kohn-Sham cDFT Hamiltonian for spin $\sigma$.
Following W&VV, we may next define the function $W\left(V_{\textit{c}}\right)$
as the evaluation of $W\left[\rho,V_{\textit{c}}\right]$ using the density
generated by the orthonormal orbitals which solve the Kohn-Sham equations
including the constraining potential
$V_{\textit{c}}w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)$ or, for the
avoidance of doubt concerning Kohn-Sham excited states, as the physical
minimum of the total-energy for a given $V_{\textit{c}}$.

###  III.1 Total-energy first derivative

We now begin to analyze the derivatives of $W\left(V_{\textit{c}}\right)$
required for the location of cDFT solutions. Partial derivatives couple only
explicit dependencies, and are sufficient for optimising the Lagrange
multiplier during the density update step, or inner loop, of self-consistent
field DFT algorithms Wu and Yang (2003); Wu and Van Voorhis (2006a). In order
to determine the character of cDFT stationary points globally, on the other
hand, we must consider total derivatives, which include orbital and density
relaxation effects. In simulations where a number of constraints are
simultaneously applied, their Lagrange multipliers are independent variables,
so that the Hessian of interest is the matrix of mixed second total
derivatives. These simulations are discussed in Section VI.

The first total derivative of the cDFT total-energy with respect to the
Lagrange multiplier, $V_{\textit{c}}$, is given by

$\displaystyle\frac{dW}{dV_{\textit{c}}}={}$
$\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\frac{\delta
W}{\delta\phi_{i\sigma}^{\ast}}\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}+\textit{c.c.}\right]+\frac{\partial
W}{\partial V_{\textit{c}}}$ $\displaystyle={}$
$\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\left(\hat{H}_{\sigma}\phi_{i\sigma}\right)\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}+\textit{c.c.}\right]$
$\displaystyle+\left(\sum_{\sigma}\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\rho^{\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}-N_{\textit{c}}\right),$
(5)

where the trace symbol Tr denotes an integral over space since the operators
are all local, and c.c. represents the complex conjugate of the preceding
term. For any value of $V_{\textit{c}}$, the orbitals that generate the
density minimizing $\Omega\left[\rho,V_{\textit{c}}\right]$ are unique up to
unitary transformations, and we are free to choose the set that diagonalize
the Hamiltonian. This allows us to simplify the latter expression, since it
guarantees that the orbital-coupling term

$\displaystyle\textrm{Tr}\left[\left(\hat{H}_{\sigma}\phi_{i\sigma}\right)\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}\right]={}$
$\displaystyle\varepsilon_{i\sigma}\textrm{Tr}\left[\phi_{i\sigma}\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}\right]$
(6) $\displaystyle={}$
$\displaystyle\varepsilon_{i\sigma}\int\phi_{i\sigma}\left(\mathbf{r}\right)\sum_{a\neq
i}\phi_{a\sigma}^{\ast}\left(\mathbf{r}\right)\;d\mathbf{r}$
$\displaystyle\times\int\frac{\phi_{a\sigma}^{\ast}\left(\mathbf{r^{\prime}}\right)\frac{dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r^{\prime}}\right)}{dV_{\textit{c}}}\phi_{i\sigma}\left(\mathbf{r^{\prime}}\right)}{\varepsilon_{i\sigma}-\varepsilon_{a\sigma}}\;d\mathbf{r^{\prime}}$

evaluates to zero by virtue of the orthonormality of $\phi_{i\sigma}$ and
$\phi_{a\sigma}$ for $a\neq i$. Here, $v_{\sigma}^{\textrm{KS}}$ is the Kohn-
Sham potential, i.e., the total effective potential which enters density-
functional perturbation theory Baroni  _et al._ (1987, 2001).

If focusing on the self-consistent field DFT inner loop as per W&VV, we may
neglect screening effects and, thereby, assert that the change in total
potential equals the external perturbation
$\delta\hat{v}^{\textrm{external}}_{\sigma}$, and then
$\delta\hat{v}_{\sigma}^{\textrm{KS}}=\delta\hat{v}_{\sigma}^{\textrm{external}}=\hat{w}^{\sigma}_{\textit{c}}\delta
V_{\textit{c}}$, whence
$dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)/dV_{\textit{c}}=\hat{w}^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)$
in the above expression. More generally, however, an account of electronic
screening of the perturbation is necessary, and such effects are encapsulated
in the inverse microscopic dielectric function defined by
$\epsilon^{-1}_{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)=dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)/dv_{\sigma^{\prime}}^{\textrm{external}}\left(\mathbf{r^{\prime}}\right)$
333We use a lunate epsilon $\epsilon$ for the microscopic dielectric function
in order to distinguish from it from the Kohn-Sham eigenvalues $\varepsilon$.
For pure (density-constrained rather than non-local) cDFT, we may write that

$\displaystyle\frac{dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)}{dV_{\textit{c}}}={}$
$\displaystyle\sum_{\sigma^{\prime}}\int\frac{dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)}{dv_{\sigma^{\prime}}^{\textrm{external}}\left(\mathbf{r^{\prime}}\right)}\frac{dv_{\sigma^{\prime}}^{\textrm{external}}\left(\mathbf{r^{\prime}}\right)}{dV_{\textit{c}}}\;d\mathbf{r^{\prime}}$
$\displaystyle={}$
$\displaystyle\sum_{\sigma^{\prime}}\int\epsilon^{-1}_{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)w^{\sigma^{\prime}}_{\textit{c}}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r^{\prime}}$
$\displaystyle\equiv{}$
$\displaystyle\left(\epsilon^{-1}w_{\textit{c}}\right)^{\sigma}\left(\mathbf{r}\right).$
(7)

For cDFT with non-local potentials, the symmetric form
$\left(\epsilon^{-1/2}w_{\textit{c}}\epsilon^{-1/2}\right)^{\sigma}\left(\mathbf{r},\mathbf{r^{\prime}}\right)$
may be used in place of the latter in order to ensure that the potential
remains Hermitian.

Figure 2: (Color online) The cDFT total-energy $W$, and its constraint,
$E_{\textit{c}}$, and DFT, $E_{\textrm{DFT}}$, components, as a function of
the Lagrange multiplier $V_{\textit{c}}$, for a charge-constrained nitrogen
molecule. The $V_{\textit{c}}$ values for the data points at which
$E_{\textit{c}}$ and $W$ attain maxima are shown with dashed vertical red and
green lines, respectively. Inset: the left-hand atom is constrained to lose
charge with respect to its ground-state population, using an on-atom
population analysis combining $s$ and $p$ orbitals. With respect to the
ground-state density, charge depleted regions are shown with a cyan charge-
difference isosurface and the charge augmented region is shown with an orange
isosurface. The unconstrained right-hand atom exhibits strong polarization,
and depletion in the lone-pair region.

Screening effects notwithstanding, all cancels to zero in the total derivative
of $W\left(V_{\textit{c}}\right)$ except for the explicit constraint
contribution on the final line of Eq. III.1, namely $\partial W/\partial
V_{\textit{c}}=C$, which then evaluates to zero when the constraint is
satisfied. Thus, the task of enforcing the constraint condition is transformed
into that of locating the stationary points of $W$ with respect to
$V_{\textit{c}}$. As an aside, this result mirrors W&VV’s Eq. 6, which was
derived in a slightly different way by invoking $\delta
W/\delta\phi^{\ast}_{i\sigma}=0$ rather than $\delta
W/\delta\phi^{\ast}_{i\sigma}=\hat{H}_{\sigma}\phi_{i\sigma}$, the numerical
distinction between which is vanishing due to Kohn-Sham orbital
orthonormality. Finally, we note that the vanishing trace in Eq. III.1 may be
partitioned into two contributions which cancel, for any value of
$V_{\textit{c}}$, by means of the expression

$\displaystyle\frac{\delta W}{\delta\phi_{i\sigma}^{\ast}}={}$
$\displaystyle\frac{\delta
E_{\textit{c}}}{\delta\phi_{i\sigma}^{\ast}}+\frac{\delta
E_{\textrm{DFT}}}{\delta\phi_{i\sigma}^{\ast}}=V_{\textit{c}}\frac{\delta
C}{\delta\phi_{i\sigma}^{\ast}}+\frac{\delta
E_{\textrm{DFT}}}{\delta\phi_{i\sigma}^{\ast}}$ $\displaystyle\Rightarrow\quad
0={}$ $\displaystyle V_{\textit{c}}\textrm{Tr}\left[\frac{\delta
C}{\delta\phi_{i\sigma}^{\ast}}\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}\right]+\textrm{Tr}\left[\frac{\delta
E_{\textrm{DFT}}}{\delta\phi_{i\sigma}^{\ast}}\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}\right].$
(8)

Thus, both the DFT energy $E_{\textrm{DFT}}$ and the constraint energy
$E_{\textit{c}}$ may individually contribute substantially to the derivatives
of $W\left(V_{\textit{c}}\right)$, as shown in Fig. 2. As a result, a
stationary total-energy $W$ with respect to the cDFT Lagrange multiplier does
not imply a stationary constraint contribution $E_{\textit{c}}$ alone, and
vice versa.

###  III.2 Total-energy second derivative by means of non-degenerate
perturbation theory

The second derivative or “curvature” of $W\left(V_{\textit{c}}\right)$, is
required to classify any stationary point, or points, at which the constraint
is satisfied. In cases where the second derivative vanishes, higher
derivatives may also be needed. We consider here the self-consistent cDFT
energy landscape, rather than the unscreened problem specific to the inner
loop of self-consistent field DFT codes Wu (2016); Wu and Yang (2003); Wu and
Van Voorhis (2006a), and the resulting curvature differs from that of W&VV’s
treatment in magnitude and, potentially, in sign. Following from Eq. III.1 and
applying the product rule for differentiation where necessary, we may write
that

$\displaystyle\frac{d^{2}W}{dV_{\textit{c}}^{2}}={}$
$\displaystyle\frac{d}{dV_{\textit{c}}}\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\left(\hat{H}_{\sigma}\phi_{i\sigma}\right)\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}+\textit{c.c.}\right]+\frac{dC}{dV_{\textit{c}}}$
$\displaystyle={}$
$\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\left(\frac{d\hat{H}_{\sigma}}{dV_{\textit{c}}}\phi_{i\sigma}\right)\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}+\textit{c.c.}\right]$
$\displaystyle+\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\left(\hat{H}_{\sigma}\frac{d\phi_{i\sigma}}{dV_{\textit{c}}}\right)\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}+\textit{c.c.}\right]$
$\displaystyle\quad\;\;+\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\left(\hat{H}_{\sigma}\phi_{i\sigma}\right)\frac{d^{2}\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}^{2}}+\textit{c.c.}\right]$
$\displaystyle\quad\;\;\quad\;\;+\sum_{\sigma}\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\frac{d\rho^{\sigma}\left(\mathbf{r}\right)}{dV_{\textit{c}}}\;d\mathbf{r}.$
(9)

This makes the contributions arising at second order in perturbation theory
explicit. These rather cumbersome terms may be circumvented by noting that the
eigencondition
$\hat{H}_{\sigma}\phi_{i\sigma}=\varepsilon_{i\sigma}\phi_{i\sigma}$ holds
continuously as we vary the parameter $V_{\textit{c}}$. As a result, the
quantity expressed in Eq. 6 vanishes for all $V_{\textit{c}}$, and so we may
write that

$\frac{d}{dV_{\textit{c}}}\textrm{Tr}\left[\left(\hat{H}_{\sigma}\phi_{i\sigma}\right)\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}\right]=0.$
(10)

Thus, all terms in Eq. III.2 numerically cancel except for the final term,
$dC/dV_{\textit{c}}$. This may then be re-written using non-degenerate (only
where applicable) first-order perturbation theory, after the present Eq. 6,
since

$\displaystyle\frac{dC}{dV_{\textit{c}}}=\frac{d}{dV_{\textit{c}}}\sum_{\sigma}\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\rho^{\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}$
$\displaystyle={}$ $\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\phi_{i\sigma}^{\ast}\left(\mathbf{r}\right)\frac{d\phi_{i\sigma}\left(\mathbf{r}\right)}{dV_{\textit{c}}}\;d\mathbf{r}+\textit{c.c.}$
$\displaystyle={}$ $\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r^{\prime}}\right)\phi_{i\sigma}^{\ast}\left(\mathbf{r}\right)\sum_{a\neq
i}\phi_{a\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}$
$\displaystyle\times\int\frac{\phi_{a\sigma}^{\ast}\left(\mathbf{r^{\prime}}\right)\frac{dv^{\textrm{KS}}_{\sigma}\left(\mathbf{r^{\prime}}\right)}{dV_{\textit{c}}}\phi_{i\sigma}\left(\mathbf{r^{\prime}}\right)}{\varepsilon_{i\sigma}-\varepsilon_{a\sigma}}\;d\mathbf{r^{\prime}}+\textit{c.c.}$
$\displaystyle={}$ $\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\sum_{a\neq
i}\frac{1}{\varepsilon_{i\sigma}-\varepsilon_{a\sigma}}$
$\displaystyle\times\left(\left(\int\phi_{a\sigma}^{\ast}\left(\mathbf{r}\right)w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\phi_{i\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}\right)\right.$
(11)
$\displaystyle\quad\times\left.\left(\int\phi_{i\sigma}^{\ast}\left(\mathbf{r^{\prime}}\right)\left(\varepsilon^{-1}w_{\textit{c}}\right)^{\sigma}\left(\mathbf{r^{\prime}}\right)\phi_{a\sigma}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r^{\prime}}\right)\right)+\textit{c.c.}.$

The latter expression reduces to W&VV’s Eq. 7 if screening effects are
neglected, that is if we set
$\varepsilon^{-1}_{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)=\delta_{\sigma\sigma^{\prime}}\delta\left(\mathbf{r}-\mathbf{r^{\prime}}\right)$.
Then, as noted by W&VV, the anti-symmetry of the summand implies both that
contributions from $a\leq{N_{\sigma}}$ cancel to zero and may be omitted, and
that the total is strictly non-positive. This condition holds in non-
degenerate, linearly-responding cases of cDFT Lagrange multiplier optimisation
carried out within the potential-update loop of self-consistent field DFT
codes, where W&VV’s result ensures the sign of the energy curvature for any
fixed Kohn-Sham potential.

Globally speaking, however, it appears that the sign of the energy curvature
cannot be inferred directly from the symmetries of Eq. III.2. The necessarily
real-valued screened weight function
$\left(\varepsilon^{-1}w_{\textit{c}}\right)^{\sigma}\left(\mathbf{r}\right)$
may locally vary, even in sign, with respect to
$w_{\textit{c}}\left(\mathbf{r}\right)$, in a complex, system-dependent
manner, typically causing an average net attenuation of the constraining
potential. Even if we may assume that
$0<\left(\varepsilon^{-1}w_{\textit{c}}\right)^{\sigma}\left(\mathbf{r}\right)<w_{\textit{c}}\left(\mathbf{r}\right)$
holds everywhere, for a particular system, and that the orbitals are filled
according to the Aufbau principle in Eq. III.2, the form of this sum offers no
guarantee regarding the sign of $d^{2}W/dV_{\textit{c}}^{2}$.

On the other hand, experience and extensive literature (see review Ref. Kaduk
_et al._ , 2012) yield observations that $d^{2}W/dV_{\textit{c}}^{2}$ is
negative for density constraints applied to a wide variety of systems. As now
we go on to numerically confirm, Eq. 10 provides that this curvature reduces,
for ground-states, to the interacting density response function of the system,
doubly integrated with $w^{\sigma}_{\textit{c}}$. On this basis, we will show
in Section V that any solutions of non-negative curvature are meta-stable
(e.g., the “CS” state of Fig. 1) or unstable (e.g., the inaccessible region in
the vicinity of “CS”) with respect to small perturbations, so that such
curvatures cannot be directly computed and plotted.

The presence of the inverse microscopic dielectric function in Eq. III.2
renders it unsuitable for use in accelerating the convergence of automated
cDFT Lagrange multiplier optimization. Even by using a finite-difference
method for Eq. 7, or indeed in the absence of screening effects, a converged
sum over unoccupied states at each cDFT Lagrange multiplier optimization step
is computationally demanding and conceptually undesirable in DFT. To overcome
these drawbacks, in the following Section we present an alternative approach
for characterizing cDFT energy curvatures, in the framework of response
theory.

##  IV cDFT Energy curvatures from the response function perspective

As an alternative to classifying the cDFT energy landscape by means of
summation over unoccupied states, in what follows we express the relevant
energy curvatures in terms of more computationally convenient integrated
response functions, which depend only on the density or occupied states. This
Section is intended to provide a comprehensive treatment of the relationship
between density response functions and energies in cDFT. We extend our
principal results to the multivariate regime of multiple simultaneous
constraints in Section VI. These response functions are convenient by-products
of automated cDFT optimization, and they are experimental observables in
simulations where the constraining potentials are those of physical fields.
For a cDFT potential representing a uniform electric field, for example,
$w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)$ has a constant gradient, and
the integrated interacting response function represents the high-frequency
electric dipole-dipole polarizability $\mathbf{\alpha}_{\infty}$. Another
example is the microscopic dielectric constant averaged over a region, which
is an important ingredient for supercell convergence acceleration Turban  _et
al._ (2016) and implicit solvation Dziedzic  _et al._ (2011); Andreussi  _et
al._ (2012) schemes. The results provided below allow such response functions
to be calculated, and even updated in self-consistent schemes.

###  IV.1 The integrated interacting response function

Figure 3: (Color online) The numerically evaluated curvature (second
derivatives) of the cDFT total-energy $W$, and its constraint,
$E_{\textit{c}}$, and DFT, $E_{\textrm{DFT}}$, components, with respect to the
Lagrange multiplier, $V_{\textit{c}}$, for the system shown in Fig. 2. Dashed
vertical lines have the same significance as in Fig. 2. The Lagrange
multiplier derivatives of the constrained occupancy $N$, which correspond to
the latter curvatures, are shown, as well as the linear-response approximation
to one of them (crosses). Evaluated values of the unscreened sum-over-states
perturbation theory expression for $d^{2}W/dV_{\textit{c}}^{2}$, following
W&VV, are also shown (solid circles).

The interacting linear response function
$\chi^{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)$
measures the response in the density due to a small change in the applied
potential. Since, by definition,

$\displaystyle d\rho^{\sigma}\left(\mathbf{r}\right)$
$\displaystyle{}=\int\chi^{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)dv^{\textrm{external}}_{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r^{\prime}}$
$\displaystyle{}=\sum_{\sigma^{\prime}}\int\chi^{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)dV_{\textit{c}}^{\sigma^{\prime}}\;d\mathbf{r^{\prime}},$
(12)

one further integration over
$w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)d\mathbf{r}$ allows us to
define the integrated interacting density response function as

$\displaystyle\chi^{\sigma\sigma^{\prime}}\equiv{}$ $\displaystyle\iint
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\chi^{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)d\mathbf{r}\;d\mathbf{r^{\prime}}$
(13) $\displaystyle={}$
$\displaystyle\frac{dN^{\sigma}}{dV^{\sigma^{\prime}}_{\textit{c}}}=\frac{d\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\rho^{\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}}{d\int
w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)V^{\sigma^{\prime}}_{\textit{c}}\;d\mathbf{r^{\prime}}}\int
w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime\prime}}\right)\;d\mathbf{r^{\prime\prime}},$

where we have made it explicit that $dV_{\textit{c}}$ is the average change in
external potential over the subspace. For constraints defined using more
general measures of the density or density matrix, such expressions may be
generalised straightforwardly by replacing the local weighting functions
$w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)$ by non-local projection
operators.

The orthonormality preservation condition expressed in Eq. 10 guarantees the
simplification of Eq. III.2 to Eq. III.2 for all values of $V_{\textit{c}}$,
and so provides that

$\displaystyle\frac{d^{2}W}{dV_{\textit{c}}^{2}}=\frac{dC}{dV_{\textit{c}}}=\frac{dN}{dV_{\textit{c}}}\equiv\chi=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}\chi^{\sigma\sigma^{\prime}},$
(14)

where the $1/2$ is specific to collinear spins. Taking the constraint
contribution alone, on the other hand, we find that the second-order response
function survives, since

$\displaystyle\frac{d^{2}E_{\textit{c}}}{dV_{\textit{c}}^{2}}={}$
$\displaystyle\frac{d^{2}\left(V_{\textit{c}}C\right)}{dV_{\textit{c}}^{2}}=\frac{d}{dV_{\textit{c}}}\left[C+V_{\textit{c}}\frac{dC}{dV_{\textit{c}}}\right]$
$\displaystyle={}$ $\displaystyle
2\frac{dN}{dV_{\textit{c}}}+V_{\textit{c}}\frac{d^{2}N}{dV_{\textit{c}}^{2}}=2\chi+V_{\textit{c}}\frac{d\chi}{dV_{\textit{c}}}.$
(15)

Combining these two results, we deduce a general result for the DFT energy
component curvature, given by

$\displaystyle\frac{d^{2}E_{\textrm{DFT}}}{dV_{\textit{c}}^{2}}={}$
$\displaystyle\frac{d^{2}\left(W-E_{\textit{c}}\right)}{dV_{\textit{c}}^{2}}=-\chi-
V_{\textit{c}}\frac{d\chi}{dV_{\textit{c}}}.$ (16)

In order to check the validity of our approach for analyzing functional
interdependencies and derivatives, and to numerically illustrate our
analytical findings, we performed a cDFT study on the nitrogen molecule shown
in Fig. 2. This serves to illustrate a case in which, if the constraint energy
$E_{\textit{c}}$ vanishes for a finite Lagrange multiplier $V_{\textit{c}}$,
the total-energy $W$ achieves a maximum with respect to $V_{\textit{c}}$. The
total and constraint energies $W$ and $E_{\textit{c}}$ exhibit different
negative curvatures, and their the difference, the DFT component
$E_{\textrm{DFT}}$, necessarily exhibits a positive curvature around the
ground-state minimum at $V_{\textit{c}}=0$ eV. For each value of the charge-
constraining Lagrange multiplier, well converged BLYP Becke (1988); Lee  _et
al._ (1988) ground-state energies and densities, with pseudized $1s$ states,
were calculated using the ONETEP linear-scaling Kohn-Sham DFT code Skylaris
_et al._ (2005). This code solves for the ground-state by optimising a minimal
set of nonorthogonal generalized Wannier functions Skylaris  _et al._ (2002)
in situ. Each of these functions is expanded in an underlying variational
plane-wave equivalent basis set and truncated within a prescribed cut-off
sphere, in this particular case to a radius of $10~{}\mathrm{a}_{0}$. This
approach has been shown to offer finite-difference linear response properties
with an accuracy matching that of conventional plane-wave DFT O’Regan  _et
al._ (2012). The constrained population was defined using the four $2s$ and
$2p$ valence pseudo-orbitals of the isolated atom 444The resulting constraint
acts on the Kohn-Sham density-matrix rather than on the density. Our
analytical findings extends to that case with minor notational changes.. In
the dimer, the resulting unconstrained ground state atomic occupancy was
approximately $6.5$ e, due to overlap between pseudo-orbitals. The charge of
one of the nitrogen atoms was constrained, with the target occupancy set to
$N_{\textit{c}}=6.0$ e.

In Fig. 3, we note a very precise numerical correspondence between the total-
energy curvature, its constraint, and DFT components, and their respectively
predicted reformulations in terms of first and second order integrated
interacting response functions, $\chi=dN/dV_{\textit{c}}$ and
$d\chi/dV_{\textit{c}}$. This serves to confirm that the orthonormality
preservation condition of Eq. 10 holds for all $V_{\textit{c}}$. The failure
of the linear-response approximation $2\chi=2dN/dV_{\textit{c}}$ (shown with
blue crosses) to $d^{2}E_{\textit{c}}/dV_{\textit{c}}$ is somewhat
discouraging for the application of root-finding algorithms on
$E_{\textit{c}}$ in order to optimize the cDFT potential, otherwise a
plausible alternative or compliment to extremizing $W$.

Calculated values for the unscreened sum-over-states perturbation theory
result for $d^{2}W/dV_{\textit{c}}^{2}$ following W&VV, corresponding to the
omission of screening in Eq. III.2, are also shown in Fig. 3 (solid circles).
For this, we generated optimized conduction band states using the method
described in Ref. Ratcliff  _et al._ , 2011, with the conduction band Wannier
function cutoff radii set to $14~{}\mathrm{a}_{0}$. This enabled us to
numerically confirm that unscreened perturbation theory does not generally
match the self-consistent total energy curvature, nor that of its constraint
or DFT energy contributions individually. The anti-symmetry of the unscreened
summand guarantees that it monotonically decreases with an increasing number
of conduction band states, so that the difference between the measured energy
curvature and the unscreened sum-over-states also monotonically increases with
the number of states.

###  IV.2 The integrated non-interacting response function and dielectric
function

In Figs 4 and 5, we illustrate the relationship between energy curvatures and
integrated density response functions with respect to the screened equivalent
of the cDFT Lagrange multiplier, that is the average change in the Kohn-Sham
potential over the constrained region. This provides a test of the magnitude
and potential importance of dielectric screening effects in the energy versus
Lagrange multiplier derivatives of self-consistent cDFT. The relevant weighted
measure of the screened potential, for pure density functionals and
constraints, is given by

$\displaystyle V^{\sigma}=\left(\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)v_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)\;d\mathbf{r}\right)\left(\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r^{\prime}}\right)^{-1}.$
(17)

Then, whereas $\chi=dN/dV_{\textit{c}}$ is the integrated interacting density
response function, we may define $\chi_{0}=dN/dV$ as its non-interacting (also
known as independent-particle) counterpart. If
$\chi^{\sigma\sigma^{\prime}}_{0}\left(\mathbf{r},\mathbf{r^{\prime}}\right)=d\rho^{\sigma}\left(\mathbf{r}\right)/dv^{\textrm{KS}}_{\sigma}\left(\mathbf{r^{\prime}}\right)$,
then

$\displaystyle\chi^{\sigma\sigma^{\prime}}_{0}={}$
$\displaystyle\frac{dN^{\sigma}}{dV^{\sigma^{\prime}}}=\frac{d\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\rho^{\sigma}\left(\mathbf{r}\right)\;d\mathbf{r}}{d\int
w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)v_{\sigma^{\prime}}^{\textrm{KS}}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r^{\prime}}}\int
w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime\prime}}\right)\;d\mathbf{r^{\prime\prime}}$
$\displaystyle\approx{}$ $\displaystyle\iint
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\chi^{\sigma\sigma^{\prime}}_{0}\left(\mathbf{r},\mathbf{r^{\prime}}\right)w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r}\;d\mathbf{r^{\prime}},$
(18)

where the final approximation, a consequence of neglected local-field effects,
becomes an equality in the special case that
$w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)dv^{\textrm{KS}}_{\sigma}\left(\mathbf{r}\right)=w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)dV$
for all $\mathbf{r}$.

Fig. 4 shows that replacing $V_{\textit{c}}$ by $V$ does not preserve any
equivalence between energy curvatures and response functions. Mixed
derivatives with respect to $V_{\textit{c}}$ and $V$, the results of which are
shown in Fig. 5, are required for consistency of screening. Firstly, returning
with the result $dW/dV_{\textit{c}}=C$ for orthonormal states, we find that

$\displaystyle\frac{d^{2}W}{dVdV_{\textit{c}}}=\frac{dC}{dV}=\frac{dN}{dV}\equiv\chi_{0}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}\chi_{0}^{\sigma\sigma^{\prime}},$
(19)

which is numerically confirmed via Fig. 5. Next, we have

$\displaystyle\frac{d^{2}E_{\textit{c}}}{dVdV_{\textit{c}}}={}$
$\displaystyle\frac{d^{2}\left(V_{\textit{c}}C\right)}{dVdV_{\textit{c}}}=\frac{d}{dV}\left[C+V_{\textit{c}}\frac{dC}{dV_{\textit{c}}}\right]$
$\displaystyle={}$
$\displaystyle\frac{dC}{dV}+\frac{dV_{\textit{c}}}{dV}\frac{dC}{dV_{\textit{c}}}+V_{\textit{c}}\frac{d^{2}C}{dVdV_{\textit{c}}}$
$\displaystyle={}$ $\displaystyle
2\frac{dN}{dV}+V_{\textit{c}}\frac{d^{2}N}{dVdV_{\textit{c}}}=2\chi_{0}+V_{\textit{c}}\frac{d\chi_{0}}{dV_{\textit{c}}},$
(20)

since the derivatives with respect to the external and internal potentials
commute. Finally, for the DFT component of the total-energy, we deduce that

$\displaystyle\frac{d^{2}E_{\textrm{DFT}}}{dVdV_{\textit{c}}}={}$
$\displaystyle\frac{d^{2}\left(W-E_{\textit{c}}\right)}{dVdV_{\textit{c}}}=-\frac{dN}{dV}-V_{\textit{c}}\frac{d^{2}N}{dVdV_{\textit{c}}}$
$\displaystyle={}$
$\displaystyle-\chi_{0}-V_{\textit{c}}\frac{d\chi}{dV}=-\chi_{0}-V_{\textit{c}}\frac{d\chi_{0}}{dV_{\textit{c}}}.$
(21)

In practice, these mixed derivatives are calculated simultaneously with the
previously detailed curvatures with respect to $V_{\textit{c}}$, by monitoring
the variation of the weighted, fully relaxed Kohn-Sham potential of Eq. 17 in
self-consistent cDFT calculations. As before, we observe a precise agreement
between the numerically evaluated energy curvatures and response functions.
This confirms Eq. 19, namely, that the averaged non-interacting (i.e.,
independent-particle) response function of DFT may be expressed as an total-
energy landscape property. We return to discuss the unscreened sum-over-states
perturbation theory results shown in Figs. 4 and 5 in Appendix A.

Figure 4: (Color online) Results of the calculations shown in Fig. 3, where
derivatives are instead taken with respect to the change in Kohn-Sham
potential averaged over the constrained region, that is the $\Delta V$ induced
by a finite $V_{\textit{c}}$, using the same weight function as used to
calculate the constrained occupancy $N$. The derivative $dN/dV$ corresponds to
the average non-interacting charge response of the constrained region, whereas
$dN/dV_{\textit{c}}$ is the interacting response. Unlike the derivatives shown
in Fig. 3, these non-interacting (i.e., bare or independent-particle)
derivatives demonstrate no equivalence. The perturbation theory data points
are as in Fig. 3.

The process of Lagrange multiplier optimization offers ready access to the
physical response properties of the constrained region, which are not limited
to those of the target state. The simplest such quantity is the subspace
inverse dielectric constant (i.e., screening factor),

$\displaystyle\epsilon^{-1}_{\sigma\sigma^{\prime}}\equiv{}$
$\displaystyle\iint
w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\epsilon^{-1}_{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r}\;d\mathbf{r^{\prime}}$
$\displaystyle\times\Big{(}\int
w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime\prime}}\right)\;d\mathbf{r^{\prime\prime}}\Big{)}^{-1}.$
(22)

This may be calculated directly using the cDFT integrated response functions,
since, by definition,

$\displaystyle dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)={}$
$\displaystyle\int\epsilon^{-1}_{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)dv_{\sigma^{\prime}}^{\textrm{external}}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r^{\prime}}$
$\displaystyle={}$
$\displaystyle\int\epsilon^{-1}_{\sigma\sigma^{\prime}}\left(\mathbf{r},\mathbf{r^{\prime}}\right)w_{\textit{c}}^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)dV_{\textit{c}}^{\sigma^{\prime}}\;d\mathbf{r^{\prime}}$
$\displaystyle\Rightarrow\quad dV^{\sigma}={}$
$\displaystyle\epsilon^{-1}_{\sigma\sigma^{\prime}}dV_{\textit{c}}^{\sigma^{\prime}}.$
(23)

From this, it is clear that $\epsilon^{-1}_{\sigma\sigma^{\prime}}$ is a
property of the constrained ground-state density and is not explicitly
dependent on unoccupied states. Next, we may apply the chain rule via the the
constrained property $N^{\sigma}$, introducing the notation
$\chi_{0}^{-1\sigma\sigma^{\prime\prime}}=dV^{\sigma}/dN^{\sigma^{\prime\prime}}$,
which provides that

$\displaystyle\epsilon^{-1}_{\sigma\sigma^{\prime}}={}$
$\displaystyle\frac{dV^{\sigma}}{dV_{\textit{c}}^{\sigma^{\prime}}}=\sum_{\sigma^{\prime\prime}}\chi_{0}^{-1\sigma\sigma^{\prime\prime}}\chi^{\sigma^{\prime\prime}\sigma^{\prime}}.$
(24)

The inverse of this quantity is the subspace-averaged dielectric constant,
neglecting local-field effects, given by

$\displaystyle\epsilon_{\sigma\sigma^{\prime}}={}$
$\displaystyle\frac{dV_{\textit{c}}^{\sigma^{\prime}}}{dV^{\sigma}}=\sum_{\sigma^{\prime\prime}}\chi_{0}^{\sigma\sigma^{\prime\prime}}\chi^{-1\sigma^{\prime\prime}\sigma^{\prime}}.$
(25)

cDFT thus provides an efficient means of estimating the dielectric constants
of spatial regions, which are central in implicit solvation Dziedzic  _et al._
(2011); Andreussi  _et al._ (2012), supercell convergence acceleration of
excitation energies Turban  _et al._ (2016), and in high-frequency optical
response, in terms of by-products of its optimization.

Figure 5: (Color online) Following from Fig. 4, the correspondence between
energy curvatures and averaged non-interacting response functions is recovered
by using mixed interacting and non-interacting second derivatives, which
restores their consistency. The sum-over-states perturbation theory (P.T.)
data points following W&VV are again shown (solid circles), together with the
“Renormalized P.T.” data points (open black circles) described in Appendix A.

##  V Stable cDFT solutions are energy maxima with respect to their Lagrange
multiplier

The response function based approach to cDFT analysis is used in this Section
to globally characterize its stationary points. W&VV have shown that a non-
negative total-energy curvature with respect to the cDFT Lagrange multiplier
is guaranteed in cases where non-degenerate perturbation theory is applicable
Wu and Van Voorhis (2005), in the absence of screening effects. This regime
may hold during cDFT optimization within the density-update loop of self-
consistent field codes, where the Kohn-Sham potential is fixed Wu and Yang
(2003); Wu and Van Voorhis (2006a). More generally, or physically, the Kohn-
Sham potential is relaxed self-consistently for each $V_{\textit{c}}$, and
sum-over-states perturbation theory becomes inconclusive as discussed in the
text surrounding Eq. III.2. We provide a general, non-perturbative proof that
stable cDFT solutions always occur at energy maxima with respect to their
Lagrange multiplier. We then extend this proof to the multivariate cDFT regime
in Section VI.

We have shown that the curvatures of the total-energy with respect to the cDFT
Lagrange multiplier $V_{\textit{c}}$ are equal to integrated density response
functions. Next, we will show that, similarly, the curvatures with respect to
the cDFT occupancy $N$, for a particular target $N_{\textit{c}}$, are equal to
the inverses of these response functions. The occupancy $N$ is not a free
parameter, and so, in order to analyze these curvatures, it is convenient to
initially take derivatives with respect to the cDFT target occupancy
$N_{\textit{c}}$, subject to the condition that $N=N_{\textit{c}}$ for a
suitable Lagrange multiplier $V_{\textit{c}}\left(N_{\textit{c}}\right)$.
Echoing Eq. III.1, assuming that the ground-state is located for all
$V_{\textit{c}}$,

$\displaystyle\frac{dW}{dN_{\textit{c}}}={}$
$\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\left(\frac{\delta
W}{\delta\phi_{i\sigma}^{\ast}}\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}+\textit{c.c.}\right)\frac{dV_{\textit{c}}}{dN_{\textit{c}}}\right]+\frac{\partial
W}{\partial N_{\textit{c}}}$ $\displaystyle={}$
$\displaystyle\sum_{\sigma}\sum_{i}^{N_{\sigma}}\textrm{Tr}\left[\left(\left(\hat{H}_{\sigma}\phi_{i\sigma}\right)\frac{d\phi_{i\sigma}^{\ast}}{dV_{\textit{c}}}+\textit{c.c.}\right)\frac{dV_{\textit{c}}}{dN_{\textit{c}}}\right]-V_{\textit{c}}$
$\displaystyle={}$ $\displaystyle-
V_{\textit{c}}\left(N_{\textit{c}}\right)\quad\Rightarrow\quad\frac{d^{2}W\left(N_{\textit{c}}\right)}{dN_{\textit{c}}^{2}}=-\frac{dV_{\textit{c}}}{dN_{\textit{c}}}.$
(26)

Since the constraint $N=N_{\textit{c}}$ is satisfied for each
$N_{\textit{c}}$, the constraint energy always vanishes and
$E_{\textrm{DFT}}\left(N\right)=W\left(V_{\textit{c}}\left(N_{\textit{c}}\right)\right)$
along the curve. Thus, we may write, for the occupancy curvature of the DFT
contribution, that

$\displaystyle\frac{d^{2}E_{\textrm{DFT}}\left(N\right)}{dN^{2}}=\frac{d^{2}W\left(V_{\textit{c}}\left(N_{\textit{c}}\right)\right)}{dN_{\textit{c}}^{2}}=-\frac{dV_{\textit{c}}}{dN}=-\chi^{-1}.$
(27)

As previously established by W&VV for the unscreened case Wu and Van Voorhis
(2005), the combination of Eqs. 14 and 27 provides that the curvature of the
total-energy with respect to Lagrange multiplier is directly related to the
curvature of the DFT energy with respect to the occupancy, namely

$\displaystyle\left(\frac{d^{2}W}{dV_{\textit{c}}^{2}}\right)^{-1}$
$\displaystyle{}=\left(\frac{dN}{dV_{\textit{c}}}\right)^{-1}=\chi^{-1}$
$\displaystyle{}=\frac{dV_{\textit{c}}}{dN}=-\frac{d^{2}E_{\textrm{DFT}}}{dN^{2}}.$
(28)

We may extend Eq. 28 to the general cDFT total-energy $W$, no longer subject
the constraint that the target occupancy is attained, by freeing the target
$N_{\textit{c}}$ after optimization of $V_{\textit{c}}$. This has no bearing
on $E_{\textrm{DFT}}$, by definition, and so Eq. 27 holds irrespective of
whether the constraint is satisfied. For $W\left(V_{\textit{c}}\right)$ with
$N$ not necessarily equal to $N_{\textit{c}}$, it is sufficient to add the
curvature of the now non-vanishing constraint energy term $E_{\textit{c}}$,
given by

$\displaystyle\frac{d^{2}E_{\textit{c}}}{dN^{2}}={}$
$\displaystyle\frac{d^{2}\left(V_{\textit{c}}C\right)}{dN^{2}}=\frac{d}{dN}\left[V_{\textit{c}}\frac{dC}{dN}+C\frac{dV_{\textit{c}}}{dN}\right]$
(29) $\displaystyle={}$ $\displaystyle
2\frac{dV_{\textit{c}}}{dN}\frac{dC}{dN}+C\frac{d^{2}V_{\textit{c}}}{dN^{2}}=2\frac{dV_{\textit{c}}}{dN}+\left(N-N_{\textit{c}}\right)\frac{d^{2}V_{\textit{c}}}{dN^{2}}.$

The total-energy curvature is provided by the sum

$\displaystyle\frac{d^{2}W}{dN^{2}}={}$
$\displaystyle\frac{d^{2}E_{\textrm{DFT}}}{dN^{2}}+\frac{d^{2}E_{\textit{c}}}{dN^{2}}$
$\displaystyle={}$
$\displaystyle\frac{dV_{\textit{c}}}{dN}+\left(N-N_{\textit{c}}\right)\frac{d^{2}V_{\textit{c}}}{dN^{2}}.$
(30)

Then, by combining Eqs. 27 to 30 and by applying the constraint condition
$N=N_{\textit{c}}$, we arrive at our central result that, valid for all cDFT
stationary points as defined,

$\displaystyle\left(\frac{d^{2}W}{dV_{\textit{c}}^{2}}\right)^{-1}$
$\displaystyle{}=\frac{1}{2}\frac{d^{2}E_{\textit{c}}}{dN^{2}}=-\frac{d^{2}E_{\textrm{DFT}}}{dN^{2}}=\frac{d^{2}W}{dN^{2}}.$
(31)

This is numerically confirmed via Figs. 3 and 6, and the inaccuracy of linear-
response approximations in Fig. 6 (shown with crosses), except at
$N=N_{\textit{c}}$, is clear.

Figure 6: (Color online) Numerically evaluated second derivatives of the cDFT
total-energy $W$, its constraint $E_{\textit{c}}$ and DFT $E_{\textrm{DFT}}$
components, with respect to the constrained occupancy $N$, for the system
shown in Fig. 2 and occupancy target $N_{\textit{c}}=6~{}\textrm{e}$. Dashed
vertical lines show the occupancy at which $W$ (green) and $E_{\textit{c}}$
(red) are maximized. The corresponding occupancy derivatives of the Lagrange
multiplier (averaged inverse screened response functions) are shown, as well
as their non-corresponding linear-response approximations (crosses).

Next, let us consider the possibility of simulating a stable constrained
ground-state at a given value of $V_{\textit{c}}$, which may be a cDFT
solution but is not necessarily so. Stability implies that the energy is
minimized with respect to the density, locally at least, and put more
precisely, that the energy is locally strictly convex. Mathematically, this
condition is represented by a positive definite matrix
$d^{2}W/d\rho^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)d\rho^{\sigma}\left(\mathbf{r}\right)$,
evaluated at a fixed $V_{\textit{c}}$, where the coordinate pairs
$\left\\{\sigma,\mathbf{r}\right\\}$ form the basis vectors. Separating the
DFT and constraint energy terms, the stability of constrained ground-states is
then determined by

$\displaystyle\left.\frac{d^{2}W}{d\rho^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)d\rho^{\sigma}\left(\mathbf{r}\right)}\right|_{V_{\textit{c}},\hat{w}_{\textit{c}}}$
$\displaystyle{}=\frac{d^{2}E_{\textrm{DFT}}}{d\rho^{\sigma^{\prime}}\left(\mathbf{r^{\prime}}\right)d\rho^{\sigma}\left(\mathbf{r}\right)},$
(32)

so that the latter is also positive definite. The underlying DFT energy
landscape therefore alone determines the stability of constrained ground-
states and thus, while cDFT may relocate the minima of $W$ within the simply
connected domains of stability of the unconstrained DFT problem, it cannot
deform those domains. This is a consequence of the linearity of
$E_{\textit{c}}$ in the density and thus its vanishing contribution to the
fixed-$V_{\textit{c}}$ curvature of $W$. The constraint may, however cause
abrupt transitions of energy minima across disconnected domains, as observed
in the “OS(-)” to “OS(+)” transition of Fig. 1 (albeit in the more general
case of a spin-dependent constraint).

In order to analyze the sign of the quantity expressed in Eq. 31, we may focus
on the DFT energy curvature with respect to the subspace occupancy,
$d^{2}E_{\textrm{DFT}}/dN^{2}$. Let us suppose that stable ground-states are
observed, in a particular constrained system, as $V_{\textit{c}}$ is
continuously varied within an interval. The local special case of Eq. 26
provides that
$dE_{\textrm{DFT}}/d\rho^{\sigma}\left(\mathbf{r}\right)=-v^{\textrm{external}}_{\sigma}\left(\mathbf{r}\right)=-w_{\textit{c}}\left(\mathbf{r}\right)V_{\textit{c}}$.
Combining this with the general properties of a positive definite matrix,
i.e., that its inverse is positive definite and that its diagonal elements are
positive, we find that

$\displaystyle 0$
$\displaystyle{}<\left(\frac{d^{2}E_{\textrm{DFT}}}{d\rho^{\sigma
2}}\right)^{-1}\left(\mathbf{r},\mathbf{r}\right)=\left(-\frac{dv_{\sigma}^{\textrm{external}}}{d\rho^{\sigma}}\right)^{-1}\left(\mathbf{r},\mathbf{r}\right)$
$\displaystyle{}=-\chi^{\sigma\sigma}\left(\mathbf{r},\mathbf{r}\right)=-\frac{d\rho^{\sigma}\left(\mathbf{r}\right)}{dv_{\sigma}^{\textrm{external}}\left(\mathbf{r}\right)}=-\frac{d\rho^{\sigma}\left(\mathbf{r}\right)}{w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)dV_{\textit{c}}}.$
(33)

Following from this, we may multiply by
$\left(w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\right)^{2}$, allowing
$w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)$ to take zero (but not at all
$\mathbf{r}$) and negative values, as is necessary, for example, when
constraining the difference in charge between two regions. Then, integrating
over all space, we arrive at the result

$\displaystyle 0$
$\displaystyle{}>\sum_{\sigma}\int\left(w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\right)^{2}\frac{d\rho^{\sigma}\left(\mathbf{r}\right)}{w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)dV_{\textit{c}}}\;d\mathbf{r}=\frac{dN}{dV_{\textit{c}}}=\chi.$
(34)

Thus, Eq. 33, which states that the local part of the interacting response
function is negative about stable ground-states, is generalized in Eq. 34 to
the cDFT response function for an arbitrary weight function, irrespective of
its profile or even its sign. Recalling Eq. 14, we arrive at the sought-after
result that $d^{2}W/dV_{\textit{c}}^{2}<0$ over stable ground-states. We note
that this result is valid up to arbitrary order in response, it is general to
degenerate and non-degenerate systems, and it accounts for screening.

Conversely, let us briefly suppose that $d^{2}W/dV_{\textit{c}}^{2}\geq 0$ for
some value of $V_{\textit{c}}$. Then, by virtue of Eq. 14, we have
$dN/dV_{\textit{c}}\geq 0$ and, combining Eq. 33, Eq. 34, and the non-
negativity of
$\left(w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\right)^{2}$, there exist
$\mathbf{r}$ and $\sigma$ for which
$w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\neq 0$ and
$0\geq\left(d^{2}E_{\textrm{DFT}}/d\rho^{\sigma
2}\right)^{-1}\left(\mathbf{r},\mathbf{r}\right)$. Then,
$\big{(}d^{2}E_{\textrm{DFT}}/d\rho^{\sigma}d\rho^{\sigma^{\prime}}\big{)}^{-1}\left(\mathbf{r},\mathbf{r^{\prime}}\right)$
may not be positive definite since it has a non-positive diagonal element.
Therefore, its inverse
$\big{(}d^{2}E_{\textrm{DFT}}/d\rho^{\sigma}d\rho^{\sigma^{\prime}}\big{)}\left(\mathbf{r},\mathbf{r^{\prime}}\right)$
may also not be positive definite, and this guarantees that the state is
unstable with respect to spin-density perturbations. Then, only meta-stable
states such as the “CS” state depicted in Fig. 1 may be observed in practice.
There are two ways in which cDFT inflection points may arise, that is, cases
where $d^{2}W/dV_{\textit{c}}^{2}=dN/dV_{\textit{c}}=0$. First, the constraint
may couple pathologically to the density or not couple to it at all. Trivial
examples of this are constraints where, for all $\sigma$ and $\mathbf{r}$,
$w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)=1$ or
$w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\rho^{\sigma}\left(\mathbf{r}\right)=0$,
which result in a vanishing response $\chi$. Second, the state may be
degenerate on a line or higher surface, and the constraint may be contrived so
as not to break this degeneracy.

Combining these arguments, we conclude that the stability of ground-states
over an interval of $V_{\textit{c}}$ both implies and requires that
$d^{2}W/dV_{\textit{c}}^{2}$ remains negative and finite at constrained
ground-states within that interval. In practice, this means that intervals of
positive curvature are not numerically observable, because unstable or meta-
stable states cannot be sampled in a continuous manner. It is however
possible, as demonstrated in Fig. 1, to observe numerical discontinuities in
the otherwise concave total-energy $W$, for fixed $N_{\textit{c}}$, when the
response, and hence energy curvature, diverges at a phase transition. Such
transitions occur at vanishing values of the quantity expressed in Eq. 31,
i.e., at inflection points of $E_{\textrm{DFT}}\left(N\right)$. For this, it
is necessary for $\chi^{\sigma\sigma}\left(\mathbf{r},\mathbf{r}\right)$ to
diverge for some $\sigma$ and $\mathbf{r}$ where
$w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\neq 0$, by virtue of Eqs. 28,
33 and 34.

##  VI Multiple constraints, the cDFT condition number, and the extension to
self-consistent field cDFT

In this section, we consider, in detail, the extent to which our principal
results may be generalized to multiple constraints, and thereafter to self-
consistent field cDFT. For the former, it is helpful to consider a vector
Lagrange multiplier $\mathbf{V}_{\textit{c}}$, acting upon a vector of
constraint functionals $\mathbf{C}=\mathbf{N}-\mathbf{N}_{\textit{c}}$,
yielding $E_{\textit{c}}=\mathbf{V}_{\textit{c}}\cdot\mathbf{C}$. The
multivariate generalization of Eq. 28 is the equation

$\displaystyle\left(\frac{d^{2}W}{d\mathbf{V}_{\textit{c}}^{2}}\right)^{-1}$
$\displaystyle{}=\left(\frac{d\mathbf{N}}{d\mathbf{V}_{\textit{c}}}\right)^{-1}=\boldsymbol{\chi}^{-1}$
$\displaystyle{}=\frac{d\mathbf{V}_{\textit{c}}}{d\mathbf{N}}=-\frac{d^{2}E_{\textrm{DFT}}}{d\mathbf{N}^{2}},$
(35)

where the negative exponent denotes matrix inversion, and the matrix is
symmetric due to the symmetry of
$\chi\left(\mathbf{r},\mathbf{r^{\prime}}\right)$ under co-ordinate exchange
in Eq. 13. We may also consider the non-interacting analogue of Eq. 35, viz.

$\displaystyle\left(\frac{d^{2}W}{d\mathbf{V}d\mathbf{V}_{\textit{c}}}\right)^{-1}$
$\displaystyle{}=\left(\frac{d\mathbf{N}}{d\mathbf{V}}\right)^{-1}=\boldsymbol{\chi}^{-1}_{0}.$
(36)

Similar expressions may be derived for the second derivatives of
$E_{\textrm{DFT}}$ and $E_{\textit{c}}$. For example, Eq. 21 becomes, noting
that the final object is a rank-three tensor,

$\displaystyle\frac{d^{2}E_{\textrm{DFT}}}{d\mathbf{V}d\mathbf{V}_{\textit{c}}}=-\boldsymbol{\chi}_{0}-\mathbf{V}_{\textit{c}}\cdot\frac{d\boldsymbol{\chi}_{0}}{d\mathbf{V}_{\textit{c}}}.$
(37)

The dielectric function may also be straightforwardly generalised in order to
describe the coupling of dielectric screening between constrained subspaces.
In the spin-degenerate or spin-averaged case, it takes the particularly simple
matrix form

$\displaystyle\boldsymbol{\epsilon}={}$
$\displaystyle\frac{d\textbf{V}_{\textit{c}}}{d\textbf{V}}=\frac{1}{2}\boldsymbol{\chi}_{0}\cdot\boldsymbol{\chi}^{-1}.$
(38)

The behaviour of multivariate cDFT optimisation depends on the nature of the
energy landscape with respect to $\textbf{V}_{\textit{c}}$, and it may be
complicated by the presence of multiple extrema or by surfaces which are
elongatated in some directions relatively to others. The curvature of the
energy landscape, about a particular value of $\textbf{V}_{\textit{c}}$, is
characterized by the eigenvalues of the Hessian matrix $\boldsymbol{\chi}$.
These must be identical for optimal extremization of
$W\left[\mathbf{V}_{\textit{c}}\right]$ by simple conjugate gradients. More
generally, the positive-valued condition number $k$, defined by

$\displaystyle
k^{2}=\sum_{IJ}\left(\left[\frac{d^{2}W}{d\mathbf{V}_{\textit{c}}^{2}}\right]_{IJ}\right)^{2}\sum_{KL}\left(\left[\frac{d^{2}W}{d\mathbf{V}_{\textit{c}}^{2}}\right]_{KL}^{-1}\right)^{2},$
(39)

where capitalized letters index constraints, gives a measure of the spread of
the Hessian eigenvalues. The optimization problem is said to be ill-
conditioned if $k$ is much larger than one. Then, conjugate gradients may be
expected to perform poorly, and Newton or quasi-Newton methods may be
considered. We may rephrase $k$ as

$\displaystyle k^{2}$
$\displaystyle{}=\sum_{IJ}\left(\left[\frac{d\mathbf{N}}{d\mathbf{V}_{\textit{c}}}\right]_{IJ}\right)^{2}\sum_{KL}\left(\left[\frac{d\mathbf{V}_{\textit{c}}}{d\mathbf{N}}\right]_{KL}\right)^{2}$
(40)
$\displaystyle{}=\sum_{IJ}\left(\left[\boldsymbol{\chi}\right]_{IJ}\right)^{2}\sum_{KL}\left(\left[\boldsymbol{\chi}^{-1}\right]_{KL}\right)^{2}$
$\displaystyle{}=\sum_{IJ}\left(\left[\left(\frac{d^{2}E_{\textrm{DFT}}}{d\mathbf{N}^{2}}\right)^{-1}\right]_{IJ}\right)^{2}\sum_{KL}\left(\left[\frac{d^{2}E_{\textrm{DFT}}}{d\mathbf{N}^{2}}\right]_{KL}\right)^{2},$

through which it is clear that $W\left[\mathbf{V}_{\textit{c}}\right]$ has the
same condition number as
$E_{\textrm{DFT}}\left[\mathbf{N}_{\textit{c}}\right]$. If optimising
$\mathbf{V}_{\textit{c}}$ during the density update step of self-consistent
field DFT codes Wu and Yang (2003); Wu and Van Voorhis (2006a), the
corresponding condition number $k_{0}$ may be calculated using the fixed-
potential equivalent of Eq. 40. Both $k$ and $k_{0}$ are amenable to
evaluation by finite differences. If the response $\boldsymbol{\chi}$ or
inverse response $\boldsymbol{\chi}^{-1}$ matrices are singular, $k$ diverges,
the problem is said to be ill-posed, and its definition must be reconsidered.
We refer the reader to Ref. Rondinelli  _et al._ , 2007 for a more general
discussion of condition numbers in the context of electronic structure.

The connection between the stability of constrained states and the negativity
of $d^{2}W/dV_{\textit{c}}^{2}$, as demonstrated in Section V, is here
generalized to the multivariate case. More precisely, the stability of
multiply-constrained ground-states both implies and requires that the matrix
$d^{2}W/d\mathbf{V}_{\textit{c}}^{2}$ has strictly negative diagonal elements.
If the matrix is strictly diagonally dominant, furthermore, stability implies
and requires that $W\left(\mathbf{V}_{\textit{c}}\right)$ is strictly concave
locally. To see this, let us assume that stable ground-states are observed as
$\mathbf{V}_{\textit{c}}$ is varied within some simply connected volume. Then,
taking the diagonal matrix elements of $d^{2}W/d\mathbf{V}_{\textit{c}}^{2}$
individually, while following Eqs. 33 and 34 for each constraint labelled $I$,
we find that

$\displaystyle 0$
$\displaystyle{}>\sum_{\sigma}\int\left(w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)\right)^{2}\frac{d\rho^{\sigma}\left(\mathbf{r}\right)}{w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)dV_{\textit{c}I}}\;d\mathbf{r}=\left[\frac{d\mathbf{N}}{d\mathbf{V}_{\textit{c}}}\right]_{II}.$
(41)

Combining this with Eq. 35, we deduce that the diagonal matrix elements of
$d^{2}W/d\mathbf{V}_{\textit{c}}^{2}$ are all negative. This is a necessary
but not sufficient condition for $d^{2}W/d\mathbf{V}_{\textit{c}}^{2}$ to be
negative definite locally. However, it becomes a sufficient condition for
negative definiteness if the symmetric matrix is also diagonally dominant
Harville (2011). Diagonal dominance is a physically reasonable, albeit system-
dependent, condition related to the extent of non-locality in
$\chi\left(\mathbf{r},\mathbf{r^{\prime}}\right)$. If it holds, stability
implies that $W\left(\mathbf{V}_{\textit{c}}\right)$ is strictly concave
locally. Additionally, if
$w^{\sigma}_{\textit{c}J}\left(\mathbf{r}\right)w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)>0$
for all $I$, $J$, $\sigma$, and $\mathbf{r}$ (such as if
$w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)>0$ for all $I$ and $\sigma$),
then

$\displaystyle 0$ $\displaystyle{}>\sum_{\sigma}\int
w^{\sigma}_{\textit{c}J}\left(\mathbf{r}\right)w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)\frac{d\rho^{\sigma}\left(\mathbf{r}\right)}{w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)dV_{\textit{c}I}}\;d\mathbf{r},$
(42)

and the Hessian is a strictly negative matrix, but that alone does not imply
its negative definiteness. Conversely, if
$d^{2}W/d\mathbf{V}_{\textit{c}}^{2}$ does not exhibit all negative diagonal
elements, as is possible if it is not negative definite, for example, then,
$\left[d\mathbf{N}/d\mathbf{V}_{\textit{c}}\right]_{II}\geq 0$ for one or more
values of $I$. Then, by the non-negativity of
$\left(w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)\right)^{2}$ for all
$I$, $\sigma$, and $\mathbf{r}$, the state is unstable, as previously shown in
the single-constraint case. Thus, we conclude that if a ground-state subject
to multiple linear, non-trivial constraints is stable with respect to
perturbations, and hence locatable by numerical cDFT optimization, this
implies and requires that the associated Hessian matrix
$d^{2}W/d\mathbf{V}_{\textit{c}}$ has a strictly negative diagonal, and that
it is negative definite in the event that it is diagonally dominant. If it is
not strictly diagonally dominant, however, we cannot rule out the possibility
of stable ground-states occurring at saddle points of the self-consistent
$W\left(\mathbf{V}_{\textit{c}}\right)$.

Finally, we consider the energy landscape associated with cDFT Lagrange
multiplier optimization within the fixed-potential inner loop of self-
consistent field DFT codes. For this, let us first assume that stable minimum
energy eigenstates are observed, for a given fixed DFT contribution
$\hat{v}^{\textrm{DFT}}_{\sigma}$ to the Kohn-Sham potential, as
$\mathbf{V}_{\textit{c}}$ is varied within some simply connected volume.
Again, we may analyze the diagonal matrix elements of
$d^{2}W/d\mathbf{V}_{\textit{c}}^{2}$ individually, but in this case
$dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)=dv_{\sigma}^{\textrm{external}}\left(\mathbf{r}\right)=w_{\textit{c}I}^{\sigma}\left(\mathbf{r}\right)dV_{\textit{c}I}^{\sigma}$,
for each $I$, and for all $\sigma$ and $\mathbf{r}$. Then, Eq. 33 may be
modified for a fixed DFT potential, to give

$\displaystyle 0$
$\displaystyle{}<\left(\frac{d^{2}E_{\textrm{DFT}}}{d\rho^{\sigma
2}}\right)^{-1}_{\hat{v}^{\textrm{DFT}}_{\sigma}}\left(\mathbf{r},\mathbf{r}\right)=-\left.\frac{d\rho^{\sigma}\left(\mathbf{r}\right)}{w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)dV_{\textit{c}I}}\right|_{\hat{v}^{\textrm{DFT}}_{\sigma}},$
(43)

whereupon Eq. 41 may be similarly adapted, yielding

$\displaystyle 0$
$\displaystyle{}>\sum_{\sigma}\int\left.\frac{w^{\sigma}_{\textit{c}I}\left(\mathbf{r}\right)d\rho^{\sigma}\left(\mathbf{r}\right)}{dV_{\textit{c}I}}\right|_{\hat{v}^{\textrm{DFT}}_{\sigma}}\;d\mathbf{r}=\left[\frac{d\mathbf{N}}{d\mathbf{V}_{\textit{c}}}\right]_{II}^{\hat{v}^{\textrm{DFT}}_{\sigma}}.$
(44)

The identity
$d^{2}\mathbf{W}/dV_{\textit{c}}^{2}=d\mathbf{N}/d\mathbf{V}_{\textit{c}}$
remains valid for a fixed $\hat{v}^{\textrm{DFT}}_{\sigma}$, since it relies
only on the orthonormality of Kohn-Sham eigenstates. Thus, the diagonal
elements of $d^{2}\mathbf{W}/dV_{\textit{c}}^{2}$ are individually negative,
and the guarantee of strict local concavity of
$W\left(\mathbf{V}_{\textit{c}}\right)$, for any given
$\hat{v}_{\sigma}^{\textrm{DFT}}$, is again conditional on a symmetric
diagonally dominant Harville (2011) $d^{2}\mathbf{W}/dV_{\textit{c}}^{2}$. If
the matrix does not exhibit a strictly negative diagonal, conversely, then
$\left[d\mathbf{N}/d\mathbf{V}_{\textit{c}}\right]_{II}\geq 0$ for one or more
values of $I$, and there exist $\mathbf{r}$ and $\sigma$ for which
$w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\neq 0$ and
$0\geq\left(d^{2}E_{\textrm{DFT}}/d\rho^{\sigma
2}\right)^{-1}\left(\mathbf{r},\mathbf{r}\right)$ at fixed
$\hat{v}_{\sigma}^{\textrm{DFT}}$. Then, the non-self-consistent
$\big{(}d^{2}E_{\textrm{DFT}}/d\rho^{\sigma}d\rho^{\sigma^{\prime}}\big{)}\left(\mathbf{r},\mathbf{r^{\prime}}\right)$
cannot be positive definite, and the state is unstable. We conclude that for a
minimum-energy eigenstate, given a fixed DFT potential and subject to one or
more linear, non-trivial constraints, to be stable with respect to
perturbations and hence locatable by SCF-type cDFT optimization, it implies
and requires the associated fixed-$\hat{v}_{\sigma}^{\textrm{DFT}}$ Hessian
matrix $d^{2}W/d\mathbf{V}_{\textit{c}}$ to have a strictly negative diagonal,
and for it to be negative definite in the event that it is diagonally
dominant. This generalizes the results of Refs. Wu and Van Voorhis, 2005,
2006a; Kaduk  _et al._ , 2012 beyond the first-order non-degenerate
perturbative regime, extending their validity to systems which exhibit non-
linear response or orbital degeneracy.

##  VII Conclusion and Summary

Constrained DFT is an flexible, potent approach that broadens the scope and
flexibility of DFT-based atomistic simulation. A growing number of software
implementations of cDFT are now appearing Kaduk  _et al._ (2012); Valiev  _et
al._ (2010); Shao  _et al._ (2015); Oberhofer and Blumberger (2009b, 2010a);
Souza  _et al._ (2013); Brooke  _et al._ (2015); R̆ezác̆  _et al._ (2012);
Kurz  _et al._ (2004); Ozaki and Kino (2005); Hutter  _et al._ (2014); Ma and
Dudarev (2015), including linear-scaling implementations designed for
application to large systems Sena  _et al._ (2011); Ratcliff  _et al._ (2015).
It is inherently parallelizable, and thus potentially suitable for use in
combination with high-throughput materials screening infrastructures Jain  _et
al._ (2011); Curtarolo  _et al._ (2013). Fundamental developments will be
required to bring cDFT into the realm of such very routine use. For
transferability and comparability between codes, for example, the
standardized, automated selection of population analysis and targeting
schemes, ideally but not necessarily based on energy considerations, would
surely be beneficial. Methods based on promolecule densities Wu  _et al._
(2007) or self-consistent Wannier functions O’Regan  _et al._ (2010) are
promising possibilities in this direction. The key findings and conclusions of
the present work are that:

  * •

For any self-consistent energy second-derivative with respect to a constrained
expectation value, there is an equivalent integrated linear response or
inverse-response function which is more convenient to calculate. This provides
a basis for future quasi-Newton or preconditioning approaches for cDFT.

  * •

A negative-diagonal self-consistent cDFT Hessian is implied by the local
stability of the system with respect to perturbations about the point of
evaluation. A Hessian lacking this property cannot be observed by finite
differences, and so may be excluded in automation. This generalizes W&VV’s
result to the self-consistent cDFT problem, the non-linear response regime,
and to degenerate orbitals.

  * •

Concave regions in the DFT energy versus spin-density landscape cannot be
explored using cDFT.

  * •

Integrated response and dielectric functions may be evaluated as by-products
of cDFT optimization, without sums over empty Kohn-Sham eigenstates.

  * •

Existing cDFT optimization schemes based on the self-consistent field
approach, which update the cDFT Lagrange multipliers with fixed DFT
potentials, may be readily adapted for the self-consistent gradients required
in direct-minimization DFT.

  * •

cDFT does not change the domains of stability of the underlying DFT energy
landscape, it moves the solutions around, or between, these domains. When
there are two or more such domains in cDFT, it is possible to observe multiple
solutions, hysteresis, and energy discontinuities at such transitions.

We expect that this work may facilitate the advanced automation of cDFT
Lagrange multiplier optimization, particularly in the high-throughput,
molecular dynamics, and linear-scaling regimes. Our general framework for
treating energy landscapes in terms of integrated response functions now
enables the extension of cDFT to new areas of atomistic and continuum
simulation.

## Acknowledgements

This work was enabled by the Royal Irish Academy – Royal Society International
Exchange Cost Share Programme (IE131505). GT acknowledges support from EPSRC
UK (EP/I004483/1 and EP/K013610/1). We acknowledge and thank the ONETEP
Developer’s Group for hosting workshops that facilitated this development, and
Qin Wu for a helpful communication and discussion.

##  Appendix A Approximate connection of the cDFT non-interacting response
function to unscreened sum-over-states perturbation theory

In this Appendix, we investigate how the non-interacting linear-response
function defined by $\chi_{0}=dN/dV=d^{2}W/dVdV_{\textit{c}}$ in Eq. 19
relates to the unscreened sum-over-states perturbation theory expression given
by Eq. 7 of Ref. Wu and Van Voorhis, 2005. The latter is obtained by replacing
$\left(\varepsilon^{-1}w_{\textit{c}}\right)^{\sigma}$ by
$w_{\textit{c}}^{\sigma}$ in Eq. III.2. It yields the linear response of the
targeted density $N$ to a change in the subspace-averaged Kohn-Sham potential
$V$ with all other aspects of the Kohn-Sham potential kept fixed. Its
evaluation for two values of $V_{\textit{c}}$, close to those at which
$E_{\textit{c}}$ and $W$ assume maxima, is shown in Fig. 5 (solid circles).
The proximity of these data points to $d^{2}E_{\textit{c}}/dVdV_{\textit{c}}$
is perhaps not entirely coincidental, but they sit rather far from their most
closely related quantity, $dN/dV$. The difference between the sum-over-states
expression and $dN/dV$ is subtle, since both are non-interacting response
functions with the same integration. It arises due to complex spatial
fluctuations in $v_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)$ both within
and without the weighted region in the case of $dN/dV$, in contrast to the
effective $\delta
v_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)=w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\delta
V$ implied in the sum-over-states expression. More precisely, $dN/dV$ is the
complete non-interacting linear-response function, whereas the sum-over-states
is the non-interacting linear-response function truncated at first-order in
perturbation theory.

Referring to Eq. III.2, the $dN/dV$ calculated using self-consistent cDFT
simplifies to the perturbation theory expression in the case that both
$dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)/dV^{\sigma}=w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)$
and the change in the average weighted Kohn-Sham potential $V^{\sigma}$ is the
wholly responsible for the change in subspace occupancy $N$, with all other
degrees of freedom fixed. The first condition and Eq. 17 together imply that

$\displaystyle 1\equiv{}$
$\displaystyle\int\left(\frac{dV^{\sigma}}{dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)}\right)\left(\frac{dv_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)}{dV^{\sigma}}\right)d\mathbf{r}$
$\displaystyle={}$ $\displaystyle\left(\int
w_{\textit{c}}^{\sigma}\left(\mathbf{r^{\prime}}\right)\;d\mathbf{r^{\prime}}\right)^{-1}\int\left(w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\right)^{2}\;d\mathbf{r},$
(45)

which would be readily satisfied only if
$w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)$ were an abrupt three
dimensional unit step function.

The second condition is more unrealistic, since we may directly vary the
Lagrange multiplier $V_{\textit{c}}$ when constraining DFT, but not, at least
directly, the average subspace Kohn-Sham potential $V^{\sigma}$. As
$V_{\textit{c}}$ is varied, self-consistent changes in
$v_{\sigma}^{\textrm{KS}}\left(\mathbf{r}\right)$ outside of the weighted
region may substantially mitigate the charge transfer due to changes within
it. Such effects are absent at first order in perturbation theory, and so it
may be expected to typically overestimate the magnitude of the non-interacting
response $dN/dV$, a situation which is exemplified in Fig. 5.

Since the differences in definition between the non-interacting linear and
sum-over-states response functions are related, albeit not entirely, to
effects outside of the weighted region, they cannot be reconciled. An
approximate reconciling renormalization may be made, however, by down-scaling
one of the two functions $w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)$ in
the unscreened analogue of Eq. III.2 (one applies the perturbation, the other
measures the charge) by a factor of $\sum_{\sigma}\left(\int
w^{\sigma}_{\textit{c}}\left(\mathbf{r}\right)\;d\mathbf{r}\right)/\left(\int\;d\mathbf{r^{\prime}}\right)$.
This scaling mimics the compensating background changes in the Kohn-Sham
potential in the realistic cDFT calculation, by reducing or “redistributing”
the change in potential $w_{\textit{c}}^{\sigma}\left(\mathbf{r}\right)\delta
V^{\sigma}$ in proportion to the system volume. The result of this
renormalization, suitably adapted for orbital-based population analysis (based
on orbital count, with the advantage that the rescaling is not extensive with
respect to the volume of the vacuum region), is shown in open black circles in
Fig. 5. We find that the unscreened sum-over-states expression is thus brought
into fair, but not exact, agreement with the non-perturbative non-interacting
response $dN/dV$. A further downscaling by a factor of the integrated inverse
microscopic dielectric function, given by Eq. IV.2, is then required to
approximately match the unscreened sum-over-states perturbation theory and
self-consistent cDFT total-energy curvatures.

## References

  * Kaduk  _et al._ (2012) B. Kaduk, T. Kowalczyk, and T. Van Voorhis, Chemical Reviews 112, 321 (2012).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Wu and Van Voorhis (2006a) Q. Wu and T. Van Voorhis, Journal of Chemical Theory and Computation 2, 765 (2006a).
  * Wu and Van Voorhis (2005) Q. Wu and T. Van Voorhis, Phys. Rev. A 72, 024502 (2005).
  * Segal  _et al._ (2007) M. Segal, M. Singh, K. Rivoire, S. Difley, T. Van Voorhis, and M. A. Baldo, Nature Materials 6, 374 (2007).
  * Kowalczyk  _et al._ (2010) T. Kowalczyk, Z. Lin, and T. V. Voorhis, The Journal of Physical Chemistry A 114, 10427 (2010).
  * Dederichs  _et al._ (1984) P. H. Dederichs, S. Blügel, R. Zeller, and H. Akai, Phys. Rev. Lett. 53, 2512 (1984).
  * Wu and Van Voorhis (2006b) Q. Wu and T. Van Voorhis, The Journal of Physical Chemistry A 110, 9212 (2006b).
  * Oberhofer and Blumberger (2009a) H. Oberhofer and J. Blumberger, The Journal of Chemical Physics 131, 064101 (2009a).
  * Sit  _et al._ (2006) P. H.-L. Sit, M. Cococcioni, and N. Marzari, Phys. Rev. Lett. 97, 028303 (2006).
  * Tim Kowalczyk  _et al._ (2011) T. Tim Kowalczyk, L.-P. Wang, and T. Van Voorhis, The Journal of Physical Chemistry B 115, 12135 (2011).
  * Behler  _et al._ (2008) J. Behler, K. Reuter, and M. Scheffler, Phys. Rev. B 77, 115421 (2008).
  * Rudra  _et al._ (2006) I. Rudra, Q. Wu, and T. Van Voorhis, The Journal of Chemical Physics 124, 024103 (2006).
  * Rudra  _et al._ (2007) I. Rudra, Q. Wu, and T. Van Voorhis, Inorganic Chemistry 46, 10539 (2007).
  * Schmidt  _et al._ (2008) J. R. Schmidt, N. Shenvi, and J. C. Tully, The Journal of Chemical Physics 129, 114110 (2008).
  * Difley  _et al._ (2008) S. Difley, D. Beljonne, and T. Troy Van Voorhis, Journal of the American Chemical Society 130, 3420 (2008).
  * Sena  _et al._ (2011) A. M. P. Sena, T. Miyazaki, and D. R. Bowler, Journal of Chemical Theory and Computation 7, 884 (2011).
  * Turban  _et al._ (2016) D. H. P. Turban, G. Teobaldi, D. D. O’Regan, and N. D. M. Hine, Phys. Rev. B 93, 165102 (2016).
  * Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * Difley and Van Voorhis (2011) S. Difley and T. Van Voorhis, Journal of Chemical Theory and Computation 7, 594 (2011).
  * Nakamura  _et al._ (2008) K. Nakamura, Y. Yoshimoto, R. Arita, S. Tsuneyuki, and M. Imada, Phys. Rev. B 77, 195126 (2008).
  * Roychoudhury  _et al._ (2016) S. Roychoudhury, C. Motta, and S. Sanvito, Phys. Rev. B 93, 045130 (2016).
  * Souza  _et al._ (2013) A. M. Souza, I. Rungger, C. D. Pemmaraju, U. Schwingenschloegl, and S. Sanvito, Phys. Rev. B 88, 165112 (2013).
  * Sau  _et al._ (2008) J. D. Sau, J. B. Neaton, H. J. Choi, S. G. Louie, and M. L. Cohen, Phys. Rev. Lett. 101, 026804 (2008).
  * Brooke  _et al._ (2015) C. Brooke, A. Vezzoli, S. J. Higgins, L. A. Zotti, J. J. Palacios, and R. J. Nichols, Phys. Rev. B 91, 195438 (2015).
  * Wu and Van Voorhis (2006c) Q. Wu and T. Van Voorhis, The Journal of Chemical Physics 125, 164105 (2006c).
  * Ding  _et al._ (2010) F. Ding, H. Wang, Q. Wu, T. Van Voorhis, S. Chen, and J. P. Konopelski, The Journal of Physical Chemistry A 114, 6039 (2010).
  * Oberhofer and Blumberger (2009b) H. Oberhofer and J. Blumberger, J. Chem. Phys. 131, 064101 (2009b).
  * Oberhofer and Blumberger (2010a) H. Oberhofer and J. Blumberger, J. Chem. Phys. 133, 244105 (2010a).
  * Kubas  _et al._ (2014) A. Kubas, F. Hoffmann, A. Heck, H. Oberhofer, M. Elstner, and J. Blumberger, J. Chem. Phys. 140, 104105 (2014).
  * Kubas  _et al._ (2015) A. Kubas, F. Gajdos, A. Heck, H. Oberhofer, M. Elstner, and J. Blumberger, Phys. Chem. Chem. Phys. 17, 14342 (2015).
  * Oberhofer and Blumberger (2010b) H. Oberhofer and J. Blumberger, Angew. Chem. Int. Ed. 49, 3631 (2010b).
  * Oberhofer and Blumberger (2012) H. Oberhofer and J. Blumberger, Phys. Chem. Chem. Phys. 14, 13846 (2012).
  * McKenna and Blumberger (2012) K. P. McKenna and J. Blumberger, Phys. Rev. B 86, 245110 (2012).
  * Blumberger and McKenna (2013) J. Blumberger and K. P. McKenna, Phys. Chem. Chem. Phys. 15, 2184 (2013).
  * Yeganeh and Van Voorhis (2010) S. Yeganeh and T. Van Voorhis, The Journal of Physical Chemistry C 114, 20756 (2010).
  * Yost  _et al._ (2014) S. R. Yost, J. Lee, M. W. B. Wilson, T. Wu, D. P. McMahon, R. R. Parkhurst, N. J. Thompson, D. N. Congreve, A. Rao, K. Johnson, M. Y. Sfeir, M. G. Bawendi, T. M. Swager, R. H. Friend, M. A. Baldo, and T. Van Voorhis, Nat. Chem. 6, 492 (2014).
  * Evans  _et al._ (2008) J. S. Evans, C.-L. Cheng, and T. Van Voorhis, Phys. Rev. B 78, 165108 (2008).
  * Anisimov  _et al._ (1991) V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
  * Meider and Springborg (1998) H. Meider and M. Springborg, Journal of Physics: Condensed Matter 10, 6953 (1998).
  * Danielsen (1986) P. L. Danielsen, Journal of Physics C: Solid State Physics 19, L741 (1986).
  * Meider and Springborg (1999) H. Meider and M. Springborg, Chemical Physics Letters 300, 339 (1999).
  * Hybertsen  _et al._ (1989) M. S. Hybertsen, M. Schlüter, and N. E. Christensen, Phys. Rev. B 39, 9028 (1989).
  * Petukhov  _et al._ (2003) A. G. Petukhov, I. I. Mazin, L. Chioncel, and A. I. Lichtenstein, Phys. Rev. B 67, 153106 (2003).
  * Wu  _et al._ (2007) Q. Wu, C.-L. Cheng, and T. Van Voorhis, The Journal of Chemical Physics 127, 164119 (2007).
  * Wu  _et al._ (2009) Q. Wu, B. Kaduk, and T. Van Voorhis, The Journal of Chemical Physics 130, 034109 (2009).
  * Kaduk and Van Voorhis (2010) B. Kaduk and T. Van Voorhis, The Journal of Chemical Physics 133, 061102 (2010).
  * Kaduk  _et al._ (2014) B. Kaduk, T. Tsuchimochi, and T. Van Voorhis, The Journal of Chemical Physics 140 (2014).
  * Jain  _et al._ (2011) A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fischer, T. Mueller, K. A. Persson, and G. Ceder, Computational Materials Science 50, 2295 (2011).
  * Curtarolo  _et al._ (2013) S. Curtarolo, G. L. W. Hart, M. B. Nardelli, N. Mingo, S. Sanvito, and O. Levy, Nat Mater 12, 191 (2013).
  * R̆ezác̆  _et al._ (2012) J. R̆ezác̆, B. Lévy, I. Demachy, and A. de la Lande, Journal of Chemical Theory and Computation 8, 418 (2012).
  * Ratcliff  _et al._ (2015) L. E. Ratcliff, L. Genovese, S. Mohr, and T. Deutsch, The Journal of Chemical Physics 142, 234105 (2015).
  * Note (1) The “curvature” is used here as a convenient shorthand for the second derivative. We do not imply the geometric curvature, which equals the second derivative only at stationary points.
  * Wu (2016) Q. Wu, personal communication (2016).
  * Wu and Yang (2003) Q. Wu and W. Yang, The Journal of Chemical Physics 118, 2498 (2003).
  * Skylaris  _et al._ (2005) C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C. Payne, The Journal of Chemical Physics 122, 084119 (2005).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Coulson and Fischer (1949) C. Coulson and I. Fischer, Philos. Mag. 40, 386 (1949).
  * Ruzsinszky  _et al._ (2005) A. Ruzsinszky, J. P. Perdew, and G. I. Csonka, The Journal of Physical Chemistry A 109, 11006 (2005).
  * Dziedzic  _et al._ (2011) J. Dziedzic, H. H. Helal, C.-K. Skylaris, A. A. Mostofi, and M. C. Payne, EPL (Europhysics Letters) 95, 43001 (2011).
  * Andreussi  _et al._ (2012) O. Andreussi, I. Dabo, and N. Marzari, The Journal of Chemical Physics 136, 064102 (2012).
  * Note (2) A single, strictly local occupancy constraint is considered in Ref. PhysRevA.72.024502. We retain these restrictions so as not to obscure the fundamental aspects under consideration. These conditions are typically lifted in practical cDFT calculations, bringing us into multivariate optimisation of constrained Kohn-Sham spin-density functional theory, which may also be non-local or orbital dependent.
  * Baroni  _et al._ (1987) S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58, 1861 (1987).
  * Baroni  _et al._ (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Note (3) We use a lunate epsilon $\epsilon$ for the microscopic dielectric function in order to distinguish from it from the Kohn-Sham eigenvalues $\varepsilon$.
  * Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  * Lee  _et al._ (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  * Skylaris  _et al._ (2002) C.-K. Skylaris, A. A. Mostofi, P. D. Haynes, O. Diéguez, and M. C. Payne, Phys. Rev. B 66, 035119 (2002).
  * O’Regan  _et al._ (2012) D. D. O’Regan, M. C. Payne, and A. A. Mostofi, Phys. Rev. B 85, 193101 (2012).
  * Note (4) The resulting constraint acts on the Kohn-Sham density-matrix rather than on the density. Our analytical findings extends to that case with minor notational changes.
  * Ratcliff  _et al._ (2011) L. E. Ratcliff, N. D. M. Hine, and P. D. Haynes, Phys. Rev. B 84, 165131 (2011).
  * Rondinelli  _et al._ (2007) J. M. Rondinelli, B. Deng, and L. D. Marks, Computational Materials Science 40, 345 (2007).
  * Harville (2011) D. Harville,  _Matrix Algebra: Exercises and Solutions_ (Springer New York, 2011).
  * Valiev  _et al._ (2010) M. Valiev, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, H. V. Dam, D. Wang, J. Nieplocha, E. Apra, T. Windus, and W. de Jong, Computer Physics Communications 181, 1477 (2010).
  * Shao  _et al._ (2015) Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kús, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock III, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. DiStasio Jr., H. Dop, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. D. Hanson-Heine, P. H. P. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.-P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, C. M. Oana, R. Olivares-Amaya, D. P. O’Neill, J. A. Parkhill, T. M. Perrine, R. Peverati, P. A. Pieniazek, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, N. Sergueev, S. M. Sharada, S. Sharmaa, D. W. Small, A. Sodt, T. Stein, D. Stück, Y.-C. Su, A. J. W. Thom, T. Tsuchimochi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, V. Vanovschi, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhou, B. R. Brooks, G. K. L. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard III, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer III, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xua, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. Van Voorhis, J. M. Herbert, A. I. Krylov, P. M. W. Gill, and M. Head-Gordon, Mol. Phys. 113, 184 (2015).
  * Kurz  _et al._ (2004) P. Kurz, F. Förster, L. Nordström, G. Bihlmayer, and S. Blügel, Phys. Rev. B 69, 024415 (2004).
  * Ozaki and Kino (2005) T. Ozaki and H. Kino, Phys. Rev. B 72, 045121 (2005).
  * Hutter  _et al._ (2014) J. Hutter, M. Iannuzzi, F. Schiffmann, and J. VandeVondele, Wiley Interdisciplinary Reviews: Computational Molecular Science 4, 15 (2014).
  * Ma and Dudarev (2015) P.-W. Ma and S. L. Dudarev, Phys. Rev. B 91, 054420 (2015).
  * O’Regan  _et al._ (2010) D. D. O’Regan, N. D. M. Hine, M. C. Payne, and A. A. Mostofi, Phys. Rev. B 82, 081102 (2010).