# Exciton ground-state energy with full hole warping structure

Roland Combescota and Shiue-Yuan Shiaub a Laboratoire de Physique, Ecole
Normale Supérieure, PSL Université, Sorbonne Université, Paris Diderot
Université, CNRS, 24 rue Lhomond, F-75005 Paris, France.
bPhysics Division, National Center for Theoretical Sciences, 10617,Taipei,
Taiwan.

(Received )

###### Abstract

Most semiconductors, in particular III-V compounds, have a complex valence
band structure near the band edge, due to degeneracy at the zone center. One
peculiar feature is the warping of the electronic dispersion relations, which
are not isotropic even in the vicinity of the band edge. When the exciton, all
important for the semiconductor optical properties, is considered, this
problem is usually handled by using some kind of angular averaging procedure,
that would restore the isotropy of the hole effective dispersion relations. In
the present paper, we consider the problem of the exciton ground-state energy
for semiconductors with zinc-blende crystal structure, and we solve it exactly
by a numerical treatment, taking fully into account the warping of the valence
band. In the resulting four-dimensional problem, we first show exactly that
the exciton ground state is fourfold degenerate. We then explore the ground-
state energy across the full range of allowed Luttinger parameters. We find
that the correction due to warping may in principle be quite large. However,
for the semiconductors with available data for the band structure we have
considered, the correction turns out to be in the $10\%-15\%$ range.

###### pacs:

03.65.-w , 31.15.-p , 71.35.-y

##  I INTRODUCTION

Optical properties of semiconductors are of utmost interest both for
fundamental and applied purposes cardona ; vurga . With respect to these
properties, excitonic excitations play a prominent role, especially in the
visible light spectrum. Hence, knowing the exciton ground-state properties in
these compounds is of fundamental importance. However, this is in general not
a simple matter theoretically, since generically the valence band of the most
relevant semiconductors is degenerate at the band edge. In the standard case
of semiconductors with cubic zinc-blende crystal structure, the states at the
valence band edge have a p-wave character which leads to a sixfold degeneracy,
when spin is taken into account. When spin-orbit coupling is strong enough to
make a twofold degenerate split-off subband far enough in energy from the
valence band edge, so that it becomes irrelevant to the excitonic structure,
one is left with a fourfold degeneracy. This is the situation we will
consider.

Unfortunately, it has been long recognized that, even in the vicinity of the
band edge, the resulting valence band dispersion relations are not isotropic,
while respecting the cubic symmetry of the crystal. As a result, the equal
energy surfaces are warped, and this whole matter makes the hole properties
difficult to handle. This complexity is usually disposed of by taking an
appropriate angular average to go back to an isotropic situation. In the
present paper, we will not make use of such a simplification, but rather
explore the exciton ground-state energy in these zinc-blende semiconductors in
full generality.

In the next section, we recall the hole effective kinetic energy Hamiltonian,
standard for these compounds. In the following section III, we write the
Schrödinger equation for the exciton under the form of an integral equation
appropriate for our study. We show in section IV in full generality that, due
to the cubic symmetry, the exciton ground state has a fourfold degeneracy
(omitting the twofold additional degeneracy coming from the conduction
electronic spin). Finally, in section V we solve numerically our integral
equation. We show that in principle the warping leads to an increase of the
exciton binding energy, which may be quite important. However, in practice,
for the tabulated semiconductors we have found, this increase is limited to a
$10\%-15\%$ range. Nevertheless, markedly higher values are possible,
considering the imprecision of the known valence band parameters data, or
other semiconductors yet to be investigated.

##  II Formalism

In the semiconductor compounds we are considering, the bottom of the
conduction band located at the zone center is non-degenerate (except for the
spin degree of freedom), but the top of the valence band located at the zone
center is, in the absence of spin-orbit coupling, threefold degenerate
(forgetting the hole spin), with a p-wave-like character for the corresponding
wave functions. As the exciton wave function is predominantly made of electron
and hole states in the vicinity of the zone center, this degeneracy has to be
taken into account in the exciton wave function. Taking spin into account,
this threefold degeneracy becomes a sixfold degeneracy. When spin-orbit
coupling is taken into account, this degeneracy is lifted into a fourfold
degeneracy, plus a twofold one. In standard compounds, these last two
degenerate electronic states are shifted to energy low enough, so that they
are irrelevant to the building of the exciton. This is the situation we will
restrict ourselves to in the following. The four-dimensional hole subspace we
now focus on has the symmetry character of a $J=3/2$ angular momentum lutt ,
resulting from the composition of the p-wave-like wave functions and of the
$1/2$ spin. For the present purpose, it is more convenient to use the
notations $|1\rangle\equiv|3/2\rangle$, $|2\rangle\equiv|1/2\rangle$,
$|3\rangle\equiv|-1/2\rangle$ and $|4\rangle\equiv|-3/2\rangle$ to denote the
standard fourfold hole basis.

Following the approach of Luttinger and Kohn lutkoh ; lutt  (with slightly
different notations note1 ) in the absence of magnetic field, the hole kinetic
energy part $H_{h}({\bf k})$ of the exciton Hamiltonian is given by

$\displaystyle H_{h}({\bf
k})=\frac{\hbar^{2}}{2m_{0}}\left(\begin{array}[]{cccc}G+\Delta&D&F&0\\\
D^{*}&G-\Delta&0&F\\\ F^{*}&0&G-\Delta&-D\\\
0&F^{*}&-D^{*}&G+\Delta\end{array}\right)\equiv\frac{\hbar^{2}}{2m_{0}}(G\,\mathbb{1}+h_{h})$
(5)

where $m_{0}$ is the vacuum electron mass (providing a mass scale), and
$G,\Delta,D$ and $F$ are defined from the hole wavevector ${\bf k}$ by

$\displaystyle G$ $\displaystyle=$ $\displaystyle\gamma_{1}k^{2}\hskip
42.67912pt\Delta=\gamma_{2}(k^{2}-3k_{z}^{2})$ (6) $\displaystyle D$
$\displaystyle=$ $\displaystyle-2\sqrt{3}\gamma_{3}(k_{x}-ik_{y})k_{z}$
$\displaystyle F$ $\displaystyle=$
$\displaystyle-\sqrt{3}\gamma_{2}(k_{x}^{2}-k_{y}^{2})+2i\sqrt{3}\gamma_{3}k_{x}k_{y}$

where $x,y$ and $z$ are the cubic axes, and $\gamma_{1}$, $\gamma_{2}$ and
$\gamma_{3}$ are the Luttinger lutt  constants characterizing the
semiconductor.

One may notice that, within a multiplicative factor, $h_{h}^{2}$ reduces to
the unit matrix

$\displaystyle\left(\frac{2m_{0}}{\hbar^{2}}H_{h}({\bf
k})-G\,\mathbb{1}\right)^{2}=h_{h}^{2}=\lambda^{2}\,\mathbb{1}$ (7)

with

$\displaystyle\lambda^{2}\equiv\Delta^{2}+|D|^{2}+|F|^{2}=4\gamma_{2}^{2}k^{4}+12(\gamma_{3}^{2}-\gamma_{2}^{2})\left[k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+k_{z}^{2}k_{x}^{2}\right]$
(8)

so that the doubly degenerate eigenvalues $E_{\pm}({\bf k})$ of $H_{h}({\bf
k})$ are given by the well-known result dkk

$\displaystyle\frac{2m_{0}}{\hbar^{2}}E_{\pm}({\bf k})=\gamma_{1}k^{2}\pm
2\Big{[}\gamma_{2}^{2}k^{4}+3(\gamma_{3}^{2}-\gamma_{2}^{2})\left[k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+k_{z}^{2}k_{x}^{2}\right]\Big{]}^{1/2}$
(9)

The square root on the right-hand side gives rise to the well-known warping,
implying in general a departure from a dispersion relation with simple
spherical symmetry. This symmetry is recovered only in the particular case
where $\gamma_{2}=\pm\gamma_{3}$. In this specific case, one finds the
standard heavy and light holes, with mass $m_{H}$ and $m_{L}$, related to our
coefficients by (assuming $\gamma_{2}>0$)

$\displaystyle\frac{m_{0}}{m_{H}}=\gamma_{1}-2\gamma_{2}\hskip
56.9055pt\frac{m_{0}}{m_{L}}=\gamma_{1}+2\gamma_{2}$ (10)

In order to make sense as a hole dispersion relation, $E_{\pm}({\bf k})$ which
gives the hole kinetic energy should naturally be positive for any ${\bf k}$.
This implies some limitations on the Luttinger coefficients, from the explicit
expression given by Eq.(9). Indeed, considering the negative branch
$E_{-}({\bf k})$, this requires the square root on the right-hand side to be
less than $\gamma_{1}k^{2}$. If we consider first the case where
$\gamma_{3}<\gamma_{2}$, the square root is maximal when the bracket
multiplying $\gamma_{3}^{2}-\gamma_{2}^{2}$ is zero, which occurs when
$k_{x}=k_{y}=0$, or $k_{y}=k_{z}=0$, or $k_{z}=k_{x}=0$, that is, when ${\bf
k}$ is on one of the cubic axes. In this case, the positivity condition is
merely

$\displaystyle\gamma_{1}>2\gamma_{2}$ (11)

In the opposite case where $\gamma_{3}>\gamma_{2}$, in order to maximize the
square root, we instead look for the maximal value of the bracket multiplying
$(\gamma_{3}^{2}-\gamma_{2}^{2})$. For fixed $k$, it is found when $\pm
k_{x}=\pm k_{y}=\pm k_{z}=k/\sqrt{3}$, that is, when ${\bf k}$ is along one of
the principal cubic diagonals. In this case, the square root is merely equal
to $\gamma_{3}k^{2}$, and the positivity condition is (assuming
$\gamma_{3}>0$)

$\displaystyle\gamma_{1}>2\gamma_{3}$ (12)

In the particular case where $\gamma_{3}=\gamma_{2}$, these conditions merely
reduce to the fact that $m_{H}$ should be positive, which is physically
obvious.

##  III Integral equation for excitonic eigenstates

The Hamiltonian for the exciton is merely obtained by adding to the hole
kinetic energy the electronic kinetic energy ${\bf p}_{e}^{2}/2m_{e}$, where
${\bf p}_{e}$ is the conduction band electronic momentum, and $m_{e}$ the
corresponding electronic band mass. We assume that the electronic dispersion
relation is isotropic ross . The electronic spin does not bring any
complication, except that naturally all the degeneracies are multiplied by
two; so, we do not indicate it explicitly. Finally, we have to include in the
Hamiltonian the electron-hole attractive Coulomb interaction $V({\bf
r}_{e}-{\bf r}_{h})$, where ${\bf r}_{e}$ and ${\bf r}_{h}$ are respectively
the electron and hole position, with $V({\bf r})=-e^{2}/(4\pi\epsilon_{sc}r)$
where $\epsilon_{sc}$ is the semiconductor dielectric constant. Since this
interaction depends only on the relative position of the hole and the
electron, it is translationally invariant, so the excitonic momentum is
conserved. In the present paper, we restrict ourselves to the case where this
momentum is zero, so that the electron and hole momentum are opposite ${\bf
p}_{e}=-{\bf p}_{h}=-\hbar{\bf k}$. Hence, the electron term merely adds a
contribution $\hbar^{2}k^{2}/2m_{e}$ to the hole kinetic energy Eq.(5), and we
are back to a one-body problem. The potential term is the Coulomb interaction,
and the kinetic energy is still given by Eq.(5), provided $\gamma_{1}$ is
replaced by

$\displaystyle{\bar{\gamma}}_{1}=\gamma_{1}+\frac{m_{0}}{m_{e}}$ (13)

Considering the complexity of the hole kinetic energy, it is more convenient
to write the excitonic Schrödinger equation in momentum space, which leads to
an integral equation. This is much easier to handle numerically than the four-
dimensional partial differential equation one would find if one worked in real
space. We denote the ${\bf k}$ representation of the four-components wave
function as $a({\bf k})\equiv\big{(}a_{1}({\bf k}),a_{2}({\bf k}),a_{3}({\bf
k}),a_{4}({\bf k})\big{)}$. This four-dimensional Schrödinger equation then
reads

$\displaystyle{\bar{H}}_{h}({\bf k})\,a({\bf k})+\int\frac{d{\bf
k}^{\prime}}{(2\pi)^{3}}\,V({\bf k}-{\bf k}^{\prime})\,a({\bf
k}^{\prime})=E\,a({\bf k})$ (14)

where $V({\bf q})=-e^{2}/(\epsilon_{sc}q^{2})$ is the Fourier transform of the
Coulomb interaction, and $E$ is the exciton energy to be solved for.
${\bar{H}}_{h}({\bf k})$ is just given by Eq.(5), except that $\gamma_{1}$ has
to be replaced by ${\bar{\gamma}}_{1}$. Note that the Coulomb interaction is
diagonal in our representation, while ${\bar{H}}_{h}({\bf k})$ is not.

It is convenient to rewrite this equation so that only dimensionless
quantities appear. The natural energy unit is the Rydberg, so we write the
exciton energy as
$E=-{\bar{E}}\,m_{0}(e^{2}/4\pi\epsilon_{sc})^{2}/(2\hbar^{2})$ where
${\bar{E}}$ is the dimensionless exciton binding energy. The corresponding
natural length unit is the Bohr radius, so we also go to a dimensionless
wavevector ${\bar{\bf k}}$ defined by ${\bf k}={\bar{\bf
k}}\,m_{0}(e^{2}/4\pi\epsilon_{sc})/\hbar^{2}$. Taking into account that
${\bar{H}}_{h}({\bf k})$ is a quadratic function of the wavevector, our
Schrödinger equation becomes

$\displaystyle\big{(}{\bar{H}}_{r}({\bar{\bf k}})+{\bar{E}}\big{)}a({\bar{\bf
k}})=\frac{1}{\pi^{2}}\int d{\bar{\bf k}}^{\prime}\frac{1}{({\bar{\bf
k}}-{\bar{\bf k}}^{\prime})^{2}}\,a({\bar{\bf k}}^{\prime})$ (15)

where we have also defined a reduced dimensionless hole kinetic energy
${\bar{H}}_{r}({\bar{\bf k}})=(2m_{0}/\hbar^{2}){\bar{H}}_{h}({\bar{\bf k}})$.

Furthermore, it is clear from Eq.(5) and Eq.(13) that the important parameters
are the ratios $\gamma_{2}/{\bar{\gamma}}_{1}$ and
$\gamma_{3}/{\bar{\gamma}}_{1}$. Setting

$\displaystyle{\bar{\gamma}_{2}}=\frac{\gamma_{2}}{{\bar{\gamma}}_{1}}\hskip
56.9055pt{\bar{\gamma}_{3}}=\frac{\gamma_{3}}{{\bar{\gamma}}_{1}}$ (16)

it is convenient to have only these parameters appearing. This requires a
further rescaling, obtained by setting ${\bar{\bf k}}=K/{{\bar{\gamma}}_{1}}$
and ${\bar{E}}={\bar{\epsilon}}/{{\bar{\gamma}}_{1}}$. This leads to

$\displaystyle\big{(}(K^{2}+{\bar{\epsilon}})\mathbb{1}+{\bar{h}}_{h}({\bf
K})\big{)}a({\bf K})=\frac{1}{\pi^{2}}\int d{\bf K}^{\prime}\frac{1}{({\bf
K}-{\bf K}^{\prime})^{2}}\,a({\bf K}^{\prime})$ (17)

where ${\bar{h}}_{h}({\bf K})$ has exactly the same expression as $h_{h}({\bf
k})$ in Eq.(5) (with ${\bf k}$ replaced by ${\bf K}$) except that, in the
explicit expressions for its elements given in Eq.(6), one has to replace
$\gamma_{2}$ and $\gamma_{3}$ by ${\bar{\gamma}_{2}}$ and ${\bar{\gamma}_{3}}$
respectively.

We note that, corresponding to our above request for the positivity of the
hole kinetic energy, we have to require that the kinetic energy for the
exciton one-body problem is always positive. Proceeding as in the above
derivation of Eq.(11) and Eq.(12), we find that our coefficients must satisfy

$\displaystyle|{\bar{\gamma}_{2}}|<\frac{1}{2}\hskip
56.9055pt|{\bar{\gamma}_{3}}|<\frac{1}{2}$ (18)

For numerical work, it is more convenient to apply to Eq.(17) the inverse of
the $4\times 4$ matrix $(K^{2}+{\bar{\epsilon}})\mathbb{1}+{\bar{h}}_{h}({\bf
K})$. Since we have noted in Eq.(7) and Eq.(8) that
${\bar{h}}_{h}^{2}=\lambda^{2}\mathbb{1}$ is proportional to the unit matrix
(with $\gamma_{2,3}$ and ${\bf k}$ replaced in $\lambda$ by
${\bar{\gamma}_{2,3}}$ and ${\bf K}$), this inverse is merely
$\big{(}(K^{2}+{\bar{\epsilon}})\mathbb{1}-{\bar{h}}_{h}({\bf
K})\big{)}/((K^{2}+{\bar{\epsilon}})^{2}-\lambda^{2})$. Hence Eq.(17) becomes

$\displaystyle a({\bf
K})=\frac{(K^{2}+{\bar{\epsilon}})\mathbb{1}-{\bar{h}}_{h}({\bf
K})}{(K^{2}+{\bar{\epsilon}})^{2}-\lambda^{2}}\frac{1}{\pi^{2}}\int d{\bf
K}^{\prime}\frac{1}{({\bf K}-{\bf K}^{\prime})^{2}}\,a({\bf K}^{\prime})$ (19)

Note that the conditions Eq.(18) correspond to the requirement that the
denominator $(K^{2}+{\bar{\epsilon}})^{2}-\lambda^{2}$ is always positive, in
particular for large $K$.

To conclude this section, it is useful to consider the particular case where
our exciton problem reduces to the case of the hydrogen atom. This corresponds
merely to the case where ${\bar{h}}_{h}({\bf K})=0$ and the integral equation
(19) reduces to a one-dimensional equation

$\displaystyle a({\bf K})=\frac{1}{\pi^{2}(K^{2}+{\bar{\epsilon}})}\int d{\bf
K}^{\prime}\frac{1}{({\bf K}-{\bf K}^{\prime})^{2}}\,a({\bf K}^{\prime})$ (20)

We readily know the solution which is the Fourier transform of the ground-
state wave function $\exp(-r/a_{0})$ of the hydrogen atom, that is (omitting
the unimportant constant prefactor) $a({\bf K})=1/(1+K^{2})^{2}$. Taking into
account that, with our reduced units, the ground-state energy is merely
${\bar{\epsilon}}=1$, one can easily check analytically that this expression
of $a({\bf K})$ is indeed solution of the integral equation (20).

It is important to note that the explicit large $K$ behaviour $a({\bf K})\sim
1/K^{4}$ of this solution, which insures the large $K$ convergence of the
integral on the right-hand side, is actually a generic feature of the general
equation (19). Indeed, making use of this convergence consistently allows us
to write in this large $K$ regime $1/({\bf K}-{\bf K}^{\prime})^{2}\simeq
1/K^{2}$ in the integral. Then, the resulting $1/K^{2}$ factor, together with
the explicit prefactor of the integral, indeed leads to $a({\bf K})\sim
1/K^{4}$ for the large $K$ behaviour of $a({\bf K})$. This behaviour is
naturally of particular interest for the numerical solution of this integral
equation.

##  IV Ground state degeneracy

One might think at first that Eq.(19) gives for the exciton ground state a
single non-degenerate solution, with a specific symmetry. This is not what
happens. Instead, we find that quite generally the exciton ground state, at
zero total momentum, has an exact fourfold degeneracy (without taking into
account the conduction electron spin degeneracy). This comes directly from the
cubic symmetry of the crystal. This is fully analogous to what would occur if
we had a complete rotational invariance. In this last case, this would
naturally imply a full degeneracy, with the eigenstates having in this four-
dimensional subspace an angular dependence given by the spherical harmonics
$Y_{\ell m}(\theta,\varphi)$ with $\ell=3/2$ and $m=(\pm 3/2,\pm 1/2)$. In the
present case, we have actually the lower cubic symmetry, but it is enough to
similarly insure the complete degeneracy of the ground state.

This result may at first look surprising if one naively considers the
situation with spherical symmetry ${\bar{\gamma}_{2}}={\bar{\gamma}_{3}}$.
Since in this case one has heavy and light holes, one could think of two
possible exciton states, made respectively with a heavy or a light hole (each
one being doubly degenerate). But this omits the fact that the Coulomb
interaction is not diagonal in this heavy-light hole representation. Indeed,
it has been shown monsean  that the proper treatment of this specific case
leads to a fourfold degeneracy. Here, we extend this result to the general
case.

To see this in a more specific way, we have to consider the symmetry
properties of the solutions. We rewrite the Schrödinger equation Eq.(17) for
the exciton, making explicit the wave function and matrix components. It
reads, with ($m,n)=1,2,3,4$

$\displaystyle\sum_{n=1}^{4}\big{(}(K^{2}+{\bar{\epsilon}})\delta_{mn}+{\bar{h}}_{mn}({\bf
K})\big{)}a_{n}({\bf K})=\frac{1}{\pi^{2}}\int d{\bf K}^{\prime}\frac{1}{({\bf
K}-{\bf K}^{\prime})^{2}}\,a_{m}({\bf K}^{\prime})$ (21)

where ${\bar{h}}_{mn}({\bf K})$ are the matrix elements of ${\bar{h}}_{h}({\bf
K})$ and $\delta_{mn}$ is the Kronecker symbol. The fact that this equation
satisfies time-reversal invariance implies that, if $a_{n}({\bf K})$ is
solution of this equation, $(-1)^{n-1}a^{*}_{5-n}({\bf K})$ is also solution
with the same energy. This transformation can be understood from the
spin-$3/2$ nature of our components. This property can be checked directly
from the explicit expression of the hole kinetic energy terms
${\bar{h}}_{mn}({\bf K})$, given by Eq.(6).

Let us first consider the symmetry properties under a $\pi/2$ rotation $R_{z}$
around our quantization axis $z$. Under this rotation, the component of the
wave function $a_{n}({\bf K})$ becomes $a_{n}(R_{z}^{-1}({\bf K}))$, that is,
$a_{n}(K_{x},K_{y},K_{z})$ becomes $a_{n}(K_{y},-K_{x},K_{z})$. For clarity,
we denote as ${\bar{a}}_{n}$ the transform of $a_{n}$, that is,
${\bar{a}}_{n}(K_{x},K_{y},K_{z})=a_{n}(K_{y},-K_{x},K_{z})$. Similarly, if we
perform an additional $\pi/2$ rotation $R_{z}$ on the wave function, amounting
to a total $\pi$ rotation, we end up with
${\bar{\bar{a}}}_{n}(K_{x},K_{y},K_{z})=a_{n}(-K_{x},-K_{y},K_{z})$. Finally,
an additional $\pi/2$ rotation leads to
${\bar{\bar{\bar{a}}}}_{n}(K_{x},K_{y},K_{z})=a_{n}(-K_{y},K_{x},K_{z})$.

However, if we want to rotate the state $|a\rangle=\sum_{n}a_{n}|n\rangle$ by
a $\pi/2$ rotation $R_{z}$, we have to take into account that our basis states
$|3/2\rangle\equiv|1\rangle$, $|1/2\rangle\equiv|2\rangle$,
$|-1/2\rangle\equiv|3\rangle$ and $|-3/2\rangle\equiv|4\rangle$ change under
rotation by acquiring a phase factor through the action of the rotation
operator $e^{-i\pi J_{z}/2\hbar}$. Disregarding the unimportant overall phase
factor $e^{-3i\pi/4}$, this leads for the various components $a_{n}$ of our
wave function to an additional phase factor $e^{i(n-1)\pi/2}=i^{n-1}$ coming
from our basis states.

If we now make in Eq.(21) the change of variables $K_{x}\to K_{y}$ and
$K_{y}\to-K_{x}$ (and similarly for ${\bf K}^{\prime}$ on the right-hand
side), $a_{n}({\bf K})$ becomes ${\bar{a}}_{n}(\bf K)$ while we see from their
explicit expressions Eq.(6) that $D({\bf K})$ becomes $i\,D({\bf K})$, $F({\bf
K})$ becomes $-F({\bf K})$ and $\Delta({\bf K})$ is unchanged. Comparing with
the original equation (21) before the change of variables, one finds the
following symmetry property: if $a_{n}$ is solution of Eq.(21), then
$i^{n-1}\,{\bar{a}}_{n}$ is also solution of this equation with the same
energy. That is, if we have a solution, the wave function obtained by a
$\pi/2$ rotation is also a solution. This physically corresponds to the fact
that our Hamiltonian is invariant under a $\pi/2$ rotation $R_{z}$. Naturally,
repeating the $\pi/2$ rotation, we find that $(-1)^{n-1}\,{\bar{\bar{a}}}_{n}$
and $(-i)^{n-1}\,{\bar{\bar{\bar{a}}}}_{n}$ are also solutions.

This situation is analogous to the simple one-dimensional case where the
Hamiltonian is invariant under the transformation $x\to-x$: if $\psi(x)$ is a
solution of the Schrödinger equation, $\psi(-x)$ is also a solution with the
same energy. So, we can build the two solutions
$\psi_{S}(x)=\psi(x)+\psi(-x)$, and $\psi_{A}(x)=\psi(x)-\psi(-x)$ which are
respectively symmetric and antisymmetric under $x\to-x$.

In the same way, in our case, if $a_{n}({\bf K})$ is a solution of the
Schrödinger equation (21), then by adding the solutions obtained by performing
rotations around the $z$ axis by $\pi/2$, $\pi$ and $3\pi/2$ respectively, we
obtain a solution with the same energy which is invariant by $\pi/2$ rotation.
This solution is explicitly

$\displaystyle
a^{(1)}_{n}=a_{n}+i^{n-1}\,{\bar{a}}_{n}+(-1)^{n-1}\,{\bar{\bar{a}}}_{n}+(-i)^{n-1}\,{\bar{\bar{\bar{a}}}}_{n}$
(22)

If we perform on this solution the change of variables $K_{x}\to K_{y}$ and
$K_{y}\to-K_{x}$, we find that $a^{(1)}_{n}$ is changed into
$(-i)^{n-1}\,a^{(1)}_{n}$, that is
${\bar{a}}^{(1)}_{n}=(-i)^{n-1}\,a^{(1)}_{n}$.

But, as possible solutions of Eq.(21), we can more generally take solutions
$a^{(p)}_{n}$, with $p=1,2,3,4$, which transform as
${\bar{a}}^{(p)}_{n}=(-i)^{n-p}\,a^{(p)}_{n}$. Generalizing the preceding case
$p=1$, they are obtained explicitly from a solution $a_{n}({\bf K})$ by

$\displaystyle
a^{(p)}_{n}=a_{n}+i^{n-p}\,{\bar{a}}_{n}+(-1)^{n-p}\,{\bar{\bar{a}}}_{n}+(-i)^{n-p}\,{\bar{\bar{\bar{a}}}}_{n}$
(23)

Performing the changes $K_{x}\to K_{y}$ and $K_{y}\to-K_{x}$, we check as for
Eq.(22) that they indeed transform as we have indicated. These are the
appropriate generalizations to the present case of choosing solutions which
are either even or odd under $x\to-x$, when the Hamiltonian is invariant under
this transformation $x\to-x$. These four solutions are orthogonal since the
scalar product of two different solutions is naturally unchanged when one
performs on ${\bf K}$ in the integral the change of variables $K_{x}\to K_{y}$
and $K_{y}\to-K_{x}$ corresponding to the $\pi/2$ rotation, while it is
changed due to the different way the two different solutions with $p\neq q$
transform under the rotation

$\displaystyle\sum_{n=1}^{4}\int d{\bf K}\,a^{(p)\,*}_{n}({\bf
K})a^{(q)}_{n}({\bf K})=\sum_{n=1}^{4}\int d{\bf
K}\,{\bar{a}}^{(p)\,*}_{n}({\bf K}){\bar{a}}^{(q)}_{n}({\bf
K})=i^{q-p}\,\sum_{n=1}^{4}\int d{\bf K}\,a^{(p)\,*}_{n}({\bf
K})a^{(q)}_{n}({\bf K})=0$ (24)

In our one-dimensional case, if we start with an even solution
$\psi(x)=\psi_{S}(x)$, the corresponding odd solution $\psi_{A}(x)$ is found
equal to zero. And conversely, if we start with $\psi_{A}(x)$, the
corresponding even solution is equal to zero. Similarly, in our case with
$\pi/2$ rotation, one finds easily that if our starting $a_{n}$ has a
$a^{(p)}_{n}$ symmetry, the resulting solutions with different symmetry
$a^{(q)}_{n}$, with $q\neq p$, are identically equal to zero.

We note that the time-reversal counterpart of $a^{(1)}_{n}$, that is
$(-1)^{n-1}a^{(1)*}_{5-n}$, transforms like $a^{(4)}_{n}$, and similarly
$a^{(2)}_{n}$ and $a^{(3)}_{n}$ are linked in the same way. A solution and its
time-reversed function can not be identical because they have different
symmetries. Hence, the time reversability automatically implies that there is
a twofold degeneracy for the exciton eigenstates (this is Kramer’s degeneracy
theorem with respect to hole states). In order to reach the conclusion that
the degeneracy is actually fourfold, we have to search for a further symmetry
of the Hamiltonian.

For this purpose let us now consider the symmetry corresponding to a $\pi/2$
rotation $R_{x}$ around the $x$ axis. With respect to the orbital variables,
the situation is analogous to the one we had above with $R_{z}$: $a_{n}({\bf
K})$ becomes $a_{n}(R_{x}^{-1}({\bf K}))$, that is, $a_{n}(K_{x},K_{y},K_{z})$
becomes ${\tilde{a}}_{n}(K_{x},K_{y},K_{z})=a_{n}(K_{x},K_{z},-K_{y})$.
However, in order to obtain the transform of our four-dimensional wave
function, we again have naturally to take into account the fact that our basis
states are strongly modified under this rotation, that is, to calculate
$e^{-i\pi J_{x}/2\hbar}|n\rangle$. This is a standard calculation, that we
will not detail here. It can be done for example by finding the evolution of
the matrix elements $f_{n}(\varphi)=\langle n|e^{-i\varphi
J_{x}/\hbar}|m\rangle$ as a function of $\varphi$. For example for $m=1$ one
writes the first-order differential equations obtained by evaluating
$\hbar\partial f_{n}(\varphi)/\partial\varphi=-i\langle n|J_{x}\,e^{-i\varphi
J_{x}/\hbar}|1\rangle$ (for example $\partial
f_{1}(\varphi)/\partial\varphi=(-i\sqrt{3}/2)f_{2}(\varphi)$ from the explicit
expression of $J_{x}$) and one integrates these equations.

One obtains in this way the transform $|n\rangle_{x}\equiv e^{-i\pi
J_{x}/2\hbar}|n\rangle$ of our basis states by the $R_{x}$ rotation as

$\displaystyle|1\rangle_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{2}}\big{[}\,|1\rangle-i\,\sqrt{3}\,|2\rangle-\sqrt{3}\,|3\rangle+i\,|4\rangle\big{]}$
(25) $\displaystyle|2\rangle_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{2}}\big{[}-i\,\sqrt{3}\,|1\rangle-|2\rangle-i\,|3\rangle-\sqrt{3}\,|4\rangle\big{]}$
$\displaystyle|3\rangle_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{2}}\big{[}-\sqrt{3}\,|1\rangle-i\,|2\rangle-|3\rangle-i\,\sqrt{3}\,|4\rangle\big{]}$
$\displaystyle|4\rangle_{x}$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{2}}\big{[}i\,|1\rangle-\sqrt{3}\,|2\rangle-i\,\sqrt{3}\,|3\rangle+|4\rangle\big{]}$

Hence, the solutions $\sum_{n}a^{(p)}_{n}|n\rangle$ are transformed into
$\sum_{m,n}{\tilde{a}}_{n}^{(p)}|m\rangle\langle m|n\rangle_{x}$, where the
matrix elements $\langle m|n\rangle_{x}$ are given by Eq.(25). These
transformed solutions are naturally also solutions of the initial Schrödinger
equation (21). If the solutions $a^{(p)}_{n}$ are fourfold degenerate, this
implies that the transformed solutions are linear combinations of the
$a^{(p)}_{n}$, that is, there exists a set of coefficients $\lambda^{(q)}_{p}$
such that

$\displaystyle\sum_{n}{\tilde{a}}_{n}^{(q)}({\bf K})\langle
m|n\rangle_{x}=\sum_{p}\lambda^{(q)}_{p}a^{(p)}_{m}({\bf K})$ (26)

These simple linear relations between the solutions of Eq.(21) are fairly
remarkable, since the solutions of Eq.(21) are not expected to be simple.
Since we have no analytical solutions of Eq.(21), we have no way to check them
in the general case.

However, there is a limiting situation where we can solve analytically Eq.(21)
and check these relations. This is the limit of large wavevectors ${\bf K}$,
corresponding physically to short distances. In this case, since on the right-
hand side of Eq.(21) the $a_{m}({\bf K^{\prime}})$ have a limited range, we
can write in the integral $1/({\bf K-K^{\prime}})^{2}\simeq 1/{\bf K}^{2}$,
which provides an explicit solution for $a_{n}({\bf K})$ in terms of the
integrals $I_{n}=\int d{\bf K}\,a_{n}({\bf K})$, by merely inverting the
matrix $K^{2}\delta_{mn}+{\bar{h}}_{mn}({\bf K})$ as we have done in Eq.(19)
(${\bar{\epsilon}}$ is negligible in this large ${\bf K}$ limit). This gives
explicitly

$\displaystyle a_{m}({\bf
K})=\sum_{n=1}^{4}\frac{K^{2}\delta_{mn}-{\bar{h}}_{mn}({\bf
K})}{K^{4}-\lambda^{2}}\frac{I_{n}}{\pi^{2}K^{2}}$ (27)

If we take the particular case of the $a^{(p)}_{n}$ solutions, we notice that,
when we make in $I_{n}^{(p)}=\int d{\bf K}\,a_{n}^{(p)}({\bf K})$ the change
of variables on ${\bf K}$ corresponding to the $R_{z}$ rotation, the integral
is naturally unchanged. On the other hand, taking into account that
$a_{n}^{(p)}({\bf K})$ is transformed into ${\bar{a}}_{n}^{(p)}({\bf
K})=(-i)^{n-p}\,a^{(p)}_{n}({\bf K})$ in this change of variables, we find
that $I_{n}^{(p)}$ is multiplied by $(-i)^{n-p}$. Hence $I_{n}^{(p)}$ is zero
unless $n=p$, so that $I_{n}^{(p)}\equiv I^{(p)}\,\delta_{n,p}$. Actually we
can even choose to have $I^{(p)}=1$ by taking an appropriate norm for the
$a^{(p)}_{n}$ solutions. Hence, in this case Eq.(27) simplifies into

$\displaystyle a^{(p)}_{m}({\bf K})=\frac{K^{2}\delta_{mp}-{\bar{h}}_{mp}({\bf
K})}{\pi^{2}K^{2}(K^{4}-\lambda^{2})}$ (28)

In particular, one can check on these explicit expressions that the symmetry
properties ${\bar{a}}^{(p)}_{n}=(-i)^{n-p}\,a^{(p)}_{n}$ are indeed satisfied.

Integrating Eq.(26) over ${\bf K}$, and noting that $\int d{\bf
K}\,{\tilde{a}}_{n}^{(q)}({\bf K})=\int d{\bf K}\,a_{n}^{(q)}({\bf
K})=\delta_{n,q}$ by change of variables, we obtain the coefficients
$\lambda^{(q)}_{p}$ in Eq.(26) as $\lambda^{(n)}_{m}=\langle m|n\rangle_{x}$,
and in particular $\lambda^{(n)}_{n}=\langle n|n\rangle_{x}$. One can then
substitute in Eq.(26) the explicit solutions $a_{n}^{(q)}({\bf K})$ given by
Eq.(28), and check after a simple but tedious calculation that it is
satisfied. For example, taking the case $q=1$ and the component $m=1$, we have
to check from Eq.(25) that

$\displaystyle{\tilde{a}}_{1}^{(1)}({\bf
K})-i\sqrt{3}{\tilde{a}}_{2}^{(1)}({\bf K})-\sqrt{3}{\tilde{a}}_{3}^{(1)}({\bf
K})+i{\tilde{a}}_{4}^{(1)}({\bf K})=a_{1}^{(1)}({\bf
K})-i\sqrt{3}a_{1}^{(2)}({\bf K})-\sqrt{3}a_{1}^{(3)}({\bf
K})+ia_{1}^{(4)}({\bf K})$ (29)

where we have omitted the common denominator $2\sqrt{2}$ occuring in all the
terms from Eq.(25). Using the definition
${\tilde{a}}_{n}(K_{x},K_{y},K_{z})=a_{n}(K_{x},K_{z},-K_{y})$, and noting
that ${\tilde{a}}_{4}^{(1)}({\bf K})=a_{1}^{(4)}({\bf K})=0$, we are left with
checking that

$\displaystyle[K^{2}-{\bar{\gamma}}_{2}(K^{2}-3K_{y}^{2})]+6i{\bar{\gamma}}_{3}(K_{x}+iK_{z})K_{y}-[3{\bar{\gamma}}_{2}(K_{x}^{2}-K_{z}^{2})+6i{\bar{\gamma}}_{3}K_{x}K_{z}]$
$\displaystyle=$ (30)
$\displaystyle[K^{2}-{\bar{\gamma}}_{2}(K^{2}-3K_{z}^{2})]-6i{\bar{\gamma}}_{3}(K_{x}-iK_{y})K_{z}$
$\displaystyle-$
$\displaystyle[3{\bar{\gamma}}_{2}(K_{x}^{2}-K_{y}^{2})-6i{\bar{\gamma}}_{3}K_{x}K_{y}]$

where we have omitted the denominator $\pi^{2}K^{2}(K^{4}-\lambda^{2})$ common
to all the terms. We see that this equality Eq.(30) is indeed satisfied, but
in a nontrivial way.

It is now clear that the degeneracy of all the four solutions $a^{(p)}_{n}$ is
necessary in order to satisfy such an identity. Indeed, if on the contrary
$a^{(1)}_{n}$ (together with its time-reversed conjugate $a^{(4)}_{n}$) and
$a^{(2)}_{n}$ (together with $a^{(3)}_{n}$) had different energies, the
transformed of $a^{(1)}_{n}$ by the rotation $R_{x}$ would have to be a linear
combination of $a^{(1)}_{n}$ and $a^{(4)}_{n}$ only, since the two other
solutions had a different energy. This would lead to relations analogous to
Eq.(29), but with only two solutions ($a^{(1)}$ and $a^{(4)}$) on the right-
hand side instead of four as in Eq.(29). This is a much more stringent
algebraic requirement, that clearly can not be satisfied. Indeed, one can
check in the large ${\bf K}$ limit that the corresponding relation is
algebraically incompatible with the explicit solution. This provides an
explicit proof that the four solutions $a^{(p)}_{n}$ are indeed degenerate
remdegnum .

##  V Numerical results

We numerically find the ground state energy (and the corresponding wave
function) by applying an iterative procedure, to this Schrödinger equation
under the form written in Eq.(19). Let us call $A$ the linear operator that
acts on the four-dimensional wave function $a({\bf K})$ on the right-hand side
of Eq.(19). Solving Eq.(19) is equivalent to finding an eigenvector of $A$
with the eigenvalue $\Lambda=1$. In the case of a very large binding energy
${\bar{\epsilon}}$, the operator $A$ goes to zero and all its eigenvalues are
quite small. Hence, none of them can be equal to $1$, and there is no state
with very large binding energy, as expected. If we decrease
${\bar{\epsilon}}$, the largest positive eigenvalue $\Lambda_{\mathrm{m}ax}$
of $A$ will grow. When this eigenvalue of $A$ reaches $1$, we will have
obtained the largest possible value for ${\bar{\epsilon}}$ corresponding to an
eigenstate. In other words, we will have found the ground-state energy.

It is easy to obtain the largest eigenvalue of $A$ by applying iteratively $A$
to some convenient starting wave function $a({\bf K})$. A practical choice for
the starting wave function is to take the first component $a_{1}({\bf K})$
equal to the hydrogen atom ground-state wave function we have discussed above,
and the other components equal essentially to zero. Indeed, iterating $n$
times is equivalent to applying the operator $A^{n}$ to $a({\bf K})$. But for
large $n$ values, $A^{n}$ is dominated by its largest eigenvalue
$(\Lambda_{\mathrm{m}ax})^{n}$ and is essentially equivalent to a projection
on the corresponding eigenvector and multiplication by
$(\Lambda_{\mathrm{m}ax})^{n}$. This allows us to identify conveniently
$\Lambda_{\mathrm{m}ax}$ and the corresponding eigenvector. Actually this
procedure works only if the spectrum of $A$ does not have nasty features, such
as closely spaced largest and second largest eingenvalues, or large negative
eigenvalues. Fortunately, in our case, this procedure happens to work quite
nicely. We have found that, in practice, typically 20 iterations, or less,
gave already a satisfactory convergence for the precision we have considered.
It is also convenient, in order to find the ground-state energy, to start from
the hydrogen atom situation ${\bar{\gamma}}_{2}={\bar{\gamma}}_{3}=0$ where
the solution is known, and to crank up ${\bar{\gamma}}_{2}$ and
${\bar{\gamma}}_{3}$ progressively. In this way, the range where the ground-
state energy lies is fairly well known at each stage of the calculation.

In the practical task of performing the integral on the right-hand side of
Eq.(19), it is much better to get rid of the Coulomb interaction term and its
singular behaviour. This is done conveniently by performing the change of
variables ${\bf K}^{\prime}={\bf K}+{\bf Q}$. In this way the integral becomes

$\displaystyle\int d{\bf K}^{\prime}\frac{1}{({\bf K}-{\bf
K}^{\prime})^{2}}\,a({\bf
K}^{\prime})=\int_{0}^{\infty}dQ\int_{0}^{\pi}d\theta\sin\theta\int_{0}^{2\pi}d\varphi\;a(K_{x}+Q\sin\theta\cos\varphi,K_{y}+Q\sin\theta\sin\varphi,K_{z}+Q\cos\theta)$
(31)

where $\theta$ and $\varphi$ are the polar and azimuthal angles for ${\bf Q}$.
However, this change of variables makes it necessary, in order to conveniently
perform numerically the integrals, to evaluate $a(K_{x},K_{y},K_{z})$ for any
values of $K_{x},K_{y}$ and $K_{z}$. This can be done by sampling a number of
$K_{x},K_{y},K_{z}$ values and infer any $a(K_{x},K_{y},K_{z})$ value by
interpolation. In practice, it is rather $(1+K^{2})^{2}\,a(K_{x},K_{y},K_{z})$
that we have interpolated, since in the case of the hydrogen atom, this
function is just a constant equal to $1$. So, in our exciton case, we do not
expect it to have strong variations. This has allowed us to use a simple
three-dimensional linear interpolation. Naturally, we have also restricted the
range of variation of our variables by the change of variables $K_{x,y,z}=\tan
t_{x,y,z}$, with the regular discretization being on the $t_{x,y,z}$’s.

It is convenient to reduce the range of the ${\bf K}$ variables by making use
of symmetries, for example planar reflexions. If we consider first the change
$K_{z}\to-K_{z}$, it produces a change of sign for $D$ in Eq.(6), all the
other terms being unchanged. As a result, one checks easily that, if
$a_{n}(K_{x},K_{y},K_{z})$ is solution of Eq.(19), then
$(-1)^{n-1}a_{n}(K_{x},K_{y},-K_{z})$ is also solution. As above, we can build
a solution that is even under this transform (and is expected to correspond to
the ground state). It satisfies
$a_{n}(K_{x},K_{y},-K_{z})=(-1)^{n-1}a_{n}(K_{x},K_{y},K_{z})$. Note that the
time-reversed solution $a_{n}^{T}({\bf K})=(-1)^{n-1}a^{*}_{5-n}({\bf K})$
transforms as the odd combination in this transform
$a_{n}^{T}(K_{x},K_{y},-K_{z})=(-1)^{n}a_{n}^{T}(K_{x},K_{y},K_{z})$.

Similarly, in the change $K_{x}\to-K_{x}$, $F$ is changed into $F^{*}$ and $D$
into $-D^{*}$. One checks that, if $a_{n}(K_{x},K_{y},K_{z})$ is solution,
then $a_{5-n}(-K_{x},K_{y},K_{z})$ is solution. Hence, we can take a solution
with the symmetry $a_{n}(-K_{x},K_{y},K_{z})=a_{5-n}(K_{x},K_{y},K_{z})$.
Finally $K_{y}\to-K_{y}$ changes $F$ into $F^{*}$ and $D$ into $D^{*}$. So, if
$a_{n}(K_{x},K_{y},K_{z})$ is solution, then
$(-1)^{n-1}a_{5-n}(K_{x},-K_{y},K_{z})$ is solution, and we can take a
solution with the symmetry
$a_{n}(K_{x},-K_{y},K_{z})=(-1)^{n-1}a_{5-n}(K_{x},K_{y},K_{z})$. Note that,
with this choice of solution satisfying simultaneously these three symmetries,
if we perform a symmetry with respect to the origin $K_{x}\to-K_{x},K_{y}\to-
K_{y}$ and $K_{z}\to-K_{z}$, we find that our solution is invariant, as we
expect from the ground-state wave function. These symmetries allow us to
restrict the range of our variables to $K_{x}\geq 0,K_{y}\geq 0,K_{z}\geq 0$.

We want to explore the whole range of ${\bar{\gamma}}_{2}$ and
${\bar{\gamma}}_{3}$ values allowed by the condition Eq.(18). However, we note
that if we change the sign of both ${\bar{\gamma}}_{2}$ and
${\bar{\gamma}}_{3}$, from Eq.(6) $\Delta$, $D$ and $F$ change sign, that is,
$h_{h}$ defined by Eq.(5) is changed into $-h_{h}$. But $h_{h}$ and $-h_{h}$
must have the same eigenvectors, and moreover have also the same doubly
degenerate eigenvalues $\pm\lambda$ as seen in Eq.(7) and Eq.(8). In other
words, going from $h_{h}$ to $-h_{h}$ amounts to merely relabeling the
eigenvectors corresponding to the two positive eigenvalues as corresponding to
the two negative eigenvalues, and conversely. It is clear that such a
relabeling, which is merely a basis change, does not change the eigenvalue
spectrum of the whole Hamiltonian we are interested in. Accordingly, in
particular, the exciton ground-state energy is unchanged when changing the
sign of both ${\bar{\gamma}}_{2}$ and ${\bar{\gamma}}_{3}$. If we have a
solution for a given sign of the ${\bar{\gamma}}$’s, we are able to obtain a
solution having the same groud-state energy for the ${\bar{\gamma}}$’s with
opposite sign.

Let us now consider what happens if we change only the ${\bar{\gamma}}_{2}$
sign, ${\bar{\gamma}}_{3}$ being unchanged. From Eq.(6), one finds that
$\Delta$ is changed into $-\Delta$, and $F$ into $-F^{*}$, while $D$ is
unchanged. Moreover, as we have seen above, in the change of variables
$K_{x}\to K_{y}$ and $K_{y}\to-K_{x}$, $D({\bf K})$ becomes $i\,D({\bf K})$,
while $F({\bf K})$ becomes $-F({\bf K})$. But we notice that a further simple
change of functions $a_{2}\to-ia_{2},a_{4}\to-ia_{4}$ is equivalent to
changing $D$ into $iD$. Combining these three changes, we have changed
$\Delta$ into $-\Delta$, $F$ into $F^{*}$ and $D$ into $-D$. But the resulting
equations for Eq.(17), or Eq.(19), are identical to the original ones provided
we make the additional change of functions $a_{1}\leftrightarrow
a_{3},a_{2}\leftrightarrow a_{4}$. Hence, we conclude that
${\bar{\gamma}}_{2}\to-{\bar{\gamma}}_{2}$ leaves the energy spectrum
invariant. Combining now the two changes of signs for the $\gamma$’s, we
obtain that the ground-state energy depends only on $|{\bar{\gamma}}_{2}|$ and
$|{\bar{\gamma}}_{3}|$, so we can restrict our study to
$0\leq{\bar{\gamma}}_{2}<1/2$, $0\leq{\bar{\gamma}}_{3}<1/2$.

Let us finally examine the limiting cases ${\bar{\gamma}}_{2}\to
1/2,{\bar{\gamma}}_{3}\to 1/2$. As we have seen below Eq.(10), when either
${\bar{\gamma}}_{2}=1/2$ or ${\bar{\gamma}}_{3}=1/2$, there are some ${\bf K}$
directions for which $K^{2}-\lambda^{2}=0$. As a result, if in Eq.(19)
${\bar{\epsilon}}$ is fixed, for large values of $K$, the denominator in the
first factor of the right-hand side will go to zero. This implies that $a({\bf
K})$ will not behave for all ${\bf K}$ directions as $1/K^{2}$, as we have
found in Eq.(20). This leads to a divergent behaviour in the integral of
Eq.(19). The natural way to escape this problem is to let ${\bar{\epsilon}}$
be very large, so this divergent behaviour of $a({\bf K})$ is pushed to ever
higher values of $K$, and one avoids an actual singularity in the limits
${\bar{\gamma}}_{2}\to 1/2,{\bar{\gamma}}_{3}\to 1/2$.

Unfortunately, an explicit handling of this behaviour is not so easy in
general because it arises from the behaviour of $a({\bf K})$ around some
specific directions of ${\bf K}$. Nevertheless when we let both
${\bar{\gamma}}_{2}\to 1/2,{\bar{\gamma}}_{3}\to 1/2$ at the same time, it is
possible to extract the divergent behaviour of ${\bar{\epsilon}}$ by a simple
rescaling. Let us set
$2\,{\bar{\gamma}}_{2}=1-\eta_{2},2\,{\bar{\gamma}}_{3}=1-\eta_{3}$, and
consider $\eta_{2}\to 0,\eta_{3}\to 0$ with $r=\eta_{3}/\eta_{2}$ fixed to a
finite value. We perform the rescaling ${\bf K}={\bar{\bf K}}/\eta_{2}$, while
at the same time setting ${\bar{\epsilon}}=\alpha/\eta_{2}$. This leads, for
the denominator on the right-hand side of Eq.(19), to
$(K^{2}+{\bar{\epsilon}})^{2}-\lambda^{2}=(2/\eta_{2}^{3})[{\bar{K}}^{4}+\alpha{\bar{K}}^{2}+3(r-1)({\bar{K}}_{x}^{2}\,{\bar{K}}_{y}^{2}+{\bar{K}}_{y}^{2}\,{\bar{K}}_{z}^{2}+{\bar{K}}_{z}^{2}\,{\bar{K}}_{x}^{2})]$.
In this rescaling, the numerator provides a factor $1/\eta_{2}^{2}$, while
another factor $1/\eta_{2}$ comes from the rescaling of the variable ${\bf
K}^{\prime}$ into ${\bar{\bf K}}^{\prime}$ in the integral. Hence, the scaling
factor $\eta_{2}$ disappears completely and we are left with the following
equation, free of singularity

$\displaystyle b({\bar{\bf
K}})=\frac{{\bar{K}}^{2}\mathbb{1}-{\bar{h}}_{h}({\bar{\bf
K}})}{2[\alpha{\bar{K}}^{2}+{\bar{K}}^{4}+3(r-1)({\bar{K}}_{x}^{2}\,{\bar{K}}_{y}^{2}+{\bar{K}}_{y}^{2}\,{\bar{K}}_{z}^{2}+{\bar{K}}_{z}^{2}\,{\bar{K}}_{x}^{2})]}\frac{1}{\pi^{2}}\int
d{\bar{\bf K}}^{\prime}\frac{1}{({\bar{\bf K}}-{\bar{\bf
K}}^{\prime})^{2}}\,b({\bar{\bf K}}^{\prime})$ (32)

where $b({\bar{\bf K}})=a({\bar{\bf K}}/\eta_{2})$. We have solved numerically
this equation for the scaling factor $\alpha$ in the energy, in the case
$r=1$, and we have found $\alpha=0.43$. This is in fair agreement with the
corresponding reduced energy ${\bar{\epsilon}}=4.77$ we have found by a direct
numerical solution of Eq.(19) for $\eta_{2}=\eta_{3}=0.1$.

Following the procedure indicated above, we have calculated the reduced
exciton ground-state energy ${\bar{\epsilon}}$ from Eq.(19) as a function of
${\bar{\gamma}}_{2}$ and ${\bar{\gamma}}_{3}$, on a grid with spacing $0.05$
for these two parameters. Rather than presenting these extensive results as a
table, we have used them to evaluate ${\bar{\epsilon}}$ for any value of
${\bar{\gamma}}_{2}$ and ${\bar{\gamma}}_{3}$ by a two-dimensional spline
interpolation, and then obtain the level lines for various values of
${\bar{\epsilon}}$ in the ${\bar{\gamma}}_{2}-{\bar{\gamma}}_{3}$ plane. They
are displayed in Fig.1. We have chosen not to smooth out the small
irregularities which are apparent along these level lines, since they result
from the imprecision of our calculations and accordingly give a fairly direct
information on it. In particular, the results for ${\bar{\epsilon}}$ close to
$1$ are very sensitive to this imprecision, and the resulting line for
${\bar{\epsilon}}=1.05$ would be pretty bad. Hence, to draw this line, we have
rather relied on a quadratic interpolation from the ${\bar{\gamma}}_{2}=0.1$
and ${\bar{\gamma}}_{3}=0.1$ results, rather than using our
${\bar{\gamma}}_{2}=0.05$ and ${\bar{\gamma}}_{3}=0.05$ results.

Figure 1: Reduced ground-state exciton binding energy ${\bar{\epsilon}}$ in
the ${\bar{\gamma}}_{2}-{\bar{\gamma}}_{3}$ plane. The level lines for
${\bar{\epsilon}}=1.05,1.1,1.2,1.3,1.5,1.7$ and $2.$ are represented, as
indicated in the figure.

The numerical results are fairly regular, with the represented level lines
being not so far from quarter circles, provided we do not go too close to the
boundaries ${\bar{\gamma}}_{2}=0.5$ and ${\bar{\gamma}}_{3}=0.5$. This means
that, roughly speaking, ${\bar{\epsilon}}$ depends only on
${\bar{\gamma}}_{2}^{2}+{\bar{\gamma}}_{3}^{2}\equiv 2{\bar{\gamma}}^{2}$.
Correspondingly, we find that ${\bar{\epsilon}}-1\simeq
2.5(1/\sqrt{1-4{\bar{\gamma}}^{2}}-1)$ provides an approximate representation
of the numerical results. In particular, it is quadratic for small
${\bar{\gamma}}$ and it diverges for ${\bar{\gamma}}\to 1/2$. It gives
${\bar{\epsilon}}\simeq 4.23$ for ${\bar{\gamma}}=0.45$, not so different from
the result from our numerical ${\bar{\epsilon}}=4.77$ mentioned above.

Regarding the size of the results, we note that taking our reduced ground-
state exciton binding energy as unity ${\bar{\epsilon}}=1$ corresponds to
making use, in the standard expression
$E=-\mu_{X}(e^{2}/4\pi\epsilon_{sc})^{2}/(2\hbar^{2})$, of a reduced mass
$\mu_{X}$ given by $1/\mu_{X}=1/m_{e}+1/\mu_{h}$, where from Eq.(10) and
Eq.(13) the reduced hole mass $\mu_{h}$ itself is taken equal to the half-sum
of the reduced heavy and light hole masses $1/\mu_{h}=(1/m_{H}+1/m_{L})/2$.
This is a frequently done approximation, for lack of better knowledge. Fig.1
shows that in principle we may have a fairly important departure from this
simple approximation.

However, we have considered in practice a few semiconductors where reliable
data exist for the valence band as well as the conduction band parameters
cardona ; kara ; vurga ; vurgap . One can still see sizeable theoretical and
experimental uncertainties in the knowledge of these band parameters kara ;
vurga . Remarkably, the strong variations among compounds in the basic
Luttinger parameters and conduction electron mass, provide a much reduced
variation in our reduced Luttinger parameters ${\bar{\gamma}}_{2}$ and
${\bar{\gamma}}_{3}$. We note that a light electronic mass increases
${\bar{\gamma}}_{1}$ from Eq.(13), which leads to a reduction of the size of
${\bar{\gamma}}_{2}$ and ${\bar{\gamma}}_{3}$. We find that
${\bar{\gamma}}_{2}$ is mostly in the range $0.1-0.15$, while
${\bar{\gamma}}_{3}$ takes slightly higher values in the range $0.1-0.2$. From
our results, the warping leads to an increase of the exciton binding energy by
about $10\%$. This should be taken into account as soon as one looks for some
precision in the exciton binding energy. The largest correction we have found
for a compound with specifically known parameters is in the case of BN vurga ,
where ${\bar{\gamma}}_{2}=0.12$ and ${\bar{\gamma}}_{3}=0.2$ leads to a $15\%$
increase from the warping in the exciton binding energy. However, it should be
kept in mind that this binding energy is quite sensitive to the parameters
values in this range. For example having
${\bar{\gamma}}_{2}={\bar{\gamma}}_{3}=0.2$ would lead to a $20\%$ increase in
the binding energy, and ${\bar{\gamma}}_{2}={\bar{\gamma}}_{3}=0.25$ would
give a $35\%$ increase. Hence, keeping in mind the uncertainties on the
valence band parameters, the effect of warping on the exciton binding energy
may happen to be quite important.

##  VI Conclusion

In this paper, we have addressed the effect on the exciton ground-state energy
of the warping of the valence band near its edge. We have specifically
considered the case of bulk exciton for semiconductors with zinc-blende
crystal structure. Assuming as usual that the split-off subband due to spin-
orbit coupling is energetically far enough to play a negligible role, we have
considered the four-dimensional problem for the hole kinetic energy. We have
shown that, due to the cubic symmetry, the exciton ground state has a fourfold
degeneracy.

We have introduced, by a simple scaling, reduced Luttinger parameters
${\bar{\gamma}}_{2}$ and ${\bar{\gamma}}_{3}$, whose absolute value are less
than $1/2$. We have studied systematically the exciton ground-state energy as
a function of these reduced Luttinger parameters, by solving numerically the
integral equation corresponding to the Schrödinger equation in Fourier
transform. We find that, compared to its standard expression in the absence of
warping (corresponding to the case ${\bar{\gamma}}_{2}={\bar{\gamma}}_{3}=0$),
the exciton binding energy can in principle increase without limitations. For
moderate increases, we have provided an approximate analytical expression for
the increase in terms of the Luttinger parameters. Going through the values of
band structure parameters, and in particular Luttinger parameters, for
tabulated semiconductors, we find a typical $10\%-15\%$ increase in the
exciton binding energy. However, this energy is fairly sensitive to the values
of the Luttinger parameters, and a slight increase in their values due to
experimental or/and theoretical uncertainties could lead to a markedly larger
increase in the exciton binding energy.

Acknowledgements: The authors are grateful to Monique Combescot for
interesting them in this warping problem and for discussions.

## References

  * (1) P.Y. Yu and M. Cardona, _Fundamentals of Semiconductors_ (Springer, Berlin, 1996).
  * (2) I. Vurgaftman, M.P. Lumb and J.R. Meyer, _Bands and Photons in III-V Semiconductor Quantum Structures_ (Oxford University Press, 2021).
  * (3) J. M. Luttinger, Phys. Rev. 102, 1030 (1956).
  * (4) J. M. Luttinger and W. Kohn, Phys. Rev. 97, 368 (1955).
  * (5) In terms of those of Ref.lutkoh , our coefficients are given by $\hbar^{2}\gamma_{1}/m_{0}=2(A+2B)/3$, $\hbar^{2}\gamma_{2}/m_{0}=(A-B)/3$ and $\hbar^{2}\gamma_{3}/m_{0}=C/3$, and we have an additional $\pi/2$ rotation around the $z$ axis. These $\gamma_{i}$ are just the coefficients introduced by Luttinger lutt .
  * (6) G. Dresselhaus, A. F. Kip and C. Kittel, Phys. Rev. 98, 368 (1955).
  * (7) The warping of the conduction band is expected to be very weak in the relevant ${\bf k}$ range for the exciton; see _e_.g. U. Rössler, Solid State Comm. 49, 943 (1984).
  * (8) Shiue-Yuan Shiau and Monique Combescot (unpublished)
  * (9) In the simple case ${\bar{\gamma}}_{3}=0$ where this $4\times 4$ problem reduces to two $2\times 2$ problems, we have checked numerically this degeneracy (within our numerical precision).
  * (10) I. Vurgaftman, J.R. Meyer and L.R. Ram-Mohan, Journal of Applied Physics 89, 5815 (2001).
  * (11) S.Zh. Karazhanov and L.C. Lew Yan Voon, Semiconductors 39, 161 (2005).


Engineering Laboratory, University of Cambridge, Trumpington Street,
Cambridge, CB2 1PZ, United Kingdom

# Preconditioners for the geometry optimisation and saddle point search of
molecular systems

Letif Mones<EMAIL_ADDRESS>Mathematics Institute, University of Warwick,
Coventry CV4 7AL, United Kingdom  Gábor Csányi  Engineering Laboratory,
University of Cambridge, Trumpington Street, Cambridge, CB2 1PZ, United
Kingdom  Christoph Ortner  Mathematics Institute, University of Warwick,
Coventry CV4 7AL, United Kingdom

##  1 Abstract

A class of preconditioners is introduced to enhance geometry optimisation and
transition state search of molecular systems. We start from the Hessian of
molecular mechanical terms, decompose it and retain only its positive definite
part to construct a sparse preconditioner matrix. The construction requires
only the computation of the gradient of the corresponding molecular mechanical
terms that are already available in popular force field software packages. For
molecular crystals, the preconditioner can be combined straightforwardly with
the exponential preconditioner recently introduced for periodic systems. The
efficiency is demonstrated on several systems using empirical, semiempirical
and ab initio potential energy surfaces.

##  2 Introduction

Geometry optimisation and transition state search are fundamental procedures
to identify important stationary points of molecules, molecular crystals and
material systems in computational chemistry. Since the evaluation of
chemically accurate ab initio potential energies and gradients are
computationally demanding, several techniques have been developed over the
last three decades to enhance the convergence of optimisation methods.

Among the most widely used are quasi-Newton methods, in particular BFGS or its
limited memory version 1, which start with a the identity as its guess for the
Hessian and update it at each iteration based on the gradient information
collected from previous steps. Initializing with Hessian guess that include
more molecule specific geometrical information can improve the convergence.
For instance, just introducing some connectivity information about a molecule
can lead to surprisingly good results2, 3.

More sophisticated methods estimate the initial Hessian guess from a surrogate
potential (e.g. force field or semiempirical potential) whose second
derivative can be obtained at low computational cost. The update of the
approximate Hessian can then be achieved either by using quasi-Newton methods
or other techniques such as DIIS 4, 5. Such strategies significantly improve
the speed of convergence in either Cartesian or internal coordinates.6

If the model Hessian is cheap to calculate (e.g., as obtained from a surrogate
model) and provides a reasonable approximation of the quantum Hessian, it may
be advantageous to recompute it at every optimisation step. Such a scheme was
introduced by Lindh et al. 7, where a model potential is constructed
consisting of quadratic terms for all distances, angles and dihedrals in the
molecule. At each geometry optimisation step the force field is constructed
such that the current conformation is its local minimum, and its “Hessian” is
computed, neglecting the dependence of the force field parameters on the
geometry. With this construction, the model “Hessian” is therefore not the
Hessian of any potential. Nevertheless, this approach yields excellent
performance, which led to its wide implementation 8.

The model “Hessian” of Lindh 7 can be considered as a preconditioner or metric
that effects a transformation to a new coordinate system where the
optimisation problem is better conditioned, hence algorithms converge more
rapidly and tend to be more robust. Geometrically, the shape of the energy
landscape becomes more isotropic. In general it is desirable that both the
construction and inversion of the preconditioner matrix are inexpensive (at
least compared to the computation of the energy and gradient). This can be
achieved by building a sparse preconditioner from simple analytical functions.
Since the preconditioner defines a metric in configuration space it needs to
be positive definite. This requirement is automatically fulfilled in the Lindh
approach7 by the use of quadratic molecular mechanical terms with equilibrium
points corresponding to the actual geometry at each step.

We recently introduced a general and effective preconditioner for geometry
optimisation and saddle search for material systems 3. This preconditioner is
determined by the local connectivity structure of atoms making both the
construction and inversion computationally inexpensive. Especially for larger
systems we observe an order of magnitude or larger reduction of the number of
optimisation steps for the preconditioned LBFGS method compared to the same
without preconditioning. Here we expand our previous work to molecules and
molecular crystals by combining it with a force field based preconditioner
inspired by the approach of Lindh et al. 7.

##  3 Methods

###  3.1 Enhancing geometry optimisation by using preconditioners

We briefly review the methodology for preconditioning geometry optimisation
and the dimer saddle search method for material systems3. For a system with
$N$ particles let $x_{k}\in\mathbb{R}^{3N}$ denote the configuration at the
$k$th iterate of an optimisation algorithm. The corresponding energy, gradient
and preconditioner are denoted by $f_{k}=f(x_{k})$, $g_{k}=\nabla f(x_{k})$
and $P_{k}\in\mathbb{R}^{3N\times 3N}\approx\nabla^{2}f(x_{k})$, respectively.
A preconditioned steepest-descent step is then given by

$x_{k+1}=x_{k}-\alpha_{k}P_{k}^{-1}g_{k}$ (1)

where $\alpha_{k}$ is the step size obtained from some line search procedure
at the $k$th iteration. If $P_{k}=I$, then (1) becomes the standard steepest
descent scheme and if $P_{k}=\nabla^{2}f(x_{k})$, then it becomes a
quadratically convergent Newton scheme. In general, different choices of
$P_{k}$ may “interpolate” between these extremes.

From an alternative point of view preconditioning can be thought of as a
coordinate transformation, where a new set of coordinates is defined as
$y_{k}:=P^{1/2}x_{k}$. The advantage of this framework is that once an
appropriate preconditioner matrix is available then any optimisation algorithm
can be modified by applying the original algorithm on the transformed
coordinates. To obtain the final form of the modified algorithm we need to
transform the variables back to the original coordinate system. As a simple
example, applying the coordinate transformation on the gradient descent
equation immediately leads to the equation of quasi Newton schemes (eq. 1):

$\displaystyle\begin{split}y_{k}&=y_{k-1}-\alpha_{k}\nabla_{y}F(y_{k})\\\
P^{1/2}_{k}x_{k}&=P^{1/2}_{k}x_{k-1}-\alpha_{k}P^{-1/2}_{k}\nabla_{x}f(x_{k})\\\
x_{k}&=x_{k-1}-\alpha_{k}P^{-1}_{k}\nabla_{x}f(x_{k})\\\ \end{split}$ (2)

where we used the fact that
$\nabla_{y}F(y_{k})=P^{-1/2}_{k}\nabla_{x}f(x_{k})$. Preconditioning popular
optimisation methods like LBFGS, conjugate gradients, FIRE9 etc. is similarly
possible.

A simple preconditioner that is effective for a wide range of materials
systems is based on the following $N\times N$ matrix 3

$L_{ij}=\begin{cases}-\mu\exp\left(-A\left(\frac{r_{ij}}{r_{\mathrm{nn}}}-1\right)\right),&i\neq
j\text{~{}and~{}}|r_{ij}|<r_{\mathrm{cut}}\\\ 0,&i\neq
j\text{~{}and~{}}|r_{ij}|\geq r_{\mathrm{cut}}\\\ -\sum_{i\neq
j}L_{ij},&i=j\end{cases}$ (3)

where $i$ and $j$ denote atomic indices and $\mu$, $A$, $r_{\mathrm{cut}}$ and
$r_{\mathrm{nn}}$ are parameters that can be user-specified or estimated
numerically. We note that (3) is a generalisation of the Laplacian matrix used
to represent undirected graphs. Given a specific connectivity defined by
$r_{\mathrm{cut}}$ and setting $A=0$ and $\mu=1$, $L_{ij}$ reduces exactly to
the Laplacian matrix. The actual $3N\times 3N$ preconditioner is simply
obtained from the corresponding $L_{ij}$ element in an isotropic manner:

$[P_{\mathrm{Exp}}]_{i+3(k-1),j+3(l-1)}=\begin{cases}L_{ij},&k=l\\\ 0,&k\neq
l\end{cases}$ (4)

where the $k$ and $l$ indices denote Cartesian components. Despite its
simplicity in capturing only geometric connectivity but no specific material
information, $P_{\rm Exp}$ provides a good model for the local curvature of
the potential energy landscape, which effectively controls ill-conditioning in
large systems. The application of $P_{\rm Exp}$ resulted in a significant
reduction of the number of optimisation steps required for several material
systems compared to the unpreconditioned LBFGS3.

###  3.2 FF-based preconditioners

Preliminary tests showed that for molecular systems such as molecules in gas
phase or molecular crystals using $P_{\rm Exp}$ still leads to a speed-up, but
a much more modest one than for material systems. The explanation for this is
that molecular systems contain a wide range of different interactions (e.g.,
pair, angle, dihedral, electrostatic, dispersive) of vastly varying stiffness
which in addition are more loosely coupled, and this creates a second source
of ill-conditioning distinct from ill-conditioning due to large system size.
Inspired by the use of internal coordinates in molecular optimisation
techniques10 and the model Hessian of Lindh et al. 7 we therefore propose a
generalisation of $P_{\rm Exp}$ that is effective also for molecular systems.

The construction of our FF-based preconditioner begins with a surrogate
potential energy function, given by a sum over internal coordinates each
describing a short-range bond in the system (distance, angle, or dihedral),

$V_{\mathrm{FF}}=\sum\limits_{\alpha}V_{\alpha}(\xi_{\alpha}(\mathbf{r})).$
(5)

The individual potential energy terms are in general simple functions of the
internal coordinates. Some examples of most typical forms are the quadratic,
Morse or torsional potentials, respectively given by

$\displaystyle V_{\mathrm{Quadratic}}(q)$
$\displaystyle=\frac{1}{2}k(q-q_{0})^{2},$ (6) $\displaystyle
V_{\mathrm{Morse}}(d)$
$\displaystyle=D_{0}(1-\exp(-\alpha(d-d_{0})))^{2},\quad\text{and}$ (7)
$\displaystyle V_{\mathrm{Torsion}}(\phi)$
$\displaystyle=\frac{1}{2}k_{\phi}(1+\cos(n\phi-\phi_{0}))$ (8)

where the corresponding parameters can be taken from standard force field
libraries.

Due to its simple functional form $H_{\rm FF}:=\nabla^{2}V_{\mathrm{FF}}$ is
cheap to compute. If only short-range bonds are taken into account then it is
also sparse, hence it is cheap to store and invert. Moreover, we expect that
$V_{\rm FF}$ gives a good qualitative approximation to the quantum potential
energy landscape, hence $H_{\rm FF}$ is a good qualitative approximation to
the quantum Hessian $\nabla^{2}f$. Therefore, $H_{\rm FF}$ satisfies all the
conditions required for a preconditioner except that it will in general be
indefinite. A conceptually straightforward but computationally expensive
approach to overcome this limitation is to enforce positivity by replacing all
eigenvalues of $H_{\rm FF}$ with their absolute values. Instead, we propose to
analytically modify the local Hessian contributions to ensure overall
positivity, resulting in a further reduction in computational cost.

The Hessian contribution from $V_{\alpha}$ is given by

$\mathbf{H}_{\alpha}=\frac{\partial^{2}V_{\alpha}}{\partial\mathbf{r}^{2}}=\underbrace{\frac{\partial\xi_{\alpha}}{\partial\mathbf{r}}\otimes\frac{\partial\xi_{\alpha}}{\partial\mathbf{r}}\frac{\partial^{2}V_{\alpha}}{\partial\xi_{\alpha}^{2}}}_{\mathbf{H}_{\alpha}^{(1)}}+\underbrace{\frac{\partial^{2}\xi_{\alpha}}{\partial\mathbf{r}^{2}}\frac{\partial
V_{\alpha}}{\partial\xi_{\alpha}}}_{\mathbf{H}_{\alpha}^{(2)}},$ (9)

where we have decomposed $H_{\alpha}$ into two terms,
$\mathbf{H}_{\alpha}^{(1)}$ and $\mathbf{H}_{\alpha}^{(2)}$. If $V_{\alpha}$
is quadratic then
$\frac{\partial^{2}V_{\mathrm{Quadratic}}}{\partial\xi^{2}}=k>0$, hence
$\textbf{H}_{\alpha}^{(1)}$ is positive semi-definite, while the sign of
$\textbf{H}_{\alpha}^{(2)}$ is ambiguous. Note however, that if the system is
at equilibrium of $V_{\alpha}$ with respect to $\xi_{\alpha}$, i.e., if
$\frac{\partial V_{\alpha}}{\partial\xi_{\alpha}}=0$, then
$\textbf{H}_{\alpha}^{(2)}=0$. This in fact is the case in the Lindh approach
7. Instead of adjusting $V_{\alpha}$ at every step such that the geometry
corresponds to its equilibrium, here we simply drop
$\textbf{H}_{\alpha}^{(2)}$ and only use $\mathbf{H}_{\alpha}^{(1)}$ to
construct the preconditioner, thus ensuring that it always stays positive
definite.

For non-quadratic contributions we expect that
$\frac{\partial^{2}V_{\alpha}}{\partial\xi_{\alpha}^{2}}>0$ for most but not
all bonds, hence we enforce positivity by replacing it with its absolute
value. This leads to the following general preconditioner for molecular
systems:

$P_{\mathrm{FF}}=\sum\limits_{\alpha}\tilde{\mathbf{H}}_{\alpha}^{(1)}=\sum\limits_{\alpha}\frac{\partial\xi_{\alpha}}{\partial\mathbf{r}}\otimes\frac{\partial\xi_{\alpha}}{\partial\mathbf{r}}\left|\frac{\partial^{2}V_{\alpha}}{\partial\xi_{\alpha}^{2}}\right|$
(10)

It is worth noting that $\frac{\partial\xi_{\alpha}}{\partial\mathbf{r}}$ is
already computed by molecular mechanics force field based MD programs since it
is required for the assembly of $\nabla V_{\alpha}$. Thus, the only new
quantity that must be computed is
$\frac{\partial^{2}V_{\alpha}}{\partial\xi_{\alpha}^{2}}$, which represents a
negligible additional computational cost.

We note that the final functional form (10) of our preconditioner is very
similar to that of Lindh et al. 7, however we arrived at it from a
fundamentally different perspective, which has several advantages. Lindh et
al.’s method was introduced for quadratic terms only, hence the force
constants have to be recomputed after each optimisation step (as the
equilibrium bond lengths, angles and dihedrals are set to the actual ones to
obtain positive semidefinite matrix). Our method can be considered as a
generalisation of their approach, allowing arbitrary functional forms of
internal coordinate dependent terms. In particular this means that the FF
parameters need not be adjusted to achieve a positive semidefinite matrix, and
incorporating different parameter sets is straightforward. Moreover, as we
will see in § 3.3 and in § 4.2 our perspective makes it easy to extend the
preconditioner construction to new situations.

###  3.3 Combining FF and Exp preconditioners

For molecular crystal, intermolecular interactions also play an important
role. This consideration led us to combine the molecular mechanics based FF
preconditioner (describing bonded interactions) with the Exp preconditioner
(4) (tuned to describe non-bonded interactions):

$P_{\mathrm{Exp+FF}}=P_{\mathrm{Exp}}+P_{\mathrm{FF}}$ (11)

$P_{\rm FF}$ is fully specified from the chosen force field $V_{\rm FF}$. To
specify $P_{\rm Exp}$ we first manually choose the parameters $A$ and
$r_{\mathrm{cut}}$ in (4) to account for the interaction between molecules3.
The remaining parameters are computed in a similar automatic manner as
described in Ref. 3, and we keep only those matrix elements of $P_{\rm Exp}$
for which the corresponding matrix element in $P_{\rm FF}$ is zero. We note
that correct scaling between $P_{\rm FF}$ and $P_{\rm Exp}$ is implicitly
ensured via the $\mu$ parameter in Eq 3.

###  3.4 Implementation details

We tested the FF and FF+Exp preconditioners on a range of optimisation and
saddle search tasks. For geometry optimisations the form of the preconditioned
LBFGS method was identical to the one we describe in ref3: At each iterate the
search direction is given by

$\displaystyle\begin{split}&{\bf input}~{}q=\nabla f(x_{k})\\\
&s_{k}=x_{k}-x_{k-1}\\\ &y_{k}=\nabla f(x_{k})-\nabla f(x_{k-1})\\\
&\rho_{k}=1/y_{k}^{T}s_{k}\\\ &\mathrm{for}\>i=k,\ldots,k-m\\\
&\qquad\alpha_{i}=\rho_{i}s_{i}^{T}q\\\ &\qquad q=q-\alpha_{i}y_{i}\\\
&\framebox{$z=P_{k}^{-1}q$}\\\ &\mathrm{for}\>i=k-m,\ldots,k\\\
&\qquad\beta_{i}=\rho_{i}y_{i}^{T}z\\\ &\qquad
z=z+(\alpha_{i}-\beta_{i})s_{i}\\\ &{\bf output}~{}p_{k}=z\\\ \end{split}$
(12)

with initial search direction $z=P_{0}^{-1}\nabla f(x_{0})$. The step length
selection is obtained by a backtracking line-search enforcing only the Armijo
condition3.

For saddle search tasks we slightly modified the superlinearly converging
dimer method 11. Dimer methods use two copies of the system with coordinates
$x_{1}$ and $x_{2}$ and a fixed separation length $l=|x_{1}-x_{k}|$ between
them. The algorithm is usually split into two alternating steps 12: (1) in the
rotation step we fix the midpoint and rotate the endpoints to approximately
align them with the lowest (negative) eigenmode of the Hessian; (2) in the
translation step we shift the dimer to maximise the energy along the dimer
direction while minimising energy in all directions perpendicular to it.

In principle, both the rotation and translation steps can be preconditioned,
however, we found that in many systems preconditioning the rotation step
results in a smaller spectral gap and hence slower convergence. Therefore, we
chose to precondition only the translation step. Our implementation employs
the conjugate gradient method using the Polak-Ribière formula:

$\displaystyle\begin{split}&{\bf input}~{}q_{k}=-(I-2v_{k}\otimes v_{k})\nabla
f(x_{k})\\\
&\framebox{$\beta=\frac{q_{k}^{T}P_{k}^{-1}(q_{k}-q_{k-1})}{q_{k-1}^{T}P_{k}^{-1}q_{k-1}}$}\\\
&s_{k}=q_{k}+\beta s_{k-1}\\\ &{\bf
output}~{}\framebox{$p_{k}=\frac{P^{-1}s_{k}}{s_{k}^{T}P^{-1}s_{k}}$}\\\
\end{split}$ (13)

and using the initial iterate $s_{0}=-(I-2v_{0}\otimes v_{0})\nabla f(x_{0})$.
For computing the step length we used the trust region radius approach
suggested by Kästner and Sherwood11, with acceptance criterion based on the
projection of the gradient of the actual step.

For molecules in gas phase $V_{\rm FF}$ is invariant under rotations and
translations, hence $P_{\rm FF}$ will be at least six fold degenerate with
zero eigenvalues for any configuration of the molecule corresponding to the
three translational and three rotational degrees of freedom. While these
degrees of freedom could in principle be fixed we found that a straightforward
solution is to simply regularise the preconditioner by replacing it with
$P\leadsto P+cI$. We found that good generic values for $c$ are 0.1 and 1.0 eV
/ Å2 for, respectively, geometry optimisations and saddle search.

###  3.5 Model systems and potentials

####  3.5.1 Organic molecules in gas phase

Three potential energy surfaces were investigated for geometry optimisations:
semiempirical PM613, DFT14, 15 and MP216. For the DFT potential we used the
the BLYP exchange-correlation functional17, 18 with the DZVP-MOLOPT basis
set19 and plane wave cutoff of 480 Ry within the Gaussian and plane waves
method (GPW) approach20 and Goedecker-Teter-Hutter (GTH) pseudopotentials21.
For the MP2 calculations the 6-31G** basis set was applied.

In the case of geometry optimisations on the PM6 surface we compared three
different force fields from which we constructed the preconditioners: the
force field of Lindh et al. (LFF)7, the universal force field (UFF)22 and the
generalised Amber force field (GAFF)23.

For transition state search we selected 7 examples from the benchmark of Baker
and Chan24 and three additional systems. The computations were performed on
the semiempirical PM6 surface13.

####  3.5.2 Molecular crystals

We compared four different optimisation schemes (unpreconditioned, only FF-
based, only Exp-based and Exp+FF-based preconditioners) on five organic
molecular crystals (systems XVIII to XXII) taken from the Organic Crystal
Structure Prediction competition of the Cambridge Crystallographic Data Centre
25, 26. We used a DFT potential energy surface with the PBE exchange-
correlation functional27 with a plane wave basis set using a cutoff energy of
800 eV and ultrasoft pseudopotentials28.

####  3.5.3 Material systems

We also tested how the force field based preconditioner works on two material
systems compared to the exponential preconditioner. We examined the
unpreconditioned and preconditioned geometry optimisation for bulk silicon and
vacancy with varying system size. The potential energy surface was the
screened Tersoff potential29, 30 and we used the universal force field (UFF)22
for building the FF-based preconditioner matrix. Next we considered bulk
tungsten and a single interstitial site in bulk tungsten, also using different
system sizes. The potential energy surface was a machine learning based
Gaussian Approximation Potential (GAP) reproducing the quality of DFT (with
PBE functional) 31. The preconditioner in this case was based on a simple
Embedded Atom Method (EAM) potential 32, 33.

####  3.5.4 Software

For the semiempirical PM6 method AmberTools1634, 35 was used. The MP2
potential surface was generated using MOLPRO8, 36. The screened Tersoff
potential was provided by Atomistica37. The DFT potentials with the BLYP and
PBE functionals were, respectively, provided by CP2K38, 39 and CASTEP40, 41,
using the QUIP interface42, 43. GAP model was called via QUIP43.

In all cases the geometry optimisation was performed from within ASE44. The
other software packages were only used to compute the energy and gradient of
the configurations. A Python implementation of the FF and Exp+FF
preconditioners with several potential forms of nonbonded terms is available
within ASE44, 45 (https://gitlab.com/molet/ase).

##  4 Results

####  4.0.1 Organic molecules in gas phase

We investigated several organic molecules’ geometry optimisations with (FF)
and without (ID) our new FF-based preconditioner on three potential energy
surfaces (PM6, DFT and MP2), using the GAFF force field for building the
preconditioner. The results are shown in Table 1 and Figure 1. Convergence
criterion of geometry optimisations was $||\nabla E||=10^{-3}$ eV Å-1 for DFT
and MP2 surfaces while for the relatively inexpensive PM6 potential we applied
a slightly tighter threshold of $||\nabla E||=10^{-4}$ eV Å-1. Depending on
the system and underlying potential we can observe a 4–10 fold decrease in the
required number of optimisation steps using our preconditioner.

System (# of atoms) | PM6 | DFT(BLYP) | MP2/6-31G
---|---|---|---
ID | FF/GAFF | ID | FF/GAFF | ID | FF/GAFF | FF/LFF | Lindh
5-nitrobenzisoxazole (16) | 89 | 25 | 119 | 63 | 71 | 27 | 31 | 32
menthone (29) | 207 | 29 | 197 | 48 | 107 | 22 | 31 | 38
alanine tripeptide (32) | 395 | 77 | 536 | 124 | 210 | 59 | 47 | 67
thc (53) | 720 | 84 | 239 | 69 |  |  |  |
heme (75) | 500 | 175 | 358 | 95 |  |  |  |
taxol (113) | 1662 | 419 |  |  |  |  |  |
16-mer polyalanine (172) | 3549 | 348 |  |  |  |  |  |
|  |  |  |  |  |  |  |
Table 1: Total number of function/gradient calls of geometry optimisation for
organic molecules in gas phase using conventional (ID) and FF-based
preconditioned (FF) LBFGS method on three different quantum chemistry
surfaces. Convergence threshold was $||\nabla E||=10^{-4}$ eV Å-1 for PM6 and
$||\nabla E||=10^{-3}$ eV Å-1 for DFT and MP2 potentials, respectively. Figure
1: Computational saving of FF/GAFF preconditioner over unpreconditioned LBFGS
optimisations for geometry optimisation of molecular systems in gas phase on
three model potentials.

For PM6 only, to highlight the correlation between performance gain and ill-
conditioning, we also computed the ratio between the condition numbers for the
unpreconditioned and preconditioned Hessians, $\kappa_{I}/\kappa_{P}$, at the
minima, where

$\kappa_{P}=\frac{\lambda_{P}^{\max}}{\lambda_{P}^{\min}}=\frac{\max\limits_{u^{T}Pu=1}u^{T}Hu}{\min\limits_{u^{T}Pu=1}u^{T}\tilde{H}u}$
(14)

where $\tilde{H}$ is a modified Hessian where the zero eigenvalue due to
symmetries are removed. In Figure 2 we observe that the computational saving
is more strongly correlated to the condition number ratio than to the system
size.

Figure 2: Correlation between the computational saving and the condition
number ratio at potential energy minima of molecular systems in gas phase on
PM6 potential.

For the three smallest systems and the PM6 surface only we also compared our
FF-based preconditioner against using the exact Hessian of the model potential
(Hessian/GAFF) and a finite-difference hessian of PM6 (Hessian/PM6) as
preconditioners for LBFGS method. Zero eigenvalues of the Hessian matrices
were shifted to a moderate positive number to avoid numerical instabilities.
The results are collected in Table 2 and Figure 3. Our FF-based clearly
preconditioner outperforms both of these variants.

System (# of atoms) | LBFGS
---|---
ID | FF/GAFF | Hessian/GAFF | Hessian/PM6
5-nitrobenzisoxazole (16) | 89 | 25 | 32 | 24
menthone (29) | 207 | 29 | 36 | 21
alanine tripeptide (32) | 395 | 77 | 109 | 110
thc (53) | 720 | 84 | 95 | 120
|  |  |  |
Table 2: Comparison of total number of function/gradient calls of geometry
optimisation of different optimisation algorithms for minimisation of small
organic molecules on PM6 surface: unpreconditioned LBFGS (ID), FF-based
preconditioned LBFGS (FF/GAFF), FF-Hessian based preconditioned LBFGS
(Hessian/GAFF), PM6-Hessian based preconditioned LBFGS (PM6/GAFF). Convergence
threshold was $||\nabla E||=10^{-4}$ eV Å-1 for all cases. Figure 3:
Computational gain of different optimisation algorithms compared to the
unpreconditioned LBFGS method for geometry optimisation of molecular systems
in gas phase on PM6 potential.

We also examined the effect of using different force fields from which to
construct the FF-based preconditioner. Beside the GAFF force field, we
investigated two general force fields: the universal force field (UFF) 22 and
a general force field introduced by Lindh et al. (LFF) 7. The results shown in
Table 3 and Figure 4 indicate that there is no significant difference between
the three force fields. We only mention that the LFF force-field includes all
possible 2, 3 and 4–body interactions 7, resulting in a dense preconditioner
matrix, which for larger systems and an efficient potential energy surface
could become a performance bottleneck. By contrast, the preconditioners based
on GAFF or UFF are sparse, hence their cost scales linearly with system size.

System (# of atoms) | ID | FF/GAFF | FF/UFF | FF/LFF
---|---|---|---|---
5-nitrobenzisoxazole (16) | 89 | 25 | 31 | 34
menthone (29) | 207 | 29 | 39 | 39
alanine tripeptide (32) | 395 | 77 | 79 | 90
thc (53) | 720 | 84 | 77 | 92
heme (75) | 500 | 175 | 219 | n.a.
taxol (113) | 1662 | 419 | 393 | 400
16-mer polyalanine (172) | 3549 | 348 | 280 | 273
|  |  |  |
Table 3: Comparison of the effect of different force field based
preconditioners for the geometry optimisation of organic molecules in gas
phase (total number function/gradient calls). Convergence threshold was
$||\nabla E||=10^{-4}$ eV Å-1 for all cases. Figure 4: Computational gain of
the FF-based preconditioner using different force fields for geometry
optimisation of molecular systems in gas phase on PM6 potential.

Finally, we tested how the different preconditioners perform when applied to
transition state search, comparing again against ID (no preconditioning) and
against the Exp preconditioner (with default parameter set and $\mu=1$; unlike
LBFGS, CG is invariant under rescaling of $\mu$). The results are collected in
Table 4 and Figure 5. Overall the Exp preconditioner does not improve
significantly over ID. We experimented with different parameters, e.g., adding
connectivity information up to the 4-body interaction, but observed no
improvements. Both FF-based preconditioners are again comparable and yield a
much improved convergence even for these relatively small systems. For
instance, the gain is already 2–3-fold for dimethyl-phosphate and tyrosine
hydrolyses.

System (# of atoms) | ID | Exp | FF/GAFF | FF/LFF
---|---|---|---|---
HCCH $\leftrightarrow$ CCH2 (4) | 20 (79) | 17 (67) | 14 (55) | 32 (124)
H2CO $\leftrightarrow$ H2 \+ CO (4) | 24 (94) | 19 (73) | 18 (69) | 18 (70)
CH3O- $\leftrightarrow$ CH2OH- (5) | 18 (69) | 18 (70) | 14 (55) | 22 (82)
vinyl alcohol $\leftrightarrow$ acetaldehyde (7) | 58 (227) | 62 (245) | 39 (154) | 46 (179)
ring opening of cyclopropyl (8) | 52 (205) | 51 (199) | 31 (121) | 35 (137)
ring opening of bicyclo[1.1.0] butane TS 1 (10) | 87 (332) | 77 (291) | 56 (207) | 52 (191)
ring opening of bicyclo[1.1.0] butane TS 2 (10) | 66 (259) | 72 (282) | 35 (135) | 33 (127)
dimethyl-phosphate + OH- TS 1 (15) | 395 (1541) | 361 (1425) | 127 (504) | 123 (489)
dimethyl-phosphate + OH- TS 2 (15) | 355 (1379) | 329 (1302) | 172 (683) | 159 (631)
tyrosine + H2O (27) | 531 (2120) | 396 (1574) | 147 (583) | 185 (728)
|  |  |  |
Table 4: Number of steps of translations (and total number of function and
gradient calls in parentheses) of saddle searches using different
preconditioned variants of superlinearly converging dimer method. Convergence
threshold was $||\nabla E||=10^{-4}$ eV Å-1 for all cases. Figure 5:
Computational gain of the Exp and FF-based preconditioners compared to the
unpreconditioned superlinearly converging dimer method for molecular systems
in gas phase on PM6 potential.

###  4.1 Molecular crystals

We compared geometry optimisation with fixed unit cells using LBFGS,
preconditioned with ID (unpreconditioned), Exp3 and FF(GAFF force field). For
Exp the nearest neighbour distance ($r_{\mathrm{nn}}$) in Eq. 4 was calculated
from the initial structure, we specified $r_{\mathrm{cut}}=2r_{\mathrm{nn}}$
and $A=3.0$. In addition, we also employed the Exp+FF preconditioner as
defined in (11). The results for different molecular crystals are shown in
Table 5 and Figure 6. As expected, Exp reduces the number of optimisation
steps compared to ID, although the improvement is significantly smaller for
molecular crystals that for material systems3. Interestingly, FF alone already
leads to a significant speed-up over both ID and Exp even though the inter-
molecular interaction is not captured well. This indicates that for molecular
crystal optimisations preconditioning based on specific intramolecular
information is crucial. The most successful method was the Exp+FF combination,
which leads to a 3-7 fold speed up even for these relatively small test
systems.

System (# of atoms) | ID | Exp | FF/GAFF | Exp+FF/GAFF
---|---|---|---|---
xxii (60) | 77 | 60 | 45 | 29
xxi (84) | 291 | 164 | 134 | 77
xx (220) | 174 | 169 | 16 | 45
xix (112) | 193 | 137 | 73 | 29
xviii (184) | 232 | 97 | 70 | 39
|  |  |  |
Table 5: Total number of function/gradient calls of geometry optimisation of
molecular crystals using different preconditioning strategies. Convergence
criterion was $\|\nabla E\|=1.0^{-3}$ eV Å-1. Figure 6: Computational gain of
the Exp, FF and Exp+FF preconditioned over unpreconditioned LBFGS for geometry
optimisation of molecular crystals on a DFT potential.

###  4.2 Material systems

Finally, it is also interesting to investigate how an FF-based preconditioner
compares against the Exp preconditioner for material systems, where Exp
performs very well3. We tested geometry optimisation of bulk silicon and a
vacancy in bulk silicon with perturbed initial conditions, with increasing
system size. The screened Tersoff potential was used as the potential energy,
while the FF preconditioner was constructed from UFF. The results are shown in
Table 6. In both cases, the FF-based preconditioner yields a clear further
speed-up over Exp for both systems.

System (# of atoms) | Bulk geometry optimisation | Vacancy geometry optimisation
---|---|---
ID | Exp | FF/UFF | ID | Exp | FF/UFF
$2\times 2\times 2$ (64) | 32 | 17 | 10 | 34 | 15 | 9
$4\times 4\times 4$ (512) | 63 | 18 | 10 | 57 | 16 | 9
$8\times 8\times 8$ (4096) | 105 | 21 | 13 | 96 | 17 | 11
$16\times 16\times 16$ (32768) | 147 | 35 | 21 | 142 | 19 | 11
|  |  |  |  |  |
Table 6: Total number of function/gradient calls geometry optimisation steps
of bulk silicon and a bulk silicon vacancy using different preconditioning
strategies. Convergence criterion was $\|\nabla E\|=1.0^{-3}$ eV Å-1.

Another test system was bulk tungsten and a single interstitial site in bulk
tungsten with perturbed initial conditions and different system sizes. The
potential energy surface was provided by a GAP model that was trained on DFT
data. The preconditioner was based on a simple EAM potential:

$V_{\mathrm{EAM}}=\sum_{i}\left[\frac{1}{2}\sum_{j\neq
i}\Phi(r_{ij})+F\left(\sum_{j\neq i}\rho(r_{ij})\right)\right]$ (15)

where $\Phi(r_{ij})$ is a pair potential, $F$ is the embedding function and
$\rho({r_{ij}})$ is the electron charge density contribution from atom $j$ to
atom $i$. Based on equation 10 our FF-based preconditioner was defined as
$P_{\mathrm{FF}}=\sum_{\alpha}\frac{\partial
r_{\alpha}}{\partial\mathbf{r}}\otimes\frac{\partial
r_{\alpha}}{\partial\mathbf{r}}\left|\frac{\partial^{2}V_{\mathrm{EAM}}}{\partial
r_{\alpha}^{2}}\right|$, where $\alpha$ runs over all $ij$ pairs. In the
actual implementation $\Phi$, $F$ and $\rho$ functions are represented by
splines so computing the corresponding curvature is fairly straightforward.

The results are presented in Table 7. For both the bulk and interstitial
systems the number of function/gradient calls of the unpreconditioned
optimisation increases with system size while the preconditioned optimisations
require almost the same number of optimisation steps to achieve the same
convergence criterion.

System (# of atoms) | Bulk geometry optimisation | Interstitial geometry optimisation
---|---|---
ID | Exp | FF/EAM | ID | Exp | FF/EAM
$4\times 4\times 4$ (64) | 21 | 11 | 9 | 62 | 42 | 33
$8\times 8\times 8$ (512) | 32 | 12 | 8 | 72 | 43 | 36
$16\times 16\times 16$ (4096) | 56 | 12 | 9 | 116 | 46 | 33
|  |  |  |  |  |
Table 7: Total number of function/gradient calls for geometry optimisation of
bulk tungsten and interstitial defect in bulk tungsten using different
preconditioning strategies. Convergence criterion was $\|\nabla E\|=1.0^{-3}$
eV Å-1.

##  5 Conclusion

We introduced a flexible preconditioner for molecular simulation based on
empirical potentials that are widely implemented in popular molecular
mechanical program packages. Our method, which can be considered a
generalisation of Lindh et al.7, decomposes the analytic Hessian of the
empirical potential and modifies individual components to ensure their
positivity. An advantage of this procedure is that it avoids the computation
of second derivatives of the collective variables (or internal coordinates).
The preconditioner yields significant improvements (at least 2 fold, and
typically 5 fold decrease in function/gradient calls compared to
unpreconditioned techniques), demonstrated thoroughly on a wide range of
systems including molecules in gas phase, molecular crystals and materials,
using different target potential energy surfaces (empirical, semiempirical and
ab initio) as well as different optimisation tasks (geometry optimisations and
saddle searches).

##  6 Acknowledgement

The authors thank Prof. Chris J Pickard for stimulating discussions. Computing
facilities were provided by the Centre for Scientific Computing of the
University of Warwick.

## References

  * Liu and Nocedal 1989 Liu, D. C.; Nocedal, J. _Mathematical Programming_ 1989, _45_ , 503–528
  * Bakken and Helgaker 2002 Bakken, V.; Helgaker, T. _The Journal of Chemical Physics_ 2002, _117_ , 9160–9174
  * Packwood et al. 2016 Packwood, D.; Kermode, J.; Mones, L.; Bernstein, N.; Woolley, J.; Gould, N.; Ortner, C.; Csányi, G. _The Journal of Chemical Physics_ 2016, _144_ , 164109
  * Császár and Pulay 1984 Császár, P.; Pulay, P. _Journal of Molecular Structure_ 1984, _114_ , 31 – 34
  * Vogel et al. 1993 Vogel, S.; Fischer, T. H.; Hutter, J.; Lüthi, H. P. _International Journal of Quantum Chemistry_ 1993, _45_ , 679–688
  * Baker 1993 Baker, J. _Journal of Computational Chemistry_ 1993, _14_ , 1085–1100
  * Lindh et al. 1995 Lindh, R.; Bernhardsson, A.; Karlström, G.; Malmqvist, P.-A. _Chemical Physics Letters_ 1995, _241_ , 423 – 428
  * Werner et al. 2012 Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. _WIREs Comput Mol Sci_ 2012, _2_ , 242–253
  * Bitzek et al. 2006 Bitzek, E.; Koskinen, P.; Gähler, F.; Moseler, M.; Gumbsch, P. _Phys. Rev. Lett._ 2006, _97_ , 170201
  * Fogarasi et al. 1992 Fogarasi, G.; Zhou, X.; Taylor, P. W.; Pulay, P. _Journal of the American Chemical Society_ 1992, _114_ , 8191–8201
  * Kästner and Sherwood 2008 Kästner, J.; Sherwood, P. _The Journal of Chemical Physics_ 2008, _128_ , 014106
  * Henkelman and Jónsson 1999 Henkelman, G.; Jónsson, H. _The Journal of Chemical Physics_ 1999, _111_ , 7010–7022
  * Stewart 2007 Stewart, J. J. P. _Journal of Molecular Modeling_ 2007, _13_ , 1173–1213
  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. _Phys. Rev._ 1964, _136_ , B864–B871
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. _Phys. Rev._ 1965, _140_ , A1133–A1138
  * Møller and Plesset 1934 Møller, C.; Plesset, M. S. _Phys. Rev._ 1934, _46_ , 618–622
  * Becke 1988 Becke, A. D. _Phys. Rev. A_ 1988, _38_ , 3098–3100
  * Lee et al. 1988 Lee, C.; Yang, W.; Parr, R. G. _Phys. Rev. B_ 1988, _37_ , 785–789
  * VandeVondele and Hutter 2007 VandeVondele, J.; Hutter, J. _The Journal of Chemical Physics_ 2007, _127_ , 114105
  * Lippert et al. 1997 Lippert, G.; Hutter, J.; Parrinello, M. _Molecular Physics_ 1997, _92_ , 477–488
  * Goedecker et al. 1996 Goedecker, S.; Teter, M.; Hutter, J. _Phys. Rev. B_ 1996, _54_ , 1703–1710
  * Rappe et al. 1992 Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. _Journal of the American Chemical Society_ 1992, _114_ , 10024–10035
  * Wang et al. 2004 Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. _Journal of Computational Chemistry_ 2004, _25_ , 1157–1174
  * Baker and Chan 1996 Baker, J.; Chan, F. _Journal of Computational Chemistry_ 1996, _17_ , 888–904
  * Bardwell et al. 2011 Bardwell, D. A.; Adjiman, C. S.; Arnautova, Y. A.; Bartashevich, E.; Boerrigter, S. X. M.; Braun, D. E.; Cruz-Cabeza, A. J.; Day, G. M.; Della Valle, R. G.; Desiraju, G. R.; van Eijck, B. P.; Facelli, J. C.; Ferraro, M. B.; Grillo, D.; Habgood, M.; Hofmann, D. W. M.; Hofmann, F.; Jose, K. V. J.; Karamertzanis, P. G.; Kazantsev, A. V.; Kendrick, J.; Kuleshova, L. N.; Leusen, F. J. J.; Maleev, A. V.; Misquitta, A. J.; Mohamed, S.; Needs, R. J.; Neumann, M. A.; Nikylov, D.; Orendt, A. M.; Pal, R.; Pantelides, C. C.; Pickard, C. J.; Price, L. S.; Price, S. L.; Scheraga, H. A.; van de Streek, J.; Thakur, T. S.; Tiwari, S.; Venuti, E.; Zhitkov, I. K. _Acta Crystallographica Section B_ 2011, _67_ , 535–551
  * Reilly et al. 2016 Reilly, A. M.; Cooper, R. I.; Adjiman, C. S.; Bhattacharya, S.; Boese, A. D.; Brandenburg, J. G.; Bygrave, P. J.; Bylsma, R.; Campbell, J. E.; Car, R.; Case, D. H.; Chadha, R.; Cole, J. C.; Cosburn, K.; Cuppen, H. M.; Curtis, F.; Day, G. M.; DiStasio Jr, R. A.; Dzyabchenko, A.; van Eijck, B. P.; Elking, D. M.; van den Ende, J. A.; Facelli, J. C.; Ferraro, M. B.; Fusti-Molnar, L.; Gatsiou, C.-A.; Gee, T. S.; de Gelder, R.; Ghiringhelli, L. M.; Goto, H.; Grimme, S.; Guo, R.; Hofmann, D. W. M.; Hoja, J.; Hylton, R. K.; Iuzzolino, L.; Jankiewicz, W.; de Jong, D. T.; Kendrick, J.; de Klerk, N. J. J.; Ko, H.-Y.; Kuleshova, L. N.; Li, X.; Lohani, S.; Leusen, F. J. J.; Lund, A. M.; Lv, J.; Ma, Y.; Marom, N.; Masunov, A. E.; McCabe, P.; McMahon, D. P.; Meekes, H.; Metz, M. P.; Misquitta, A. J.; Mohamed, S.; Monserrat, B.; Needs, R. J.; Neumann, M. A.; Nyman, J.; Obata, S.; Oberhofer, H.; Oganov, A. R.; Orendt, A. M.; Pagola, G. I.; Pantelides, C. C.; Pickard, C. J.; Podeszwa, R.; Price, L. S.; Price, S. L.; Pulido, A.; Read, M. G.; Reuter, K.; Schneider, E.; Schober, C.; Shields, G. P.; Singh, P.; Sugden, I. J.; Szalewicz, K.; Taylor, C. R.; Tkatchenko, A.; Tuckerman, M. E.; Vacarro, F.; Vasileiadis, M.; Vazquez-Mayagoitia, A.; Vogt, L.; Wang, Y.; Watson, R. E.; de Wijs, G. A.; Yang, J.; Zhu, Q.; Groom, C. R. _Acta Crystallographica Section B_ 2016, _72_ , 439–459
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Vanderbilt 1990 Vanderbilt, D. _Phys. Rev. B_ 1990, _41_ , 7892–7895
  * Tersoff 1986 Tersoff, J. _Phys. Rev. Lett._ 1986, _56_ , 632–635
  * Pastewka et al. 2013 Pastewka, L.; Klemenz, A.; Gumbsch, P.; Moseler, M. _Phys. Rev. B_ 2013, _87_ , 205410
  * Szlachta et al. 2014 Szlachta, W. J.; Bartók, A. P.; Csányi, G. _Phys. Rev. B_ 2014, _90_ , 104108
  * Daw and Baskes 1984 Daw, M. S.; Baskes, M. I. _Phys. Rev. B_ 1984, _29_ , 6443–6453
  * Zhou et al. 2004 Zhou, X. W.; Johnson, R. A.; Wadley, H. N. G. _Phys. Rev. B_ 2004, _69_ , 144113
  * Case et al. 2016 Case, D. A.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E.; III,; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Luchko, T.; Luo, R.; Madej, B.; Mermelstein, D.; Merz, K. M.; Monard, G.; Nguyen, H.; Nguyen, H. T.; Omelyan, I.; Onufriev, A.; Roe, D. R.; Roitberg, A.; Sagui, C.; Simmerling, C. L.; Botello-Smith, W. M.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; Xiao, L.; Kollman, P. A. _AMBER 16_ ; 2016
  * 35 http://www.ambermd.org
  * Werner et al. 2015 Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; P., M. S.; Celani,; Györffy, W.; Kats, D.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M. MOLPRO, version 2015.1, a package of ab initio programs. 2015; https://www.molpro.net
  * 37 http://www.atomistica.org
  * VandeVondele et al. 2005 VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. _Computer Physics Communications_ 2005, _167_ , 103 – 128
  * cp 2 http://www.cp2k.org
  * Segall et al. 2002 Segall, M. D.; Lindan, P. J. D.; Probert, M. J.; Pickard, C. J.; Hasnip, P. J.; Clark, S. J.; Payne, M. C. _Journal of Physics: Condensed Matter_ 2002, _14_ , 2717
  * 41 http://www.castep.org
  * Csányi et al. 2007 Csányi, G.; Winfield, S.; Kermode, J.; Comisso, A.; Vita, A. D.; Bernstein, N.; ; Payne, M. C. _IoP Comput. Phys. Newsletter_ 2007, _Spring_
  * 43 http://www.libatoms.org
  * Bahn and Jacobsen 2002 Bahn, S. R.; Jacobsen, K. W. _Computing in Science Engineering_ 2002, _4_ , 56–66
  * Larsen et al. 2017 Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Jensen, P. B.; Kermode, J.; Kitchin, J. R.; Kolsbjerg, E. L.; Kubal, J.; Kaasbjerg, K.; Lysgaard, S.; Maronsson, J. B.; Maxson, T.; Olsen, T.; Pastewka, L.; Peterson, A.; Rostgaard, C.; Schiøtz, J.; Schütt, O.; Strange, M.; Thygesen, K. S.; Vegge, T.; Vilhelmsen, L.; Walter, M.; Zeng, Z.; Jacobsen, K. W. _Journal of Physics: Condensed Matter_ 2017, _29_ , 273002