# Human-machine cooperation for semantic
feature listing

Kushin Mukherjee, Siddharth Suresh & Timothy Rogers
Department of Psychology
University of Wisconsin-Madison
Madison, WI 53706, USA
{ kmukherjee2, siddharth.suresh<EMAIL_ADDRESS>
###### Abstract

Semantic feature norms — lists of features that concepts do and do not possess
— have played a central role in characterizing human conceptual knowledge, but
require extensive human labor. Large language models (LLMs) offer a novel
avenue for the automatic generation of such feature lists, but are prone to
significant error. Here, we present a new method for combining a learned model
of human lexical-semantics from limited data with LLM-generated data to
efficiently generate high-quality feature norms.

Introduction. A central goal in cognitive science is to characterize human
knowledge of concepts and their properties. Many have used human-generated
feature lists as norms for establishing the structural relationship between
concepts in the human mind (McRae et al., 2005; Devereux et al., 2014; De
Deyne et al., 2008; Buchanan et al., 2019), but this requires extensive human
labor. Large language models (LLMs) have recently shown impressive
capabilities when generating properties of objects (Hansen & Hebart, 2022) or
answering questions(Ouyang et al., 2022; Brown et al., 2020; Hoffmann et al.,
2022; Chowdhery et al., 2022; Wei et al., 2021) and thus suggest an avenue for
more efficient characterization of human knowledge structures, but even state-
of-the-art models can routinely fail on many common-sense questions of fact.
GTP3-davinci, for instance, will deny that alligators are green, while
asserting that they can be used to suck dust up from surfaces. Thus, human
effort can generate high-quality norms, but with prohibitive costs, while LLMs
can produce norms with little human effort, but with considerably less
accuracy. This paper considers whether human and machine effort can combine to
efficiently estimate high-quality semantic feature vectors. Our approach
leverages two key observations. First, the concept-by-feature matrices
collected in norming experiments are typically high in dimension but low in
rank—thus a low-dimension matrix decomposition based on a representative
subset of concepts can be used to estimate how subsets of features covary with
one another across concepts. Second, open-source LLMs optimized for question-
answering, such as Google’s FLAN-T5 XXL (Wei et al., 2021), can be used to
generate “guesses” about which concepts possess which features. While these
guesses may not be fully accurate, they can be be used to infer where a new
concept might reside within the low-dimension similarity space computed from
human-generated features. The predicted coordinates of the new item can then
be re-composed to the original human feature space to make inferences about
its features. This idea suggests the following procedure:

  1. 1.

Collect human feature norms for n concepts representative of a domain and
assemble in concept-by-feature matrix $H$.

  2. 2.

Compute a singular value decomposition of $H$, retaining the matrix of the
first $d$ singular values times the first $d$ right singular vectors.

  3. 3.

For a new concept, use FLAN-T5 XXL to answer questions about which properties
are ‘true’ for the concept to create a binary vector of LLM-generated
features, $h$.

  4. 4.

Fit a logistic regression model to predict the new concept’s $d$ left singular
vector coordinates using the $d$ right singular vectors as predictors and $h$
as the target.

  5. 5.

Multiply out the predicted left singular vectors through the rest of the SVD
to get predictions about features of the new item.

Refer to Figure 1 A. for a visual schematic of the procedure. We evaluated
this approach using a well-known pre-existing dataset of semantic feature
norms De Deyne et al. (2008).

Method. The human norm dataset included 129 animal concepts evaluated on 764
features, and 166 tools evaluated on 1262 features. Each feature corresponded
to a property generated by human participants when asked to ’list all the
features they knew to be true of each concept’. For each domain separately, 4
human raters evaluated whether a given feature was generally true of a given
concept (e.g. “Yes or no, is the property is green true of the concept
alligator”). We threshold the values in these matrices such that only features
that all 4 raters agreed to be true were treated as valid ’trues’. This gave
us a set of ground truth binary human-generated feature matrices. Next, we
used the XXL variant of FLAN-T5, a 11B parameter language model accessed
through the huggingface transformers library (Wolf et al., 2019), to perform
the same feature verification task for every cell in both matrices (see
Appendix). We then treated human/machine feature comparisons as a signal
detection problem in which the human norms provided the ground truth (signal
present/absent) against which the machine responses were evaluated as hits,
misses, false-alarms, or correct rejections. From these data we computed
$d^{\prime}$ as the measure of agreement between human and machine matrices.
This metric was computed separately for the animal and tool datasets using
just the FLAN-verified feature vectors. The central question was whether the
matrix completion-based data-augmentation approach described above improves
the match between machine and human feature vectors.

Results. Figure 1 A illustrates the general workflow. To test if our matrix
completion technique improved alignment between human and machine features, we
iteratively held-out one concept at a time and computed the SVD and regression
models on the remaining concepts. We then estimated the features for the held-
out concept and computed its $d^{\prime}$ relative to ground-truth human
features. We found that our matrix completion technique led to higher
$d^{\prime}$ values relative to using the raw features generated by the LLM
for both animal ($t=15.86$, $p<.001$) and tool ($t=3.79$, $p<.001$) concepts
(Figure 1B.).

Figure 1: (A.) Our procedure for estimating features for new concepts. A human
feature-listing matrix is decomposed using (SVD) into low-rank matrices of
principal components (right) and weights to linearly combine them to estimate
the features for each item (left). For a new concept, logistic regression
models estimate weights to map the principal components to LLM-generated
features for the concept. Multiplying the estimated weights with the principal
components improves the alignment between humans and machines. (B.)
$d^{\prime}$ distributions for raw LLM features and LLM features generated
from the process in (A). (C.) The difference in $d^{\prime}$ values using raw
LLM features and the method in (A) with different amounts of the original
dataset withheld.

For this method to be efficient, it should require little data to fit the SVD
model. To estimate how much data would be needed for our matrix completion
method to show meaningful gains, we held out 10% - 90% of the original human
data when fitting the SVD model and estimated the differences in $d^{\prime}$
when using our method vs. raw LLM feature lists. We found statistically
significant improvements using our method while holding out up to 70% of the
original data, with much greater improvements for animals than tools (Figure 1
C.).

Conclusion. We evaluated a new method for combining human effort and LLM
behavior to more efficiently generate semantic feature norms. Such norms have
been critical for understanding the nature and structure of human conceptual
knowledge, but prior efforts have required extensive human labor. Our method
yields higher-quality vectors than those produced by LLMs alone, with
potentially huge savings in human effort–opening the door for mapping a much
broader set of concepts in future.

#### URM Statement

The authors acknowledge that at least one key author of this work meets the
URM criteria of ICLR 2023 Tiny Papers Track.

## References

  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901, 2020.
  * Buchanan et al. (2019) Erin M Buchanan, Kathrene D Valentine, and Nicholas P Maxwell.  English semantic feature production norms: An extended database of 4436 concepts.  _Behavior Research Methods_ , 51:1849–1863, 2019.
  * Chowdhery et al. (2022) Aakanksha Chowdhery, Sharan Narang, Jacob Devlin, Maarten Bosma, Gaurav Mishra, Adam Roberts, Paul Barham, Hyung Won Chung, Charles Sutton, Sebastian Gehrmann, et al.  Palm: Scaling language modeling with pathways.  _arXiv preprint arXiv:2204.02311_ , 2022.
  * De Deyne et al. (2008) Simon De Deyne, Steven Verheyen, Eef Ameel, Wolf Vanpaemel, Matthew J Dry, Wouter Voorspoels, and Gert Storms.  Exemplar by feature applicability matrices and other dutch normative data for semantic concepts.  _Behavior research methods_ , 40:1030–1048, 2008.
  * Devereux et al. (2014) Barry J Devereux, Lorraine K Tyler, Jeroen Geertzen, and Billi Randall.  The centre for speech, language and the brain (cslb) concept property norms.  _Behavior research methods_ , 46:1119–1127, 2014.
  * Hansen & Hebart (2022) Hannes Hansen and Martin N Hebart.  Semantic features of object concepts generated with gpt-3.  _arXiv preprint arXiv:2202.03753_ , 2022.
  * Hoffmann et al. (2022) Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al.  Training compute-optimal large language models.  _arXiv preprint arXiv:2203.15556_ , 2022.
  * McRae et al. (2005) Ken McRae, George S Cree, Mark S Seidenberg, and Chris McNorgan.  Semantic feature production norms for a large set of living and nonliving things.  _Behavior research methods_ , 37(4):547, 2005\.
  * Ouyang et al. (2022) Long Ouyang, Jeff Wu, Xu Jiang, Diogo Almeida, Carroll L Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Ray, et al.  Training language models to follow instructions with human feedback.  _arXiv preprint arXiv:2203.02155_ , 2022.
  * Wei et al. (2021) Jason Wei, Maarten Bosma, Vincent Y Zhao, Kelvin Guu, Adams Wei Yu, Brian Lester, Nan Du, Andrew M Dai, and Quoc V Le.  Finetuned language models are zero-shot learners.  _arXiv preprint arXiv:2109.01652_ , 2021.
  * Wolf et al. (2019) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, et al.  Huggingface’s transformers: State-of-the-art natural language processing.  _arXiv preprint arXiv:1910.03771_ , 2019.

##  Appendix A Appendix

###  A.1 FLAN-T5 Feature verification

We used a standardized prompt for querying FLAN-T5 about whether a given
feature was applicable for each concept in the animal and tool concept sets.
Table 1 shows an example of a prompt to query whether the feature ’has two
eyes’ is valid for the concept ’Dolhpin’. We prefaced the question with the
same example of a positive and negative example for the concept ’book’. We
repeated this procedure for every concept x feature combination in the norm
dataset collected by De Deyne et al. (2008)

Q: Is the property [is_female] true for the concept [book]?
---
A: False
Q: Is the property [can_be_digital] true for the concept [book]
A: True
In one word True/False, answer the following question
Q: Is the property [has_two_eyes] true for Dolphins?
A: $<$mask$>$
Table 1: Example Prompt.

###  A.2 $d^{\prime}$ as an alignment metric

Since both the ground-truth human feature matrices and LLM-generated feature
matrices had binary entries, the problem of human-machine comparison can be
described in terms of signal detection theory. That is, if we treat the human-
matrix as being the source of ’signal’, the predictions of 1s and 0s in the
machine matrix can be — (1) Hits if the cell in the matrix was 1 for both the
human and machine matrix, (2) Misses if the cell in the machine matrix is 0
and 1 in the human matrix, (3) False alarms if the cell in the machine matrix
is 1 and 0 in the human matrix, and (4) Correct rejections if the cell in both
the machine and human matrices is 0. The number of hits, misses, false alarms,
and correct rejections can be tallied to compute hit rate (HR) and false-alarm
rate (FAR) as follows -

$HR=\frac{hits}{hits+misses}$ (1) $FAR=\frac{false\ alarms}{false\
alarms+correct\ rejections}$ (2)

Finally, $d^{\prime}$ is computed as

$d^{\prime}=z(H)-z(FAR),$ (3)

where $z$ is the inverse of the cumulative distribution function (CDF) of the
standard normal distribution $\mathcal{N}(0,1)$. A higher $d^{\prime}$
indicates a greater degree of alignment between human and machine features.

Figure 2: Profile of singular values for both human feature matrices
visualized as scree plots. As can be observed, these matrices were low-rank,
thus when fitting our feature verification model, we used a rank-10
decomposition when computing the SVD. We did find minor variations in
performance as a function of changing the rank of the decomposition, but
overall as long as the rank was around 10, which corresponds to the elbow of
the scree plots, the reconstructed values were accurate. Figure 3:
$d^{\prime}$ distributions for different holdout percentages. As the amount of
data used to fit the matrix-completion model is reduced, the more the model-
based distribution overlaps with the pure LLM distribution (shown in blue).
While greater gains are seen for the animal concepts, FLAN is already quite
good at estimating the structure of tool concepts. Figure 4: Human feature-
ratings, FLAN-based LLM generated feature-ratings, and features generated
using our matrix-completion model using a leave-one-out approach for (A)
animal concepts and (B) tool concepts . Yellow cells correspond to 1s and
purple cells correspond to 0s. The LLM tends to provide many false positives
as is indicated by the greater proportion of yellow cells. Our method helps to
’clean’ the LLM generated matrix to bring it into closer alignment with human
ratings.


# Self-interaction corrected SCAN functional for molecules and solids in the
numeric atom-center orbital framework

Sheng Bi  The NOMAD Laboratory at the FHI of the Max-Planck-Gesellschaft and
IRIS-Adlershof of the Humboldt-Universität zu Berlin, Faradayweg 4-6, D-14195
Berlin-Dahlem Germany  Department of Chemistry, Fudan University, Shanghai
200433, People’s Republic of China  Research Center for Intelligent
Supercomputing, Zhejiang Lab, Hangzhou 311100, P. R. China  Christian Carbogno
<EMAIL_ADDRESS>The NOMAD Laboratory at the FHI of the
Max-Planck-Gesellschaft and IRIS-Adlershof of the Humboldt-Universität zu
Berlin, Faradayweg 4-6, D-14195 Berlin-Dahlem Germany  Igor Ying Zhang
<EMAIL_ADDRESS>Department of Chemistry, Fudan University,
Shanghai 200433, People’s Republic of China  MOE Key Laboratory of
Computational Physical Sciences, Shanghai Key Laboratory of Bioactive Small
Molecules  Matthias Scheffler  The NOMAD Laboratory at the FHI of the Max-
Planck-Gesellschaft and IRIS-Adlershof of the Humboldt-Universität zu Berlin,
Faradayweg 4-6, D-14195 Berlin-Dahlem Germany

###### Abstract

Semilocal density-functional approximations (DFAs), including the state-of-
the-art SCAN functional, are plagued by the self-interaction error (SIE).
While this error is explicitly defined only for one-electron systems, it has
inspired the self-interaction correction method proposed by Perdew and Zunger
(PZ-SIC), which has shown promise in mitigating the many-electron SIE.
However, the PZ-SIC method is known for its significant numerical instability.
In this study, we introduce a novel constraint that facilitates self-
consistent localization of the SIC orbitals in the spirit of Edmiston-
Ruedenberg orbitals [Rev. Mod. Phys. 35, 457 (1963)]. Our practical
implementation within the all-electron numeric atom-centered orbitals code
FHI-aims guarantees efficient and stable convergence of the self-consistent
PZ-SIC equations for both molecules and solids. We further demonstrate that
our PZ-SIC approach effectively mitigates the SIE in the meta-GGA SCAN
functional, significantly improving the accuracy for ionization potentials,
charge-transfer energies, and band gaps for a diverse selection of molecules
and solids. However, our PZ-SIC method does have its limitations. It can not
improve the already accurate SCAN results for properties such as cohesive
energies, lattice constants, and bulk modulus in our test sets. This
highlights the need for new-generation DFAs with more comprehensive
applicability.

##  I Introduction

Kohn-Sham (KS) density-functional theory (DFT) [1, 2] is currently the most
widely used electronic-structure method across various scientific disciplines
[3, 4, 5]. To a large extent, this is due to the favorable balance between
computational accuracy and efficiency. The key contribution to this favorable
balance is the existence of density-functional approximations (DFAs) to the
exact exchange-correlation (XC) functional [2, 5]. In this context, the non-
empirical SCAN functional proposed by Sun, Ruzsinszky, and Perdew in 2015 [6]
is particularly promising. Compared to earlier semi-local DFAs, including the
local-density approximations (LDAs), the generalized gradient approximations
(GGAs), and other meta-GGAs [7, 8, 9], the SCAN meta-GGA functional yields
consistent and notable improvements in describing many chemical and physical
properties of both molecules and solids [10, 11, 12]. This can be attributed
to its compliance with 17 constraints, which the exact XC functional must
satisfy. However, like all other semi-local DFAs, the SCAN functional is not
immune to the self-interaction error (SIE) [13, 14, 15]. This limitation
accounts for its inability to reliably describe the ionization potential,
specifically the energy of the highest occupied molecular orbital, as well as
charge-transfer properties and band gaps in insulators and semiconductors [16,
17].

Formally, the SIE arises from the ground of the XC functional within the KS-
DFT framework, i.e., the Hartree term 111Hartree atomic units are used
throughout the paper.:

$\displaystyle
E_{\text{H}}[n]:=\frac{1}{2}\iint\frac{n(\textbf{r})n(\textbf{r}^{\prime})}{|\textbf{r}-\textbf{r}^{\prime}|}d^{3}\textbf{r}^{\prime}d^{3}\textbf{r}\quad,$
(1)

which captures the most part of electron-electron Coulomb interaction in the
system. However, since the density $n(\textbf{r})$ encompasses all electrons,
the Hartee term also introduces non-physical Coulomb interactions where
individual electrons interact with themselves. In one-electron systems, this
spurious interaction is precisely counterbalanced by the Fock term within the
Hartree-Fock (HF) theory. In other words, the HF method is free from one-
electron self-interactions. However, it is well-documented that the HF
approach is not enough to cancel out these non-physical self-interactions in
multi-electron systems [19]. Within the KS-DFT framework, these non-physical
interactions are expected to be precisely canceled out by the exact XC
functional.

Unfortunately, this cancellation is insufficiently achieved for all existing
semi-local DFAs. The residual self-interaction implies the SIE of a given DFA,
leading to incorrect electron delocalization [13, 20]. Obviously, a single
orbital’s XC energy must cancel its self-interacted Coulomb interaction. This
analysis has motivated the development of so-called self-interaction
corrections (SIC), initially proposed by Perdew and Zunger in the early 1980s
(PZ-SIC) [15]:

$\displaystyle E^{\text{PZ-
SIC}}[\\{n_{i}\\}]:=-\sum_{i}^{N_{\text{e}}}\left(E_{\text{xc}}^{\text{DFA}}[n_{i}]+E_{\text{H}}[n_{i}]\right)\quad.$
(2)

Here, $N_{e}$ represents the total number of electrons in the system, and
$\\{n_{i}(\textbf{r})=|\phi_{i}(\textbf{r})|^{2}\\}$ denotes the densities of
occupied single-electron orbitals, referred to as SIC orbitals
$\\{\phi_{i}(\textbf{r})\\}$ henceforth (for further details, see Section II).
Accordingly, $E_{\text{xc}}^{\text{DFA}}[n_{i}]$ is the DFA-specific single-
particle XC energy evaluated for the SIC orbital density $n_{i}(\textbf{r})$.
Meanwhile, $E_{\text{H}}[n_{i}]$ is the Coulomb self-interaction energy of
$n_{i}(\textbf{r})$, calculated using the Hartree equation (Eq. (1)). By
incorporating these one-electron correction terms into the semi-local DFA XC
energy, the PZ-SIC scheme effectively eliminates the SIE of a given DFA in
one-electron systems. Meanwhile, it has been shown that the PZ-SIC scheme can
substantially mitigate many-electron SIE when the appropriate SIC orbitals are
employed [21, 22, 23].

While SIC orbitals coincide with KS orbitals in one-electron systems, they
generally diverge in many-electron systems. Theoretically, the optimal SIC
orbitals should be determined by variationally minimizing the total energy of
the self-interaction corrected DFA (SIC-DFA) [15]. This minimization must
adhere to the total-density constraint, which stipulates that the sum of all
occupied SIC orbital densities should yield the total electron density
$n(\textbf{r})=\sum_{i}^{N_{\text{e}}}n_{i}(\textbf{r})$ (for further details,
see Section II). Moreover, as initially proposed by Pederson, Heaton, and Lin,
the variational minimization of the SIC-DFA total energy not only produces the
PZ-SIC equations for KS orbitals but also imposes the so-called orbital
potential constraint on SIC orbitals [24]. Similar to the standard KS
equations, PZ-SIC equations should be solved iteratively in a self-consistent
manner because they depend on both KS and SIC orbitals (also see Section II
for more details).

In practice, the PZ-SIC equations often yield multiple self-consistent
solutions for SIC orbitals, all satisfying the total density and orbital
potential constraints. As a result, while the correct SIC orbitals that
minimize the SIC-DFA total energy are among these multiple self-consistent
solutions, the specific solution obtained is highly sensitive to the choice of
initial guess.

In the seminal work of PZ-SIC, Perdew and Zunger observed that localized
orbitals often provide lower SIC-DFA total energies compared to their
delocalized (canonical) counterparts, such as standard KS orbitals [15]. This
observation led to frequently characterizing SIC orbitals as “localized
orbitals”. This insight has inspired several empirical guidelines for
initializing SIC orbitals [25]. A variety of localization methods have been
used to generate localized initializations [26], including Pipek-Mezey [27],
Edmiston-Ruedenberg [28], von Niessen [29], Foster-Boys [30], and the fourth-
moment [31] methods. Lehtola, Head-Gordon, and Jónsson [26] systematically
examined various initial SIC orbitals for atoms in the first three rows of the
periodic table. Remarkably, the SIC total energies associated with these
different sets of orbitals can diverge by more than 0.5 eV. Consequently,
identifying the correct SIC orbitals remains a significant challenge, even for
atoms and simple diatomic molecules [25, 24].

More recently, Pederson, Ruzsinszky, and Perdew noted the topological
similarity between Fermi orbitals with Löwdin orthogonalization (FLOs) and
previously employed localized orbitals. Therefore, they suggested constructing
the SIC orbitals using FLOs [32, 33, 34, 35], giving rise to a modified PZ-SIC
algorithm known as FLO-SIC. This innovative approach effectively filters out
the FLO-like SIC orbitals from the multiple self-consistent solutions,
mitigating the issue of multiple solutions. FLO-SIC has demonstrated
considerable promise for atoms and molecules [36, 37, 38], with particular
success observed in systems containing $\pi$ bonds [36, 39].

For extended materials, finding the correct SIC orbitals becomes even more
serious [40]. To capture the localized nature of SIC orbitals, Heaton,
Harrison, and Lin suggested representing them in terms of real-space Wannier
orbitals [22] under periodic boundary conditions (PBCs). However, it does not
thoroughly solve the issue of multiple solutions. Different sets of Wannier
orbitals can satisfy the total density and orbital-potential constraints.
Consequently, discrepancies arise; for instance, the SIC-LDA band gap of the
Ar crystal reported by Heaton, Harrison, and Lin differs from that of Stengel
and Spaldin as much as 0.7 eV [22, 41]. Large uncertainties are reported for
the band gaps of bulk Silicon [41, 23] and transition-metal mono-oxides (MnO,
FeO, CoO, and NiO) [42, 43]. The SIC orbitals of those results were all
filtered by satisfying the constraints of total density and orbital potential
but using different localization definitions for the Wannier orbitals.

More recently, in 2020, Shinde et al. introduced the use of conceptually
related FLOs for filtering out the Wannier orbitals (WFLO) [44]. As a result,
the WFLO-SIC approach identifies proper SIC orbitals for solids, leading to a
significant improvement in the band gap of the semi-local PBE functional [45].
Additionally, alternative localization strategies have an enormous potential
to establish the SIC orbitals [26, 23], including the maximally localized
Wanner orbitals (MLWOs) [46] or the Edmiston-Ruedenberg (E-R) type orbitals
[28]. The E-R orbitals are similar to MLWOs in molecular systems [47] and are
essentially identical to MLWOs in extended periodic systems [48]. However,
these works evaluate the SIC-DFA total energies in a non-self-consistent
manner. To our knowledge, there is no report on the self-consistent SIC-DFA
implementation that explicitly incorporates MLWOs or E-R localization.

Besides the SIC problems discussed above, there is another potential source of
inaccuracy in most of the discussed PZ-SIC studies, i.e., the usage of semi-
local DFAs at the level of LDA or GGA [14, 49, 50, 35, 39, 22, 22]. Since the
SIC is inherently linked to the employed DFA via Eq. (2), it is yet unclear to
which extent the observed trends are related to the SIC or to the DFA. This is
further substantiated by the fact that the combination of PZ-SIC and the
state-of-the-art meta-GGA SCAN has shown some promising results [32, 51, 52,
53]. Consequently, there is an urgent need for practical, self-consistent SIC
implementations compatible with the SCAN functional for both molecules and
solids.

In this work, we overcome these challenges by building upon the total-density
constraint of Perdew and Zunger and the orbital potential constraint of
Pederson, Heaton, and Lin. We introduce a novel third constraint, the orbital
density-potential constraint, to PZ-SIC equations. Inspired by the E-R
orbitals [28], our proposed constraint aims to maximize spatial localization
of the SIC orbitals under the orbital potential constraint. As demonstrated in
the subsequent sections, this approach facilitates efficient convergence of
the PZ-SIC equations for both molecular and extended systems.

We demonstrate the efficacy of our proposed methodology through its
implementation in the all-electron numeric atom-center orbital (NAO) package
FHI-aims [54, 55]. Besides applying our PZ-SIC approach for semi-local DFAs
such as LDA and GGA, we primarily focus on PZ-SIC calculations with the meta-
GGA SCAN functional (SIC-SCAN). For the latter, we observe significant
improvements in the highest occupied orbital energies when compared to the
ionization potential of atoms from H to Ar, as well as for 18 molecules taken
from reference [56]. Furthermore, the band gaps of various solids, including
Si, MnO, diamond-C, MgO, and LiF, are calculated using SIC-SCAN, yielding
comparable results to $GW$ calculations. The combination of PZ-SIC and SCAN
offers enhanced accuracy in electronic structure calculations while
maintaining a moderate computational burden, particularly for periodic solids.

The paper is organized as follows: Section II provides a derivation of the
standard PZ-SIC formula and introduces the three constraints used to solve the
PZ-SIC equations, both for finite and periodic systems. In particular, we
analyze and discuss the physical motivation behind the orbital density-
potential constraint by using the carbon atom, the Helium dimer (He2), the
Helium crystal, and the periodic polyethylene chain. In section III, we assess
the performance of PZ-SIC in conjunction with SCAN by examining various
metrics such as SIC energy, ionization potential, band gap, and cohesive
energy for a group of molecules and solids. Conclusions and outlook are
presented in section IV.

##  II Theory

With the PZ-SIC concept (Eq. (2)), the total energy of a self-interaction
corrected DFA, termed SIC-DFA, is written as [15]

$\displaystyle E^{\text{SIC-
DFA}}[n,\\{n_{i}\\}]:=E^{\text{DFA}}[n]+E^{\text{PZ-SIC}}[\\{n_{i}\\}]\quad.$
(3)

Here $E^{\text{DFA}}[n]$ is the total energy of a given DFA in the standard KS
framework [2]:

$\displaystyle
E^{\text{DFA}}[n]:=T_{\text{s}}[n]+E_{\text{ext}}[n]+E_{\text{H}}[n]+E^{\text{DFA}}_{\text{xc}}[n]\quad,$
(4)

where $T_{\text{s}}$ is the kinetic energy of the KS non-interacting
electrons. $E_{\text{ext}}[n]$ and $E_{\text{H}}[n]$ are the external
potential energy and the Hartree energy, respectively, both of which are
explicit functionals of the electron density $n(\textbf{r})$.
$E^{\text{DFA}}_{\text{xc}}[n]$ is the XC energy of the DFA.

In the standard KS-DFT scheme, the KS equations are obtained via a variational
minimization of the DFA energy $E^{\text{DFA}}[n]$ (Eq. (4)), whereby the
conservation of the number of electrons $N_{e}$ has to be ensured. This is
equivalent to enforcing orthonormalization on the canonical occupied KS
orbitals $\\{\psi_{l}\\}$ [15], leading to the Euler equation

$\displaystyle\delta$
$\displaystyle\left[E^{\text{DFA}}-\sum_{ij}^{N_{e}}\epsilon_{ij}\left(\langle\psi_{i}|\psi_{j}\rangle-\delta_{ij}\right)\right]=0\quad.$
(5)

Because the DFA energy $E^{\text{DFA}}[n]$ remains invariant under the unitary
transformations of the canonical occupied KS orbitals $\\{\psi_{l}\\}$, any
orthonormalized sets of single-particle orbitals $\\{\phi_{i}\\}$, derived
through a unitary transformation from $\\{\psi_{l}\\}$, satisfies Eq. (5) and
yields the minimal DFA energy $E^{\text{DFA}}[n]$. It is worth pointing out
that while the unitary transformation preserves the total density

$\displaystyle\sum_{i=1}^{N_{e}}\left|\phi_{i}(\textbf{r})\right|^{2}=\sum_{i=1}^{N_{e}}\left|\psi_{i}(\textbf{r})\right|^{2}=n(\textbf{r})\quad,$
(6)

it does not necessarily maintain the individual orbital densities
$\\{n_{i}(\textbf{r})\\}$

$\displaystyle
n_{i}(\textbf{r})=\left|\phi_{i}(\textbf{r})\right|^{2}\neq\left|\psi_{i}(\textbf{r})\right|^{2}\quad.$
(7)

The variational minimization of the SIC-DFA energy $E^{\text{SIC-DFA}}$ (Eq.
(LABEL:eq:SIC-Energy)) is more complicated because the PZ-SIC contribution
$E^{\text{PZ-SIC}}[\\{n_{i}\\}]$—the second part of Eq. (LABEL:eq:SIC-
Energy)—is not invariant under unitary transformations of the occupied KS
orbitals, $\\{\psi_{l}\\}$. Consequently, there exists an alternative
orthonormalized set of single-electron orbitals $\\{\phi_{i}\\}$, which
preserves the total density (Eq. (6)). This set can further minimize the PZ-
SIC contribution and thus the overall SIC-DFA energy. These orbitals,
$\\{\phi_{i}\\}$, referred to as SIC orbitals, are derived from unitary
transformations of the KS orbitals $\\{\psi_{l}\\}$, which are not necessarily
identical to the KS orbitals [40]. Therefore, the effort to minimize
$E^{\text{SIC-DFA}}$ (Eq. (LABEL:eq:SIC-Energy)) necessitates an additional
orthonormalization constraint on the occupied SIC orbitals $\\{\phi_{i}\\}$
[26, 25], leading to the generalized Euler equation

$\displaystyle\delta$ $\displaystyle\left[E^{\text{SIC-
DFA}}-\sum_{ij}^{N_{e}}\epsilon_{ij}\left(\langle\psi_{i}|\psi_{j}\rangle-\delta_{ij}\right)\right.$
(8)
$\displaystyle\left.-\sum_{ab}^{N_{e}}\lambda_{ab}\left(\langle\phi_{a}|\phi_{b}\rangle-\delta_{ab}\right)\right]=0\quad.$

###  II.1 PZ-SIC One-Electron Equations

Unfortunately, the fact that two different sets of orbitals, i.e.,
$\\{\phi_{i}\\}$ and $\\{\psi_{i}\\}$, enter Eq. (8) and that these orbitals
are inherently coupled complicates the solution of this equation. In practice,
the problem is approached by solving two distinct problems: one to determine
the KS orbitals $\\{\psi_{l}\\}$ and one for the SIC orbitals
$\\{\phi_{i}\\}$. To ensure self-consistency, these two sets of equations are
solved iteratively.

The PZ-SIC one-electron equations for the KS non-interacting systems are

$\displaystyle\hat{h}^{\text{SIC-DFA}}\psi_{l}=\epsilon_{l}\psi_{l}\quad.$ (9)

Here, $\hat{h}^{\text{SIC-DFA}}$ is the KS non-interacting Hamiltonian for the
SIC-DFA

$\displaystyle\hat{h}^{\text{SIC-
DFA}}=\hat{t}_{\text{s}}+\hat{v}_{\text{ext}}+\hat{v}_{\text{H}}+\hat{v}_{\text{xc}}+\hat{v}^{\text{SIC}}\quad,$
(10)

which includes the kinetic energy operator $\hat{t}_{\text{s}}$, the external
potential $\hat{v}_{\text{ext}}$, the Hartree potential $\hat{v}_{\text{H}}$,
the XC potential $\hat{v}_{\text{xc}}$, and the SIC operator
$\hat{v}^{\text{SIC}}$. Compared to the standard KS Hamiltonian, Eq. (10)
additionally features the SIC operator

$\displaystyle\hat{v}^{\text{SIC}}\psi_{l}$ $\displaystyle=\frac{\delta
E^{\text{PZ-SIC}}[\\{n_{i}\\}]}{\delta
n}\cdot\psi_{l}=\sum_{i}^{N_{e}}\hat{v}_{i}^{1e\text{SIC}}|\phi_{i}\rangle\langle\phi_{i}|\psi_{l}\rangle\quad.$
(11)

Here, the one-electron SIC potential $\\{\hat{v}^{1e\text{SIC}}_{i}\\}$ in the
SIC operator is the negative of summing the single-particle Hartree potentials
$\hat{v}_{\text{es}}[n_{i}]$ and the DFA-specific single-particle XC
potentials $\hat{v}_{\text{xc}}[n_{i}]$, evaluated for the SIC orbital
densities $n_{i}(\textbf{r})$:

$\displaystyle\hat{v}^{1e\text{SIC}}_{i}=-(\hat{v}_{\text{H}}[n_{i}]+\hat{v}_{\text{xc}}[n_{i}])$
(12)

It is important to note that the SIC operator $\hat{v}^{\text{SIC}}$ depends
on the SIC orbitals $\\{\phi_{i}\\}$. Consequently, the computations of the
SIC operator in this work occur within the framework of the generalized Kohn-
Sham (gKS) theory [57, 58, 59] in analogy to the gKS derivation of the meta-
GGA SCAN potential [6]. (Please refer to Ref. [15, 24, 14] for more detailed
derivations). Given a proper initial guess for $\\{\phi_{i}\\}$, the KS
orbitals $\\{\psi_{i}\\}$ can be obtained by solving the PZ-SIC equations (Eq.
(9)). From this, the ground-state density is derived as Eq. (6), which in turn
yields the minimal SIC-DFA total energy (Eq. (LABEL:eq:SIC-Energy))
corresponding to the initial guess $\\{\phi_{i}\\}$ [26]. In practical
scenarios, the inherent challenge of multiple solutions in the standard PZ-SIC
scheme manifests as a sensitivity of the self-consistent SIC orbitals to this
initial guess. This phenomenon is further exemplified by the cases of the
Carbon atom, Helium dimer, and Helium crystal in Sec. II.2 and Sec. II.3.

###  II.2 PZ-SIC Constraints

In the pursuit of minimizing the SIC-DFA total energy using the self-
consistent PZ-SIC scheme, it is imperative to iteratively refine the SIC
orbitals $\\{\phi_{i}\\}$ based on the KS orbitals $\\{\psi_{l}\\}$ from the
prior iteration (Eq. (9)). This updating process for the SIC orbitals must
adhere to specific physical constraints.

Constraint 1. Total density constraint (TDC):
The SIC orbitals $\\{\phi_{i}\\}$ are enforced to produce the same density
$n(\textbf{r})$ as the canonical KS orbitals $\\{\psi_{l}\\}$:

$\displaystyle\sum_{i}^{N_{e}}|\phi_{i}(\textbf{r})|^{2}=n(\textbf{r})=\sum_{l}^{N_{e}}|\psi_{l}(\textbf{r})|^{2}\quad.$
(13)

Technically speaking, this implies that the occupied SIC orbitals must be
related to the occupied KS orbitals via a unitary transformation [15]:

$\displaystyle\phi_{i}=\sum^{N_{\text{e}}}_{l}T_{il}\psi_{l}\quad,$ (14)
$\displaystyle\delta_{ml}=\sum_{i}^{N_{e}}T^{*}_{im}T_{il}\quad.$

Under this constraint and under the assumption that the KS orbitals
$\\{\psi_{l}\\}$ and the total density $n(\textbf{r})$ from the preceding
iteration remain unchanged, the generalized Euler equation (Eq. (8)) yields
the one-electron SIC potential equations:

$\displaystyle\hat{v}^{1e\text{SIC}}_{i}\phi_{i}=\sum_{a}\lambda_{ai}\phi_{a}\quad,$
(15)

which lead to the Pederson constraint that follows in Eq. (16) (Please refer
to Ref. [24] for more detailed derivations).

Constraint 2. Orbital potential constraint (OPC):

$\displaystyle\langle\phi_{m}|\hat{v}_{m}^{1e\text{SIC}}|\phi_{n}\rangle=\langle\phi_{m}|\hat{v}_{n}^{1e\text{SIC}}|\phi_{n}\rangle\quad.$
(16)

For the KS orbitals $\\{\psi_{l}\\}$ from the preceding iteration, the minimal
SIC-DFA total energy (Eq. (LABEL:eq:SIC-Energy)) is provided by the SIC
orbitals that satisfies this orbital potential constraint. The resulting SIC
orbitals are frequently (but not always) localized [25].

In summary, to find the minimal SIC-DFA energy via the variational
minimization scheme, we must solve the PZ-SIC one-electron equations self-
consistently under Constraints 1 and 2 (TDC and OPC, Eqs. (13) and (16)). The
well-documented challenge of multiple solutions arises from the fact that the
$E^{\text{SIC-DFA}}$ functional (Eq. (LABEL:eq:SIC-Energy)) possesses many
local minima. As a result, PZ-SIC outcomes are highly sensitive to the initial
guess of SIC orbitals [26].

As discussed in the introduction, the benefits of localization procedures have
been established within the framework of PZ-SIC [15]. This is a logical
outcome since more localized densities generally result in stronger self-
interaction, leading to a most significant reduction in the total energy
within the SIC-DFA scheme. Moreover, localization procedures, as exemplified
by FLO, can reduce the number of unitary transformations (Eq. (14)), thereby
facilitating the self-consistent PZ-SIC convergence. Inspired by these
localized methods, we propose the construction of SIC orbitals by self-
consistently enforcing additional strong localization. Note that different
procedures and approaches exist for this purpose [25, 23, 26]. In this work,
we choose the Edmiston-Ruedenberg (E-R) restriction [28], which, by
definition, maximizes the integrals over the squared densities of the
individual SIC orbitals

$\displaystyle\max\int
d\textbf{r}\sum_{i}n^{2}_{i}(\textbf{r})\Leftrightarrow\min\int
d\textbf{r}\sum_{i\neq j}n_{i}(\textbf{r})n_{j}(\textbf{r})\quad,$ (17)

which is equivalent to minimizing the overlap of different SIC orbital
densities [29], and thus maximizing the sum of the single-particle Hartree
energy $\sum_{i}^{N_{e}}E_{\text{H}}[n_{i}]$ [28]. For this reason, Perdew and
Zunger suggested in the original PZ-SIC paper that E-R localized orbitals are
an appropriate choice for the PZ-SIC framework as also confirmed by Pederson
and Lin [25].

In order to incorporate the E-R restriction (Eq. (17)) as a constraint into
the self-consistent PZ-SIC framework, we follow the spirit and the derivation
of the OPC to obtain the self-consistently localized SIC Euler equation

$\displaystyle\delta\left[\int
d^{3}\textbf{r}\sum_{i}n^{2}_{i}(\textbf{r})-\eta(\delta E^{\text{SIC-
DFA}})\right]=0\quad.$ (18)

This variation leads to a new constraint, see the following Eq. (19).

Constraint 3. Orbital density-potential constraint (ODPC):

$\displaystyle\langle\phi_{m}|n_{m}\hat{v}_{n}^{1e\text{SIC}}|\phi_{n}\rangle=\langle\phi_{m}|\hat{v}_{m}^{1e\text{SIC}}n_{n}|\phi_{n}\rangle\quad,$
(19)

which can be formulated as a set of localization equations (Please see
Supplementary Sec. I for a detailed derivation of these relationships). The
herein-introduced constraint based on the E-R restriction imposes an even
stronger emphasis on localization than the orbital-potential constraint. As a
consequence, it helps to steer the solution of the PZ-SIC equations towards
the correct global minimum $E^{\text{SIC-DFA}}$, especially when multiple
solutions associated with nonlocal minima are possible if only the two first
constraints are enforced. First, we explain and illustrate this for two simple
systems, He2 and C. Below, we will then further substantiate this by extended
benchmarks for solids, see Sec. II.3.

Figure 1: Depictions of the KS orbital densities for the He2 molecule in (a)
and their corresponding SIC orbital densities in (b). For the C atom, the KS
orbital densities are presented in (c), and the SIC densities in (d). The
(SIC-)SCAN calculations in FHI-aims were used to assess both sets of orbitals,
KS and SIC, employing tight tier-1 NAO basis sets and default output settings.
The results were visualized using the open-source tool Jmol [60] with its
predefined settings.

As a first example, we consider the closed-shell dihelium molecule He2 with a
bond length of 1.1 Å. Here, $1\sigma_{g}$ and $1\sigma_{u}$ are the fully
occupied KS orbitals featuring two electrons each. These occupied KS orbitals
are delocalized over the atoms, cf. Fig. 1(a). As derived in Supplementary
Sec. II, there are two sets of orbitals that fulfill both TDC and OPC. These
are the canonical KS orbitals $\\{1\sigma_{g;u}\\}$ and the more localized
orbitals $\\{\phi_{A;B}\\}$, see Fig. 1(b), that can be obtained via the
unitary transformation:

$\displaystyle\phi_{A;B}=\frac{1}{\sqrt{2}}[1\sigma_{g}\pm 1\sigma_{u}]\quad.$
(20)

More specifically, Constraint 1 (TDC, Eq. (13)) is fulfilled via

$\displaystyle n(\textbf{r})$
$\displaystyle=2|1\sigma_{g}(\textbf{r})|^{2}+2|1\sigma_{u}(\textbf{r})|^{2}$
(21)
$\displaystyle=2|\phi_{A}(\textbf{r})|^{2}+2|\phi_{B}(\textbf{r})|^{2}\quad,$

and Constraint 2 (OPC, Eq. (16)) via:

$\displaystyle\langle
1\sigma_{g}|\hat{v}^{1e\text{SIC}}_{g}|1\sigma_{u}\rangle=\langle
1\sigma_{g}|\hat{v}^{1e\text{SIC}}_{u}|1\sigma_{u}\rangle\quad,$ (22)
$\displaystyle\langle\phi_{A}|\hat{v}^{1e\text{SIC}}_{A}|\phi_{B}\rangle=\langle\phi_{A}|\hat{v}^{1e\text{SIC}}_{B}|\phi_{B}\rangle\quad.$

This means that the SIC orbitals cannot be uniquely determined via Constraints
1 and 2 (TDC and OPC), even in this relatively simple case, in which
$\\{\phi_{A;B}\\}$ are the obviously better, more localized choice.

This hurdle is overcome by Constraint 3 (ODPC, Eq. (19)), which is only
fulfilled by the set of $\\{\phi_{A;B}\\}$, but not by $\\{1\sigma_{g;u}\\}$:

$\displaystyle\langle
1\sigma_{g}|n_{g}\hat{v}^{1e\text{SIC}}_{u}|1\sigma_{u}\rangle\neq\langle
1\sigma_{g}|\hat{v}^{1e\text{SIC}}_{g}n_{u}|1\sigma_{u}\rangle\quad,$ (23)
$\displaystyle\langle\phi_{A}|n_{A}\hat{v}^{1e\text{SIC}}_{B}|\phi_{B}\rangle=\langle\phi_{A}|\hat{v}^{1e\text{SIC}}_{A}n_{B}|\phi_{B}\rangle\quad.$

As expected, the $\\{\phi_{A;B}\\}$ orbitals then also yield a SIC-SCAN total
energy that is about 5.29 eV lower than that of the $\\{1\sigma_{g;u}\\}$
orbitals 222 The SIC-SCAN calculations for He2 were carried out with FHI-aims
using the tight tier-1 NAO basis set.. This is in line with the original
intention of the constraint to enforce localization and showcases how ODPC is
a useful and often necessary complement to PZ-SIC theory.

Second, we showcase a similar but less trivial example, the SIC-SCAN total
energy calculations for the carbon atom with the electronic ground state
configuration $1s2s(2p)^{2}$. The atomic orbitals $\\{1s,2s,2p_{x},2p_{y}\\}$
shown in Fig. 1(c) satisfy both TDC and OPC [25] but not ODPC. However, all
three constraints can be simultaneously met by the following set of orbitals
obtained from a unitary transformation of the atomic orbitals:

$\displaystyle\phi_{1}=1s$ (24)
$\displaystyle\phi_{20}=\frac{1}{\sqrt{3}}2s-\frac{\sqrt{2}}{\sqrt{3}}2p_{x}$
$\displaystyle\phi_{21}=\frac{1}{\sqrt{3}}2s+\frac{1}{\sqrt{6}}2p_{x}+\frac{1}{\sqrt{2}}2p_{y}$
$\displaystyle\phi_{22}=\frac{1}{\sqrt{3}}2s+\frac{1}{\sqrt{6}}2p_{x}-\frac{1}{\sqrt{2}}2p_{y}\quad.$

As shown in Fig. 1, these hybridized orbital densities ($|\phi_{20}|^{2}$ and
$|\phi_{21}|^{2}$) break the spherical harmonic symmetry and are thus more
localized, i.e., they exhibit a smaller overlap compared to the atomic orbital
densities ($|2s|^{2}$ and $|2p_{x}|^{2}$). Accordingly, these hybridized
orbitals defined in Eq. (24) result in a SIC-SCAN energy that is about 0.86 eV
lower than that obtained using atomic orbitals.

###  II.3 Periodic Systems

In periodic boundary conditions (PBCs), the PZ-SIC one-electron equations have
to be solved for multiple k-points in the first Brillouin zone

$\displaystyle\hat{h}^{\text{SIC-
DFA}}\psi_{l,\textbf{k}}=\epsilon_{l,\textbf{k}}\psi_{l,\textbf{k}}\quad.$
(25)

The KS orbitals $\\{\psi_{l,\textbf{k}}(\textbf{r})\\}$ are the generalized
Bloch orbitals located in the first Brillouin zone. Due to the localized
nature of the SIC orbitals, Heaton, Harrison, and Lin proposed that it is more
convenient to express the SIC orbitals $\phi_{i,\textbf{I}}(\textbf{r})$ in
terms of Wannier orbitals $\\{\Psi_{l,\textbf{L}}(\textbf{r})\\}$ [22].
Accordingly, the SIC orbital $\phi_{i,\textbf{I}}(\textbf{r})$ centered in the
Ith unit cell is expressed as:

$\displaystyle\phi_{i,\textbf{I}}(\textbf{r})=\sum_{l}^{N_{e}}\sum_{\textbf{L}}^{N_{L}}T_{i\textbf{I},l\textbf{L}}\Psi_{l,\textbf{L}}(\textbf{r})\quad.$
(26)

Here, $N_{e}$ is the number of electrons in a unit cell, and $N_{L}$ is the
unit cell number (also the k-point number) in the supercell. The Wannier
orbitals $\\{\Psi_{l,\textbf{L}}(\textbf{r})\\}$ can be generated via a
Fourier transformation of the KS orbitals
$\\{\psi_{l,\textbf{k}}(\textbf{r})\\}$ and thus span the same space as the KS
orbitals,

$\displaystyle\Psi_{l,\textbf{L}}(\textbf{r})$
$\displaystyle=\Psi_{l}(\textbf{r}-\textbf{R}_{\textbf{L}})$ (27)
$\displaystyle=\frac{1}{\sqrt{N_{L}}}\sum_{\textbf{k}}^{N_{L}}\exp[-i\textbf{k}\cdot\textbf{R}_{\textbf{L}}]\psi_{l,\textbf{k}}(\textbf{r})\quad.$

Here, $\\{\textbf{R}_{\textbf{L}}\\}$ are the lattice vectors, with the
subscript $\textbf{L}=[L_{1},L_{2},L_{3}]$ being the index of a given lattice
vector.

In such a Wannier representation, TDC reads

$\displaystyle\frac{1}{N_{L}}\sum_{i\textbf{I}}^{N_{e}N_{L}}|\phi_{i,\textbf{I}}(\textbf{r})|^{2}=n(\textbf{r})=\frac{1}{N_{L}}\sum_{l,\textbf{k}}^{N_{e}N_{L}}|\psi_{l,\textbf{k}}(\textbf{r})|^{2}\quad.$
(28)

Similarly, OPC and ODPC are expressed as

$\displaystyle\langle\phi_{m,\textbf{M}}|\hat{v}_{m,\textbf{M}}^{1e\text{SIC}}|\phi_{n,\textbf{N}}\rangle=\langle\phi_{m,\textbf{M}}|\hat{v}_{n,\textbf{N}}^{1e\text{SIC}}|\phi_{n,\textbf{N}}\rangle$
(29)

and

$\displaystyle\langle\phi_{m,\textbf{M}}|n_{m,\textbf{M}}\hat{v}_{n,\textbf{N}}^{1e\text{SIC}}|\phi_{n,\textbf{N}}\rangle=\langle\phi_{m,\textbf{M}}|\hat{v}_{m,\textbf{M}}^{1e\text{SIC}}n_{n,\textbf{N}}|\phi_{n,\textbf{N}}\rangle\quad.$
(30)

In this case, $\hat{v}^{1e\text{SIC}}_{i,\textbf{I}}=\delta E^{\text{PZ-
SIC}}[\\{n_{i,\textbf{I}}\\}]/\delta n_{i,\textbf{I}}$ is the one-electron SIC
potential under PBCs, and $n_{i,\textbf{I}}=|\phi_{i,\textbf{I}}|^{2}$ is the
associated SIC orbital density (Please refer to Ref. [22, 44] and
Supplementary Sec. I for more detailed derivations).

Figure 2:  (a) Schematic representation of the hcp Helium crystal’s unit cell
on the left; while the depiction of two unit cells incorporating periodic
boundary conditions on the right. Geometry was taken from Ref. [62]. (b) The
energy profile of $E^{\text{PZ-SIC}}$ (Eq. (2)) with respect to the rotation
angles ($\theta_{1},\theta_{2}$). (c) The absolute errors (AE) defined in Eq.
(31) pertaining to OPC (Constraint 2) and (d) ODPC (Constraint 3),
respectively. (e) Combined AE of both OPC and ODPC for all possible orbitals
in the $1\times 1\times 2$ Helium unit cell. These orbitals are unitarily
transformed from the $1\Sigma_{g,\textbf{I}},1\Sigma_{u,\textbf{I}}$ orbitals
using rotation angles $\theta_{1}$ and $\theta_{2}$. In the subfigures (c-e),
the regions in white indicate the solutions for OPC, ODPC, and their
combinations, respectively. Both SCAN and SIC-SCAN calculations in FHI-aims
employed PBCs and tight tier-1 NAO basis sets.

As mentioned in the introduction, the SIC multiple-solution problem is known
to be more serious for solids [22, 42, 43, 63, 41]. To illustrate this and to
showcase how the introduction of ODPC helps to alleviate this issue, we here
discuss crystal He, a prototypical case with a multiple-solution problem in
standard SIC approaches. As there are two He atoms in each unit cell, see Fig.
2(a), the occupied KS orbitals in the Wannier representation can be marked as
$1\Sigma_{g,\textbf{I}}$ and $1\Sigma_{u,\textbf{I}}$. Naturally, these
orbitals have the same shapes as the occupied molecular orbitals in the
closed-shell dihelium He2 molecule, see Fig. 1-(a). All orbitals
$\\{(\Phi_{1,\textbf{I}},\Phi_{2,\textbf{I}})\\}$ that satisfy TDC can be
expressed via unitary transformations of the Wannier KS orbitals
($1\Sigma_{g,\textbf{I}},1\Sigma_{u,\textbf{I}}$). In this particular simple
case, two distinct angles $\theta_{I=1,2}\in[0,\pi/2]$ can be used to
characterize this unitary transformation, as derived in detail in
Supplementary Sec. III. Accordingly, the original Wannier KS orbitals
($1\Sigma_{g,\textbf{I}},1\Sigma_{u,\textbf{I}}$) are obtained by choosing the
angles to be $(\theta_{1}=0,\theta_{2}=\pi/4)$.

To visualize the influence of the actual choice of
$\\{(\Phi_{1,\textbf{I}},\Phi_{2,\textbf{I}})\\}$, we here show $E^{\text{PZ-
SIC}}$ as an energy surface as function of $\theta_{I=1,2}$ in Fig. 2-(b). The
global $E^{\text{PZ-SIC}}$ minimum of 0.32 eV is obtained for the SIC orbitals
$\Phi_{1;2,\textbf{I}}$ with $(\theta_{1}=\pi/4,\theta_{2}=\pi/2)$. Respective
SIC orbitals at these angles give
$\frac{1}{\sqrt{2}}[1\Sigma_{g,\textbf{I}}\pm 1\Sigma_{u,\textbf{I}}]$, which,
not too surprisingly, localize at He atoms and have similar shapes as in the
case of the non-periodic, closed-shell He2 molecule, see Fig. 1-(b). At
variance with Fig. 2-(b), which is produced using all three constraints, Fig.
2-(c) and Fig. 2-(d) show what happens if only TDC and OPC viz. TDC and ODPC
are used, respectively. Additionally, Fig. 2-(e) shows, once more, the results
of using all three constraints. To better visualize the differences in these
cases, we here locate the solutions of constraints in the individual plots by
using the absolute error (AE) of OPC or ODPC:

$\displaystyle\text{AE}=\left\\{\begin{aligned}
&\text{max}\\{|\langle\phi_{m,\textbf{M}}|\hat{v}_{m,\textbf{M}}^{1e\text{SIC}}-\hat{v}_{n,\textbf{N}}^{1e\text{SIC}}|\phi_{n,\textbf{N}}\rangle|\\}\\\
&\text{max}\\{\\\
&~{}~{}~{}|\langle\phi_{m,\textbf{M}}|n_{m,\textbf{M}}\hat{v}_{n,\textbf{N}}^{1e\text{SIC}}-\hat{v}_{m,\textbf{M}}^{1e\text{SIC}}n_{n,\textbf{N}}|\phi_{n,\textbf{N}}\rangle|\\}\\\
\end{aligned}\right.\quad.$ (31)

This reveals that the standard SIC method, which only accounts for TDC and
OPC, is multi-valued since the three solutions at (0, $\pi/4$), ($\pi/2$,
$\pi/4$), and ($\pi/4$, $\pi/2$) satisfy TDC and OPC, as shown in Fig. 2-(c).
Similarly, using only TDC and ODPC also yields multi-valued solutions at
($\pi/4$, 0.0) and ($\pi/4$, $\pi/2$), as shown in Fig. 2-(d). The solution
becomes uniquely single-valued only when all three constraints are employed,
as shown in Fig. 2-(e). In this case, only one solution, i.e., ($\pi/4$,
$\pi/2$), fulfills all constraints. This position also corresponds to the
minimal SIC energy shown in Fig. 2-(a).

Figure 3:  (a) Primitive unit cell of the periodic polyethylene chain,
$($C2H${}_{4})_{n}$. (b) Schematic view of the periodic polyethylene chain
with two unit cells along the chain direction. Geometry was taken from Ref.
[64]. (c) PZ-SIC energies based on the SCAN functional (Eq. (2)). Different
initial guesses of SIC orbitals were generated with respect to the mixing
factor $a$ defined in Eq. (32). All calculations were performed using tight
tier-1 NAO basis sets in FHI-aims.

As a last example, we turn our attention to the polyethylene chain,
$($C2H${}_{4})_{n}$. Its geometry is depicted in Fig. 3-(a) and (b). In such a
complex system, the SIC orbitals cannot be exhaustively explored using a low-
dimensional parametric model, unlike the simpler model systems previously
discussed. Nonetheless, we here investigate the dependence of the SIC
procedure by employing a set of initial guesses for the SIC orbitals, which
are grounded in physical reasoning. These initial guesses are formulated by
mixing the atomic orbitals $\varphi_{i}$ centered at the atomic position
$\textbf{R}^{ac}_{i}$ and the Wannier orbitals
$\Psi^{\text{SCAN}}_{i,\textbf{L}}$ obtained from SCAN calculations without
including SIC:

$\displaystyle\phi^{\text{init}}_{i,\textbf{I}}=$
$\displaystyle\frac{1}{A_{i}}\left[\right.a\varphi_{i}(\textbf{r}-\textbf{R}^{ac}_{i}-\textbf{R}_{\textbf{I}})$
(32)
$\displaystyle+(1-a)\Psi^{\text{SCAN}}_{i}(\textbf{r}-\textbf{R}_{\textbf{I}})\left.\right]\quad.$

Here, $A_{i}$ denotes a normalization factor, and $a$ is a mixing factor that
ranges from $0$ to $1$. Accordingly, turning up $a$ from 0 to 1, the initial
guess $\phi^{\text{init}}_{i,\textbf{I}}$ can be varied from having the pure
atomic-orbital character to being a Wannier orbital, i.e., a KS molecular
orbitals associated with the unit cell I in real space. This allows tuning the
localization of the initial guess, since, as discussed in the case of He2, the
atom orbitals that surround the atomic centers can be expected to be more
localized than the Wannier orbitals [63]. As shown in Fig. 3-(c), different
solutions are found for different values of $a$, with energies varying by
about 0.3 eV at most, if only TDC and OPC are taken into account. This
indicates that the problem has multiple local minima. Let us emphasize that
even when starting from fully localized orbitals $a=0$, TDC and OPC are not
converging to the lowest energy solution. However, this doesn’t apply when
ODPC is additionally imposed because, in this case, all initial guesses
converge to the same lowest-energy solution.

##  III Results

We have implemented the aforementioned PZ-SIC method within the FHI-aims
package [54]. Detailed aspects of the implementation can be found in
Supplementary Sec. IV. We here present benchmark results for a range of finite
and periodic systems, including fundamental atoms, molecules, and bulk solids.
It is worth noting that there are two distinct approaches for representing SIC
orbitals in finite systems [65]. One can either use a real-valued unitary
transformation matrix $T$ or a complex one (Eq. (14)). While both real and
complex transformation matrices yield identical total densities, the resulting
SIC energies can vary [66]. This section will present results for both real
and complex SIC orbitals when dealing with finite systems. For periodic
systems, where KS orbitals are inherently complex, we will exclusively
consider complex SIC orbitals.

###  III.1 Finite Systems

####  III.1.1 Ionization Potential

The ionization potential (IP) measures the capability of an element to
participate in chemical reactions that necessitate ion formation or electron
transfer. IP is typically calculated via the total energy difference between
the neutral molecule and the corresponding ion. Within the framework of (g)KS
approach, as applied in specific DFAs, the IP is theoretically equivalent to
the negative of the energy of the highest occupied molecular orbital (HOMO),
as delineated in Ref. [67, 68]. However, this relationship often encounters
deviations due to the errors, most notably the many-electron SIE. While gKS
outperforms KS in orbital energy calculations [69], these errors still lead to
significant discrepancies in the IP and predicted HOMO energies using (g)KS
approaches [68]. In other words, this inconsistency affects the accuracy of IP
estimations derived from HOMO energies. Consequently, the precision of
calculated IPs based on HOMO energies remains a crucial metric for evaluating
the effectiveness of various electronic-structure theory methodologies,
including SIC [56].

The HOMO energies for atoms ranging from H to Ar, as calculated by LDA, PBE,
SCAN, and SIC-SCAN, are presented in Fig. 4, with experimental values [70]
provided for reference. The pronounced SIE in LDA and PBE [71] leads to a
notable underestimation of the IPs. The performance of SCAN offers a marginal
improvement. The mean absolute percentage errors (MAPE) are 42.03% for LDA,
41.75% for PBE, and 38.5% for SCAN, substantiating that SCAN also suffers
under SIEs. However, with SIC-SCAN, there is a notable improvement in the
predictions. As shown in Fig. 4-top, the MAPE is reduced to 5.02% when using
real SIC orbitals and further drops to 3.95% for complex SIC orbitals. This is
in line with the findings in Ref. [65], where the use of complex SIC orbitals
yielded superior results compared to their real counterparts.

Figure 4:  Performance of various DFAs in describing the ionization potentials
for the atoms ranging from H to Ar. SIC-SCAN calculations using real and
complex SIC orbitals are marked as SIC-SCAN(real) and SIC-SCAN(complex),
respectively. (top) Relative errors of $-\epsilon_{\text{HOMO}}$ energies to
IPs, while the $-\epsilon_{\text{HOMO}}$ values are given in (bottom). All
calculations were performed using tight tier-1 NAO basis sets. The
experimental IPs were taken from the NIST database [70].

Fig. 5 showcases the performance of various DFAs in predicting negative HOMO
energies ($-\epsilon_{\text{HOMO}}$) for a set of 18 molecules. As observed,
LDA, PBE, and SCAN calculations tend to underestimate the
$-\epsilon_{\text{HOMO}}$ energies significantly, whereas the SIC approach
considerably improves the results. Specifically, all absolute percentage
errors of SIC-SCAN(complex) are under 23%, and its MAPE is $12.8\%$, which is
substantially lower than the SCAN error ($33.5\%$). We note, however, that the
PZ-SIC method tends to slightly over-correct the SCAN calculations, resulting
in the $-\epsilon_{\text{HOMO}}$ energies consistently higher than the
experimental IPs. To contextualize these results, $G_{0}W_{0}$@PBE, as
reported in Ref. [56], offers IP predictions with a mere $3\%$ MAPE.

Figure 5:  (top) Performance of various DFAs in describing the ionization
potentials of 18 molecules. The mean absolute percentage errors (MAPEs) of
HOMO energies are given in (bottom). Calculations of LDA, PBE, SCAN, and SIC-
SCAN methods were carried out in FHI-aims using tight tier-1 NAO basis sets.
$G_{0}W_{0}$@PBE results were taken from Ref. [56], which were also calculated
in FHI-aims but utilized def2-QZVP basis sets. The experimental vertical
ionization energies, taken from the NIST database [70], are reported for
comparison.

####  III.1.2 Energy Curve and Broken Symmetry

As the SIE tends to delocalize electrons, semi-local DFAs typically struggle
to describe charge transfer processes. A typical example is that semi-local
DFAs tend to overly preserve the ground state symmetry of dissociating neutral
heterodimers by producing fractionally charged fragments. For instance, SCAN
incorrectly predicts a positive fractional charge on the H atom in the
dissociation of the molecule H-F (Fig. 6), hence giving too high energies in
the dissociation limit. While Kim et al. [72] highlighted that a meticulous
selection of orbital occupations can address this issue, the hybrid functional
PBE0 [73] outperforms both PBE and SCAN with self-consistent calculations.
This superiority is attributed to incorporating a fraction (0.25) of Fock
exchange in PBE0, which reduces the SIE in PBE. Nevertheless, PBE0 cannot
correctly predict the change distribution at distances between 1.5-2.5 Å, as
illustrated in Fig. 6-bottom. In contrast, SIC reproduces comparable charge
analyses of references and successfully maintains electron-neutral atoms in
the dissociation limit. Here, accurate results from coupled cluster theory
with single, double, and perturbative triple excitations (CCSD(T)) are used as
references. As shown, the energy curve of SIC-SCAN(complex) matches the
CCSD(T) curve in both equilibrium (around $r\approx 0.8$ Å) and stretched-bond
($r>0.8$ Å) regions. For comparison, the energy curve of SIC-SCAN(real)
deviates from the one of CCSD(T) in the equilibrium region.

Figure 6:  (top) HF dissociation energy curves of various methods with the
zero-energy level set to the total energy of isolated atoms/ions. (bottom)
Mulliken charge analyses on the H atom along the dissociation. All
calculations were performed in FHI-aims. Tight tier-1 NAO basis sets were
employed for PBE, PBE0, SCAN and SIC-SCAN incorporating real and complex SIC
orbitals, while the CCSD(T) results utilize cc-pVTZ Gaussian-type basis sets.
Figure 7:  Dissociation curves of H${}_{2}^{+}$ (a), He${}_{2}^{+}$ (b), and
Ne${}_{2}^{+}$ (c), with the zero-energy reference set to the total energy of
isolated atoms/ions. The mirror symmetry ensures no charge transfer in these
two dissociations, as demonstrated by the Mulliken charge analyses. All
calculations were performed in FHI-aims. Tight tier-1 NAO basis sets were
employed for SCAN and SIC-SCAN incorporating real and complex SIC orbitals,
while the CCSD(T) results utilize cc-pVTZ Gaussian-type basis sets. The
Hartree-Fock (HF) method, being exact for one-electron systems, provided the
H${}_{2}^{+}$ dissociation reference curve.

Fig. 7 illustrates how the SIE in SCAN calculations affects open-shell
configurations, e.g., H${}^{+}_{2}$, He${}^{+}_{2}$ and Ne${}^{+}_{2}$. For
H${}^{+}_{2}$, which contains only one electron, SIC exactly offsets the SIE
in the SCAN calculation. However, for the He${}^{+}_{2}$ and Ne${}^{+}_{2}$
dissociations with more electrons, this exact cancellation does not occur. As
the dissociation behaviors of He${}^{+}_{2}$ and Ne${}^{+}_{2}$ are similar,
we take the Ne${}^{+}_{2}$ system with more electrons as an example to analyze
in the following. The conclusions drawn for Ne${}^{+}_{2}$ are applicable to
He${}^{+}_{2}$ as well. Nevertheless, SIC, whether using real or complex SIC
orbitals, significantly improves the SCAN curve. The resulting SIC-SCAN curves
align closely with CCSD(T) calculations, indicating that the PZ-SIC scheme
successfully rectifies the majority of the many-electron SIEs in the SCAN
calculations of the Ne${}^{+}_{2}$ dissociation. As shown in Fig. 7, both SCAN
and SIC-SCAN, without symmetry breaking, predict the same charge distribution,
with each Ne atom carrying a fractional charge of 0.5 (Ne0.5). It is also
worth noting the combined total energy of Ne and Ne+1, which is taken as the
reference for the Ne${}^{+}_{2}$ dissociation curve in Fig. 7. Regardless of
the chosen SIC orbitals (either real or complex), SIC-SCAN closely reproduces
the degeneracy between the configurations of Ne${}^{+}\cdots$Ne and
Ne${}^{0.5}\cdots$Ne0.5, with a deviation of less than 0.3 eV. This
observation suggests that SIC-SCAN can provide an accurate description of
fractional charged systems, which has been widely used to discuss the SIE in
the conventional DFAs [74].

###  III.2 Periodic Systems

####  III.2.1 Performance for Crystal Properties

Given that our SIC-SCAN implementation effectively addresses the multiple
solution issue in both finite and period systems, we have extended our
attention to the crystal properties, including cohesive energy, bulk modulus,
and lattice constant.

Figure 8:  (a) Lattice constants at equilibrium (Å), bulk modulus (GPa) and
cohesive energies (eV) of the diamond Si, diamond C, rocksalt LiH, fcc Ne, and
zincblende AlP, calculated using SCAN and SIC-SCAN calculations. (b) Absolute
errors in the calculated lattice constants (Å), bulk modulus (GPa), and
cohesive energies (eV) for SCAN and SIC-SCAN. All calculations were carried
out in FHI-aims using tight tier-1 NAO basis sets and the k-grid setting of
$4\times 4\times 4$. The experimental reference values were collected from the
main-group test set [75, 76].

We evaluate the cohesive energies of a fundamental set of materials, including
diamond Si, diamond C, rocksalt LiH, fcc Ne, and zincblende AlP, as presented
in Fig. 8-(a). Additionally, we computed the lattice constant and bulk modulus
for this test set, with results also displayed in Fig. 8-(a). The relevant
absolute errors (AEs) are plotted in Fig. 8-(b) for comparative analysis.

Consistent with findings in previous studies [77, 52], our calculations
suggest that the SIC contribution tends to deteriorate the accuracy of SCAN
predictions for cohesive energies and bulk modulus. However, the impact of SIC
on SCAN’s lattice constant predictions is relatively minor. As a typical
system of van der Waals (vdW) interaction, the fcc Ne poses a challenge for
standard DFAs. Our results suggest that the SCAN method tends to overestimate
the vdW interaction, with its calculated cohesive energy being double the
experimental value and its lattice constant being 0.37 Å shorter. Given that
the PZ-SIC scheme is irrelevant to addressing the SCAN error for vdW
interactions, SIC-SCAN does not improve the terrible performance of SCAN in
predicting the lattice constant of fcc Ne. The observed improvements in the
cohesive energy and bulk modulus of the fcc Ne might be due to a fortuitous
offset between the inappropriate SIC contribution and SCAN’s overestimation of
vdW interactions.

####  III.2.2 Band Gaps

The fundamental band gap is a ground-state property, which is defined as the
difference between the ionization energy and the electron affinity [67]. For
periodic solids, it has been demonstrated that, for any DFA within the KS or
gKS framework, the fundamental band gap is equivalent to the energy difference
between the lowest unoccupied orbital (LUMO) and HOMO, providing that the
given DFA potential is continuous upon introducing an electron or hole [67,
69]. This finding is consistent across both the KS and gKS methods. The well-
documented underestimation of the band gap in semi-local DFAs, including LDA,
PBE, and SCAN, is largely attributed to significant SIEs present in these
approximations [2, 67]. In this work, we evaluate the performance of various
methods in predicting band gaps for a selection of 10 fundamental solids. This
includes covalent crystals (such as Si, C), ionic crystals (such as NaCl,
LiF), and Mott insulators composed of transition-metal oxides (MnO). The
experimental band gaps are taken as the reference. It’s worth emphasizing that
although electron-phonon interactions can influence experimental band gaps
[78, 79, 80], the remarkable agreement between the band gaps determined by the
$G_{0}W_{0}$ method and experimental data suggests that, at least for the
solids investigated here, the impact of this interaction is probably less
significant compared to the SIE in DFAs. As illustrated in Fig. 9, the mean
absolute percentage errors (MAPEs) of LDA and PBE are 45.4 % and 43.2 %,
respectively, substantially higher than the $G_{0}W_{0}$ value (4.7 %). SCAN
outperforms LDA and PBE with a lower MAPE (31.5 %), as the step outside KS to
gKS is to enhance orbital energies. Unfortunately, However, Figure 9 reveals
that SCAN still has poor accuracy in predicting band gaps, consistent with
previous findings [17]. Our results suggest that applying SIC to SCAN can
reduce the band gap error, lowering the MAPE to 18.5 %. We acknowledge that
SIE may not be the sole source of error for SCAN, but it significantly
contributes to inaccuracies in band gap calculations.

MnO is classified as a Mott insulator characterized by strong electron
interactions [81, 82]. Accurately predicting the band gap of MnO is a
challenge for conventional semi-local DFAs [83]. Our results confirm that LDA,
PBE, and SCAN methods predict the band gaps for MnO around 0.60 eV, a
significant underestimation compared to the experimental value of 3.70 eV. A
common remedy is the DFT$+U$ method, which introduces an empirical electron
repulsion “$U$” [84, 85]. For instance, GGA$+U$ gives a band gap of
approximately 4 eV when applying a “$U$” of 5.25 eV via the linear-response
approach [86]. Moreover, the self-consistent $GW_{0}$ approach, based on the
LDA$+U$ wave functions, produces a band gap of 3.32 eV [87]. In contrast, our
research indicates that the PZ-SIC scheme can also effectively improve the
SCAN description for MnO’s band gap. The SIC-SCAN method calculates a band gap
of 3.73 eV, which is impressively close to the experimental value, deviating
by a mere 0.03 eV. This result aligns with previous studies, suggesting that
PZ-SIC captures the strong electron interactions often overlooked in standard
semi-local DFAs [15, 83].

Figure 9:  Calculated band gaps using LDA, PBE, SCAN, (sd-)SIC-SCAN, and $GW$
against experimental counterparters. The graph features a black diagonal line
representing the perfect agreement to the experiment. Mean absolute percentage
errors (MAPEs) are collected below the graph. Calculations of LDA, PBE, SCAN,
and (sd-)SIC-SCAN were performed in FHI-aims. A k-grid setting of $5\times
5\times 5$ was used together with tight tier-1 NAO basis sets. Geometries were
taken from the experiments. Results for $G_{0}W_{0}$@PBE and experimental data
were collected from Refs. [88, 89].

To further improve the calculated band gaps in the KS-DFT framework, a scaled-
down PZ-SIC (sd-SIC) variant has been proposed [90, 44, 91]. Compared to the
standard PZ-SIC scheme (Eq. 2), the pivotal adjustment in sd-SIC is to
introduce a set of orbital-density specific scaling factors, denoted as
$\\{X_{i}\\}$, into the self-interaction correction

$\displaystyle E^{\text{sd-SIC-DFA}}=E^{\text{DFA}}+E^{\text{sd-SIC}}$ (33)
$\displaystyle E^{\text{sd-
SIC}}[\\{n_{i}\\}]=-\sum_{i}^{N_{\text{e}}}X_{i}\left(E_{\text{xc}}^{\text{DFA}}[n_{i}]+E_{\text{H}}[n_{i}]\right)\quad.$

In addition to determining the scaling factor for each orbital, adopting a
uniform scaling factor across all orbital densities ($X_{i}=X$) presents a
straightforward and effective strategy, albeit it takes on empirical impacts.
For instance, it has been observed that a scaling factor of $0.4$ enhances
thermochemistry [90], while a factor of $0.2$ yields desirable band gaps in
both silicon and diamond [92]. Furthermore, the coeﬀicient $X=2/3$, as
proposed by Shinde et al. [44] in 2020, to enhanced the SIC-PBE results for
fundamental solids. Building upon this, our study investigates the impact of
this coefficient on the SIC-SCAN band gaps. We discovered that a scaling
factor of 0.60 yields optimal results for SIC-SCAN. These findings are
comprehensively discussed in Supplementary Section V. The outcomes are
illustrated in Fig. 9, where they are denoted as scaled-down SIC-SCAN (sd-SIC-
SCAN). In agreement with the previous research by Shinde et al. [44], we find
that the scaled-down algorithm with a constant scaling factor enhances the
accuracy of the calculated band gaps. The sd-SIC-SCAN approach outperforms the
SIC-SCAN method, bringing down the MAPE to 12.9%. However, it is crucial to
highlight that the standard PZ-SIC method is exact for one-electron systems,
which, unfortunately, is not the case in the sd-SIC method.

####  III.2.3 Crystalline trans-Polyacetylene

Figure 10:  (a) Experimental unit cell of the TPA structure [93]. Carbon atoms
are depicted in gray, while hydrogen atoms are in white. (b) and (c) are the
right and bottom side views, respectively. (d) The schematic view of a single
(CH)x chain with the dimerization parameter $\Delta z$, representing the
difference between the distance of single-bond and double-bond carbon atoms
along the chain direction $z_{s}$ and $z_{d}$. (e) Total energy curves as a
function of the dimerization parameter $\Delta z$, computed using the SCAN,
(sd-)SIC-SCAN, PBE0, and RPA@PBE methods. All calculations were executed in
FHI-aims. Vertical red dashed lines in the subfigure (e) mark the experimental
dimerization parameter ($\Delta z\approx 0.047$ Å) [94]. The calculated total
energies for $\Delta z=0$ were set as the zero-energy reference. Markers
pinpointed the optimized dimerization for each method. A dense k-grid setting
of $14\times 12\times 16$ was used together with tight tier-1 NAO basis sets
for these calculations. As the dimerization factor varies, the combined length
of $z_{s}+z_{d}$ remained equal to the lattice vector length $|c|$. A vertical
height $h$ of $0.70$ Å–from the middle carbon atom to the connection of the
other two carbon atoms–was maintained.

We then turn to a challenging problem for periodic systems, the broken
symmetry in crystalline polyacetylene (CH)x. The polyacetylene chain has two
configurations, named cis and trans. In 1957, Ooshika found that the trans
configuration is a stable configuration for the long-chain polyacetylene (TPA)
[95, 96, 97], where a sequence of alternating single and double carbon bonds
is formed along the chain, as illustrated in Fig. 10-(a-d). The dimerization
parameter $\Delta z$ describes the length difference of these bonds along the
chain’s direction ($z_{s}$ and $z_{d}$):

$\displaystyle\Delta z=z_{s}-z_{d}\quad.$ (34)

Here, $z_{s}$ and $z_{d}$ are vertical distances from the central carbon atom
on both sides. X-ray studies on crystalline TPA identified a dimerization of
$\Delta z\approx 0.047$ Å [94]. The underlying mechanism of symmetry break is
yet to be clarified. On the one hand, the formation of alternating bonds and
charge transfer between carbons has been explained by using the celebrated 1D
Su-Schrieffer-Heeger model, emphasizing the phonon/electron-phonon interaction
[98]. On the other hand, the electron-electron interaction could be another
possible driving force for the induced dimerization [99, 100, 73].

Fig. 10-(e) displays the TPA energy curves as a function of the dimerization
factor $\Delta z$, which are calculated using the SCAN, SIC-SCAN, PBE0, and
RPA@PBE methods. In this context, “RPA@PBE” denotes the random phase
approximation based on PBE density and orbitals[101]. The standard hybrid
functional PBE0 [73], which incorporates a portion of the Hartree-Fock
exchange with $\alpha=0.25$, predicts a dimerization of 0.044 Å. This aligns
closely with the experimental findings (0.047 Å) and the RPA@PBE calculation
(0.045Å) [101]. Our analysis indicates that the calculated dimerization is
sensitive to the $\alpha$ value. When $\alpha$ is increased to 0.50, the
resulting dimerization parameter rises to 0.059 Å, overshooting the
experimental value by 0.012 Å.

While all the methods investigated here predict a symmetry-broken trans
configuration with $\Delta z>0$, the dimerization factor $\Delta z$ determined
by the SCAN method is 0.026 Å, significantly shorter than the experimental
measurement of 0.047 Å. In agreement with a previous study [102], our results
confirm that the PZ-SIC scheme enhances the accuracy of conventional semi-
local DFAs for this system. SIC-SCAN predicts a more extended dimerization
parameter of around 0.036 Å. While the calculated dimerization parameter
remains unchanged when applying a uniform scaling factor $X=0.60$ via the
scaled-down SIC formula (Eq. (33)), the sd-SIC-SCAN method considerably
diminishes the energy well compared to the SIC-SCAN method.

Figure 11: Calculated band gaps of crystalline trans-polyacetylene using the
SCAN, PBE0 and SIC-SCAN methods. All calculations were performed in FHI-aims.
The experimental geometry was taken from Ref. [94]. A dense k-grid setting of
$14\times 12\times 16$ was used together with the tight tier-1 NAO basis sets.
The blue dash line and orange line presented the direct and indirect
experimental band gaps, respectively, which were taken from Ref. [94].

Finally, we turn our attention to the band gap of crystalline TPA, using its
experimental geometry from reference [94]. Direct and indirect adsorption gaps
of TPA have been identified by Fincheret al., which are 1.4 eV and 1.1 eV,
respectively [94]. Fig. 11 presents the direct and indirect band gaps as
calculated by various methods. The standard hybrid PBE0 with $\alpha=0.25$
slightly underestimates the indirect band gap at 0.95 eV. In contrast, using
$\alpha=0.50$ with PBE0 leads to a notable overestimation. As expected, due to
the serious SIE, SCAN significantly underestimates the band gap, erroneously
predicting identical values of 0.12 eV for both direct and indirect gaps. The
SIC-SCAN method, on the other hand, rectifies the difference between direct
and indirect gaps but greatly overestimates the indirect band gap to 4.50 eV.
Introducing a uniform scaling factor, $X$, alleviates the overcorrection of
PZ-SIC, but the calculated indirect band gap of sd-SIC-SCAN remains
considerably larger than the experimental findings.

Crystalline TPA represents a challenging periodic system to study. Numerous
factors, including lattice constants [103], can influence its calculated
properties. As illustrated in Fig. 10 and Fig. 11, significant discrepancies
are evident in energy and band gap values across various computational
methods. The widely used hybrid functional PBE0, with a default alpha of
$0.25$, predicts dimerization that closely aligns with the results of RPA@PBE.
However, increasing the alpha value to 0.5 adversely affects PBE0’s
performance. It’s worth noting that PBE0 outperforms the (SIC-)SCAN method in
predicting these TPA properties, highlighting that while mitigating the SIE is
crucial, it is not the only factor to consider when analyzing crystalline TPA
using density functional theory.

##  IV Conclusion

The state-of-the-art, non-empirical meta-GGA SCAN functional provides a
significant and consistent improvement over the non-empirical GGA PBE across a
range of chemical interactions in both molecules and solids. However, it is
not immune to the pervasive self-interaction error. In this study, we
implement and generalize the self-interaction correction algorithm proposed by
Perdew and Zunger (PZ-SIC), to alleviate this error within the SCAN
functional. Specifically, we introduce the orbital density-potential
constraint, which, in conjunction with the total density constraint and the
orbital potential constraint, facilitates the self-consistent localization of
the SIC orbitals, leading to an accurate PZ-SIC energy. Our findings indicate
that the orbital density-potential constraint can narrow down multiple
potential SIC orbitals to a singular optimal choice in systems like C, He2,
hcp-Helium, and the CH-chain.

Our PZ-SIC approach markedly enhances the performance of SCAN, as evidenced by
the highest occupied molecular orbital (HOMO) energies
$\epsilon_{\text{HOMO}}$ of H-Ar atoms and a set of 18 molecules. Built upon
the accurate description of orbital energies, our PZ-SIC approach yields
commendable energy curves and corrected charge transfers in the dissociations
of H-F and Ne-Ne+. When tested on 10 foundational solids, SIC-SCAN
consistently outperforms SCAN in band gap predictions. However, for other
crystalline properties, such as cohesive energy, bulk modulus, and lattice
constants, SIC-SCAN does not fare as well as SCAN for our test set, which
includes diamond Si, diamond C, rocksalt LiH, fcc Ne, and zincblende AlP. We
then delve into the ground state of crystalline polyacetylene, a renowned
challenge for periodic systems that remains elusive for semi-local DFT
calculations like the SCAN functional. Our findings indicate that the PZ-SIC
methodology adeptly identifies the symmetry-breaking in crystalline
polyacetylene, emphasizing the role of electron-electron interactions.

## supplementary Information

The following files are available free of charge:

  * •

supplementary.pdf, containing additional information needed for detailed
discussions.

  * •

data.zip, containing geometries and electronic-structure energy results
produced in this project and also available on the NOMAD repository with the
DOI: 10.17172/NOMAD/2023.11.28-1 .

###### Acknowledgements.

Sheng thanks Zhenkun Yuan, Hagen-Henrik Kowalski and Florian Knoop for
inspiring discussions. This work received funding from the European Union’s
Horizon 2020 Research and Innovation Programme (grant agreement No. 951786,
the NOMAD CoE), and the ERC Advanced Grant TEC1P (No. 740233). This work was
partially supported by the National Natural Science Foundation of China (No.
21973015, No. 22125301) and Innovative Research Team of High-Level Local
universities in Shanghai, and a Key Laboratory Program of the Education
Commission of Shanghai Municipality (ZDSYS14005).

## References

  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Draxl and Scheffler [2020] C. Draxl and M. Scheffler, “Big data-driven materials science and its fair data infrastructure,” in  _Handbook of Materials Modeling_ (Springer, Cham, 2020) pp. 49–73.
  * Carbogno  _et al._ [2022] C. Carbogno, K. S. Thygesen, B. Bieniek, C. Draxl, L. M. Ghiringhelli, A. Gulans, O. T. Hofmann, K. W. Jacobsen, S. Lubeck, J. J. Mortensen, M. Strange, E. Wruss, and M. Scheffler, NPJ Comput. Mater. 8 (2022).
  * Teale and et al. [2022] A. Teale and et al., ChemRxiv: 2022-13j2v (2022).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Perdew [1985] J. P. Perdew, Phys. Rev. Lett. 55, 1665 (1985).
  * Ghosh and Parr [1986] S. K. Ghosh and R. G. Parr, Phys. Rev. A 34, 785 (1986).
  * Perdew  _et al._ [1999] J. P. Perdew, S. Kurth, A. Zupan, and P. Blaha, Phys. Rev. Lett. 82, 2544 (1999).
  * Zhang  _et al._ [2017] Y. Zhang, J. Sun, J. P. Perdew, and X. Wu, ArXiv: 1702.04817 (2017).
  * Peng  _et al._ [2016] H. Peng, Z. Yang, J. P. Perdew, and J. Sun, Phys. Rev. X 6, 041005 (2016).
  * Sun  _et al._ [2016] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, M. L. Klein, and J. P. Perdew, Nat. Chem. 8, 831 (2016).
  * Mori-Sánchez  _et al._ [2006] P. Mori-Sánchez, A. J. Cohen, and W. Yang, J. Chem. Phys. 125, 201102 (2006).
  * Perdew [1990] J. P. Perdew, Adv. Quantum Chem. 21, 113 (1990).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * Cohen  _et al._ [2012] A. J. Cohen, P. Mori-Sánchez, and W. Yang, Chem. Rev. 112, 289 (2012).
  * Isaacs and Wolverton [2018] E. B. Isaacs and C. Wolverton, Phys. Rev. Mater. 2, 063801 (2018).
  * Note [1] Hartree atomic units are used throughout the paper.
  * Cohen  _et al._ [2008] A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008).
  * Mori-Sánchez  _et al._ [2008] P. Mori-Sánchez, A. J. Cohen, and W. Yang, Phys. Rev. Lett. 100, 146401 (2008).
  * Vydrov  _et al._ [2007] O. A. Vydrov, G. E. Scuseria, and J. P. Perdew, J. Chem. Phys. 126, 154109 (2007).
  * Heaton  _et al._ [1983] R. A. Heaton, J. G. Harrison, and C. C. Lin, Phys. Rev. B 28, 5992 (1983).
  * Bylaska  _et al._ [2004] E. Bylaska, K. Tsemekhman, and H. Jonsson, in  _APS March Meeting Abstracts_, Vol. 2004 (2004) p. L38.004.
  * Pederson  _et al._ [1984] M. R. Pederson, R. A. Heaton, and C. C. Lin, J. Chem. Phys. 80, 1972 (1984).
  * Pederson and Lin [1988] M. R. Pederson and C. C. Lin, J. Chem. Phys. 88, 1807 (1988).
  * Lehtola  _et al._ [2016] S. Lehtola, M. Head-Gordon, and H. Jónsson, J. Chem. Theory Comput. 12, 3195 (2016).
  * Pipek and Mezey [1989] J. Pipek and P. G. Mezey, J. Chem. Phys. 90, 4916 (1989).
  * Edmiston and Ruedenberg [1963] C. Edmiston and K. Ruedenberg, Rev. Mod. Phys. 35, 457 (1963).
  * von Niessen [1972] W. von Niessen, J. Chem. Phys. 56, 4290 (1972).
  * Foster and Boys [1960] J. M. Foster and S. F. Boys, Rev. Mod. Phys. 32, 300 (1960).
  * Høyvik  _et al._ [2012] I.-M. Høyvik, B. Jansik, and P. Jørgensen, J. Chem. Phys. 137, 224114 (2012).
  * Pederson  _et al._ [2014] M. R. Pederson, A. Ruzsinszky, and J. P. Perdew, J. Chem. Phys. 140, 121103 (2014).
  * Pederson [2015] M. R. Pederson, J. Chem. Phys. 142, 064112 (2015).
  * Hahn  _et al._ [2017] T. Hahn, S. Schwalbe, J. Kortus, and M. R. Pederson, J. Chem. Theory Comput. 13, 5823 (2017).
  * Yang  _et al._ [2017] Z. Yang, M. R. Pederson, and J. P. Perdew, Phys. Rev. A 95, 052505 (2017).
  * Pederson and Baruah [2015] M. R. Pederson and T. Baruah, Adv. At. Mol. Opt. Phys. , 153 (2015).
  * Sharkas  _et al._ [2018] K. Sharkas, L. Li, K. Trepte, K. P. K. Withanage, R. P. Joshi, R. R. Zope, T. Baruah, J. K. Johnson, K. A. Jackson, and J. E. Peralta, J. Phys. Chem. A 122, 9307 (2018), pMID: 30412407.
  * Withanage  _et al._ [2021] K. P. K. Withanage, P. Bhattarai, J. E. Peralta, R. R. Zope, T. Baruah, J. P. Perdew, and K. A. Jackson, J. Chem. Phys. 154, 024102 (2021).
  * Pederson  _et al._ [2016] M. R. Pederson, T. Baruah, D. you Kao, and L. Basurto, J. Chem. Phys. 144, 164117 (2016).
  * Goedecker and Umrigar [1997] S. Goedecker and C. J. Umrigar, Phys. Rev. A 55, 1765 (1997).
  * Stengel and Spaldin [2008] M. Stengel and N. A. Spaldin, Phys. Rev. B 77, 155106 (2008).
  * Szotek  _et al._ [1993] Z. Szotek, W. M. Temmerman, and H. Winter, Phys. Rev. B 47, 4029 (1993).
  * Arai and Fujiwara [1995] M. Arai and T. Fujiwara, Phys. Rev. B 51, 1477 (1995).
  * Shinde  _et al._ [2020] R. Shinde, S. S. Yamijala, and B. M. Wong, J. Phys.: Condens. Matter 33, 115501 (2020).
  * Perdew  _et al._ [1996a] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996a).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Pauling [1941] L. Pauling, J. Am. Pharm. Assoc. 30, 30 (1941).
  * Miyake and Aryasetiawan [2008] T. Miyake and F. Aryasetiawan, Phys. Rev. B 77, 085122 (2008).
  * Vydrov and Scuseria [2004] O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. 121, 8187 (2004).
  * Klüpfel  _et al._ [2012] S. Klüpfel, P. Klüpfel, and H. Jónsson, J. Chem. Phys. 137, 124102 (2012).
  * Perdew  _et al._ [2014] J. P. Perdew, A. Ruzsinszky, J. Sun, and K. Burke, J. Chem. Phys. 140, 18A533 (2014).
  * Yamamoto  _et al._ [2019] Y. Yamamoto, C. M. Diaz, L. Basurto, K. A. Jackson, T. Baruah, and R. R. Zope, J. Chem. Phys. 151, 154105 (2019).
  * Sharkas  _et al._ [2020] K. Sharkas, K. Wagle, B. Santra, S. Akter, R. R. Zope, T. Baruah, K. A. Jackson, J. P. Perdew, and J. E. Peralta, Proc. Natl. Acad. Sci. 117, 11283 (2020).
  * Blum  _et al._ [2009] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Comput. Phys. Commun. 180, 2175 (2009).
  * Ren  _et al._ [2012a] X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter, and M. Scheffler, New J. Phys. 14, 053020 (2012a).
  * van Setten  _et al._ [2015] M. J. van Setten, F. Caruso, S. Sharifzadeh, X. Ren, M. Scheffler, F. Liu, J. Lischner, L. Lin, J. R. Deslippe, S. G. Louie, C. Yang, F. Weigend, J. B. Neaton, F. Evers, and P. Rinke, J. Chem. Theory Comput. 11, 5665 (2015).
  * Seidl  _et al._ [1996] A. Seidl, A. Görling, P. Vogl, J. A. Majewski, and M. Levy, Phys. Rev. B 53, 3764 (1996).
  * Zahariev  _et al._ [2013] F. Zahariev, S. S. Leang, and M. S. Gordon, J. Chem. Phys. 138, 244108 (2013).
  * Görling and Levy [1997] A. Görling and M. Levy, J. Chem. Phys. 106, 2675 (1997).
  * Jmol development team, Jmol, Available at http://jmol.sourceforge.net/ [2016] Jmol development team, Jmol, Available at http://jmol.sourceforge.net/,  version 14.6.4 (2016).
  * Note [2] The SIC-SCAN calculations for He2 were carried out with FHI-aims using the tight tier-1 NAO basis set.
  * Schuch  _et al._ [1958] A. F. Schuch, E. R. Grilly, and R. L. Mills, Phys. Rev. 110, 775 (1958).
  * Pemmaraju  _et al._ [2007] C. D. Pemmaraju, T. Archer, D. Sánchez-Portal, and S. Sanvito, Phys. Rev. B 75, 045101 (2007).
  * Shang  _et al._ [2017] H. Shang, C. Carbogno, P. Rinke, and M. Scheffler, Comput. Phys. Commun. 215, 26 (2017).
  * Klüpfel  _et al._ [2011] S. Klüpfel, P. Klüpfel, and H. Jónsson, Phys. Rev. A 84, 050501 (2011).
  * Shahi  _et al._ [2019] C. Shahi, P. Bhattarai, K. Wagle, B. Santra, S. Schwalbe, T. Hahn, J. Kortus, K. A. Jackson, J. E. Peralta, K. Trepte, S. Lehtola, N. K. Nepal, H. Myneni, B. Neupane, S. Adhikari, A. Ruzsinszky, Y. Yamamoto, T. Baruah, R. R. Zope, and J. P. Perdew, J. Chem. Phys. 150, 174102 (2019).
  * Zheng  _et al._ [2011] X. Zheng, A. J. Cohen, P. Mori-Sánchez, X. Hu, and W. Yang, Phys. Rev. Lett. 107, 026403 (2011).
  * Atalla  _et al._ [2016] V. Atalla, I. Y. Zhang, O. T. Hofmann, X. Ren, P. Rinke, and M. Scheffler, Phys. Rev. B 94, 035140 (2016).
  * Yang  _et al._ [2016] Z.-h. Yang, H. Peng, J. Sun, and J. P. Perdew, Phys. Rev. B 93, 205205 (2016).
  * R. Johnson, Computational Chemistry Comparison and Benchmark Database [2002] R. Johnson, Computational Chemistry Comparison and Benchmark Database,  NIST Standard Reference Database 101 (2002).
  * Godby  _et al._ [1986] R. W. Godby, M. Schlüter, and L. J. Sham, Phys. Rev. Lett. 56, 2415 (1986).
  * Kim  _et al._ [2015] M.-C. Kim, H. Park, S. Son, E. Sim, and K. Burke, J. Phys. Chem. Lett. 6, 3802 (2015), pMID: 26722874\.
  * Perdew  _et al._ [1996b] J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105, 9982 (1996b).
  * Li  _et al._ [2017] C. Li, X. Zheng, N. Q. Su, and W. Yang, Natl. Sci. Rev. 5, 203 (2017).
  * Zhang  _et al._ [2019] I. Y. Zhang, A. J. Logsdail, X. Ren, S. V. Levchenko, L. Ghiringhelli, and M. Scheffler, New J. Phys. 21, 013025 (2019).
  * Wang  _et al._ [2021] Y. Wang, Y. Li, J. Chen, I. Y. Zhang, and X. Xu, JACS Au 1, 543 (2021).
  * Xiao  _et al._ [2013] B. Xiao, J. Sun, A. Ruzsinszky, J. Feng, R. Haunschild, G. E. Scuseria, and J. P. Perdew, Phys. Rev. B 88, 184103 (2013).
  * Mahan [1990] G. D. Mahan,  _Many-Particle Physics_ (Springer US, 1990).
  * Zacharias  _et al._ [2020] M. Zacharias, M. Scheffler, and C. Carbogno, Phys. Rev. B 102 (2020).
  * Brown-Altvater  _et al._ [2020] F. Brown-Altvater, G. Antonius, T. Rangel, M. Giantomassi, C. Draxl, X. Gonze, S. G. Louie, and J. B. Neaton, Phys. Rev. B 101, 165102 (2020).
  * de Boer and Verwey [1937] J. H. de Boer and E. J. W. Verwey, Proc. Phys. Soc. 49, 59 (1937).
  * Hubbard and Flowers [1963] J. Hubbard and B. H. Flowers, Proc. R. Soc. A: Math. Phys. Eng. Sci. 276, 238 (1963).
  * Brandow [1977] B. Brandow, Adv. Phys. 26, 651 (1977).
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * Jiang  _et al._ [2010] H. Jiang, R. I. Gomez-Abal, P. Rinke, and M. Scheffler, Phys. Rev. B 82, 045108 (2010).
  * Floris  _et al._ [2011] A. Floris, S. de Gironcoli, E. K. U. Gross, and M. Cococcioni, Phys. Rev. B 84, 161102 (2011).
  * Jiang [2018] H. Jiang, Phys. Rev. B 97, 245132 (2018).
  * Shishkin and Kresse [2007] M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007).
  * Surh  _et al._ [1995] M. P. Surh, H. Chacham, and S. G. Louie, Phys. Rev. B 51, 7464 (1995).
  * Vydrov  _et al._ [2006] O. A. Vydrov, G. E. Scuseria, J. P. Perdew, A. Ruzsinszky, and G. I. Csonka, J. Chem. Phys. 124, 094108 (2006).
  * Aquino  _et al._ [2020] F. W. Aquino, R. Shinde, and B. M. Wong, J. Comput. Chem. 41, 1200 (2020).
  * Bylaska  _et al._ [2006] E. J. Bylaska, K. Tsemekhman, and F. Gao, Phys. Scr. 2006, 86 (2006).
  * Kahlert  _et al._ [1987] H. Kahlert, O. Leitner, and G. Leising, Synth. Met. 17, 467 (1987).
  * Fincher  _et al._ [1982] C. R. Fincher, C. E. Chen, A. J. Heeger, A. G. MacDiarmid, and J. B. Hastings, Phys. Rev. Lett. 48, 100 (1982).
  * Atala  _et al._ [2013] M. Atala, M. Aidelsburger, J. T. Barreiro, D. Abanin, T. Kitagawa, E. Demler, and I. Bloch, Nat. Phys. 9, 795 (2013).
  * Ooshika [1957a] Y. Ooshika, J. Phys. Soc. Jpn. 12, 1238 (1957a).
  * Ooshika [1957b] Y. Ooshika, J. Phys. Soc. Jpn. 12, 1246 (1957b).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979).
  * Hirata  _et al._ [1998] S. Hirata, H. Torii, and M. Tasumi, Phys. Rev. B 57, 11994 (1998).
  * Becke [1993] A. D. Becke, J. Chem. Phys. 98, 1372 (1993).
  * Ren  _et al._ [2012b] X. Ren, P. Rinke, C. Joas, and M. Scheffler, J. Mater. Sci. 47, 7447 (2012b).
  * Körzdörfer  _et al._ [2012] T. Körzdörfer, R. M. Parrish, J. S. Sears, C. D. Sherrill, and J.-L. Brédas, J. Chem. Phys. 137, 124305 (2012).
  * Knuth [2015] F. Knuth,  _Strain and Stress: Derivation, Implementation, and Application to Organic Crystals_ , Ph.D. thesis, Freie Universität Berlin (2015).