# iOI: an Iterative Orbital Interaction Approach for Solving the Self-
Consistent Field Problem

Zikuan Wang  Qingdao Institute for Theoretical and Computational Sciences,
Institute of Frontier and Interdisciplinary Science, Shandong University,
Qingdao, Shandong 266237, P. R. China  Wenjian Liu<EMAIL_ADDRESS>Qingdao
Institute for Theoretical and Computational Sciences, Institute of Frontier
and Interdisciplinary Science, Shandong University, Qingdao, Shandong 266237,
P. R. China

###### Abstract

An iterative orbital interaction (iOI) approach is proposed to solve, in a
bottom-up fashion, the self-consistent field problem in quantum chemistry.
While it belongs grossly to the family of fragment-based quantum chemical
methods, iOI is distinctive in that (1) it divides and conquers not only the
energy but also the wave function, and that (2) the subsystems sizes are
automatically determined by successively merging neighboring small subsystems
until they are just enough for converging the wave function to a given
accuracy. Orthonormal occupied and virtual localized molecular orbitals are
obtained in a natural manner, which can be used for all post-SCF purposes.

keywords: iterative orbital interaction, self-consistent field, bottom-up
solver, fragment localized molecular orbitals, modified projected hybrid
orbitals, automatic fragmentation and merging, iterative matrix block-
diagonalization

##  1 Introduction

The molecular orbitals (MO) of a chemical system resulting from a self-
consistent field (SCF) calculation can always be localized as long as there
exists a gap between the occupied and virtual spaces.1, 2 Because of this, the
electronic structure of a local part of a gapped large molecule (macromolecule
or molecular cluster) is not much affected by those parts that are far away in
distance. Such locality (also called “nearsightedness”3) is fundamental to all
low-scaling quantum chemical methods. The issue is how to materialize such
locality in a bottom-up fashion (i.e., without invoking the expensive top-down
localization1, 2 of the global occupied and virtual canonical MOs). As a
matter of fact, the best answer to this was already provided by chemists in
the early days of chemistry: as the basic building units, functional groups or
fragments in general reflect best the locality of chemical systems. As such,
various fragment-based quantum chemical methods have been developed in the
last decades, including mixed quantum mechanics/molecular mechanics (QM/MM)4,
5, 6, 7, multilayer QM/QM8, 9, 10, 11, divide-and-conquer (DC)12, 13, 14, 15,
embedding16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, and
orbital-32, 33, 34, 35, 36, 37, 38, 39, 40, 41 and energy-based42, 43, 44, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61 fragmentation
schemes, which differ mainly in how to divide the entire system into fragments
and how to conquer their mutual interactions, so as to achieve high accuracy
yet with low cost. For instance, in the energy-based fragmentation approaches,
the electronic energy and energy-related properties of the entire system are
obtained simply by assembling the fragmental values in one way or another62,
63, 64, 65, 66, 67, 68, 69, 70. The accuracy can further be improved greatly
by invoking a low-level treatment of the entire system71 or by using a
superposition of multiple fragmentation topologies72. Of course, instead of
making direct use of such chemical intuitions, one can also exploit locality
by purely mathematical techniques such as prescreening73, 74, 75, 76.

It is not surprising that the above fragment-based methods can also be used to
describe low-lying excited states that are localized in a small region of the
system77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89. Low-lying
delocalized excited states (excitons or transitions between delocalized
orbitals) can be accessed either by linear combinations of local states90, 91,
92, 93, 94, 95, 96 or by enlarged buffer regions97. Short-range charge
transfer (CT) effects can also be accounted for deliberately to improve those
states dominated by local excitations98. Given such successes, it must be
admitted that such descriptions of excited states all rely on some information
known in advance. It is often the case that the so-calculated excitation
energies are sufficiently good but transition properties that are sensitive to
the tails of the wave functions are still problematic. Moreover, genuine CT
excitations are missed by construction by fragment-based methods, needless to
say high-lying states, the radial extensions of which are usually far beyond
the intuitively partitioned fragments. For such types of states, one should
generate converged localized molecular orbitals (LMO) in a bottom-up fashion
(to stay within the fragmentation picture)1 and make proper use of both the
energetic and spatial localities of particle-hole pairs (the real basis of
post-SCF approaches) for describing excited states of a given energy range99.

Stimulated by the iterative configuration interaction (iCI) approach100, 101,
102, 103 for solving the many-electron Schrödinger equation, we propose here
an iterative orbital interaction (iOI) approach for solving the Hartree-
Fock/Kohn-Sham (HF/KS) equation in a bottom-up fashion. Since the HF/KS
equation represented in a finite basis can, at each iteration, be viewed as a
one-electron full CI problem, the same strategy employed by iCI, i.e.,
projecting the Hamiltonian onto a small set of physically oriented basis
functions (linear combinations of determinants or configuration state
functions, the number of which grows along macroiterations), can also be used
in iOI, in which the “physically oriented basis functions” are just linear
combinations of primitive fragment LMOs (pFLMO) from subsystem
calculations104, 105. By merging successively neighboring small fragments into
larger ones along the macroiterations, the pFLMOs are improved steadily and
many of them are essentially converged. They provide good initial guess for
the final global SCF calculation such that the number of SCF iterations is
reduced greatly. In short, iOI belongs grossly to the family of fragment-based
methods but is distinctive in that it divides and conquers both the energy and
wave function (orbitals). Orthonormal occupied and virtual LMOs are obtained
in a natural manner104, 105, which can be used for all post-SCF purposes. The
algorithm will be detailed in Sec. 2. Some showcases are then provided in Sec.
3 to reveal the efficacy of iOI.

##  2 iOI

The aim here is to determine the converged, orthonormal LMOs $\\{\psi_{p}\\}$
of a gapped molecule in a bottom-up fashion. As usual, they are expanded in
terms of atom-centered Gaussian basis functions $\\{\chi_{\mu}\\}$,

$\psi_{p}(\bm{r})=\sum_{\mu}^{M_{p}}\chi_{\mu}(\bm{r})C_{\mu p}[V_{eff}].$ (1)

The questions are how to identify a priori those atomic orbitals (AO)
$\\{\chi_{\mu}\\}$ that are most relevant for an LMO $\psi_{p}(\bm{r})$
centered at $\bm{G}_{p}=\langle\psi_{p}|\bm{r}|\psi_{p}\rangle$ with a radius
of $R_{p}$, and how to determine the corresponding expansion coefficients
$\mathbf{C}_{p}$ that depend on the effective potential $V_{eff}(\bm{r})$
under the constraints of the Brillouin conditions. Since every basis function
$\chi_{\mu}$ centered at $\bm{G}_{\mu}$ also has a finite radius $R_{\mu}$ of
distribution, it is clear that only those functions with distances
$|\bm{G}_{\mu}-\bm{G}_{p}|$ smaller than $R_{\mu}+R_{p}$ contribute to the
expansion (1). However, such information cannot be employed to truncate a
priori the expansion (1) since neither the center $\bm{G}_{p}$ nor the radius
$R_{p}$ of $\psi_{p}$ is known in advance. Nonetheless, a good guess on the
center $\bm{G}_{p}$ of $\psi_{p}$ stems naturally from the concept of
functional group or fragment in general. Therefore, the issue becomes how to
merge the initially partitioned fragments in a way that is just enough to
determine the right basis functions for expanding the set of LMOs localized on
a given fragment. As for the determination of the expansion coefficients
$\mathbf{C}_{p}$, the short-range exchange-correlation part of the effective
potential $V_{eff}(\bm{r})$ is taken into account automatically during the
merging process, whereas the long-range electrostatic interaction from the
environment can be accounted for explicitly by means of, e.g., the multipolar
expansion technique106, 107. Nevertheless, when all LMOs of the system are to
be targeted, it is more convenient and accurate to account for the long-range
interactions by a few global SCF iterations.

The above analysis suggests the following algorithm for the generation of the
LMOs of the whole system:

  1. (1)

Partition the whole system into nonoverlapping functional groups (e.g. CH3,
CH2, C=O, SH, C-Cl, etc.), each of which is obtained by cutting the single
bonds between main group elements other than hydrogen, halogens, and rare
gases;

  2. (2)

Initialize the macroiteration count $\tilde{m}$ to 0. Join the functional
groups defined in Step (1) to form $N_{\mathcal{F}}^{(\tilde{m})}$
nonoverlapping fragments $\mathcal{F}_{L}^{(\tilde{m})}$
($L\in[1,N_{\mathcal{F}}^{(\tilde{m})}]$) according to chemical intuition, so
that every atom in the whole system belongs to one and only one fragment (see
Fig.1 for an illustration). The nonoverlapping fragments are called here
primitive fragments (pFrag), which form precisely the whole system and are
typically of 10 to 30 atoms for the best cost-performance ratio;

Figure 1: Fragmentation of a short peptide. (a) Illustration of the iOI
merging process, where the buffer regions and the PHO (projected hybrid
orbital) boundary atoms are designated by blue and red circles, respectively.
(b) Definitions of $\mathcal{F}_{L}^{(\tilde{m})}$,
$\mathcal{B}_{L}^{(\tilde{m})}$, $\mathcal{C}_{L}^{(\tilde{m})}$, and
$\mathcal{P}_{L}^{(\tilde{m})}(\cdot)$.

  3. (3)

Cap each pFrag with parts of the system that are sufficiently close to the
pFrag (see Appendix A for the detailed algorithm). The atoms added in this
step to a pFrag $\mathcal{F}_{L}^{(0)}$ are collectively called buffer
$\mathcal{B}_{L}^{(0)}$;

  4. (4)

Saturate the dangling bonds to form closed subsystems $\mathcal{S}_{L}^{(0)}$.
Although a natural choice for this purpose is to use hydrogen atoms as link
atoms as in our previous work104, 105, it appears more convenient to use the
PHO (projected hybrid orbital)108 type of boundary atoms. Specifically, the
atoms directly bonded to buffer $\mathcal{B}_{L}^{(0)}$ are taken as boundary
atoms. Since $\mathcal{B}_{L}^{(0)}$ has been chosen in such a way that the
boundary atoms are always $p$-block element atoms, there always exist one $s$
orbital and three $p$ orbitals for each boundary atom, which can be hybridized
to give four hybrid orbitals (see Appendix B for the modification of the
PHO108). Only one of the hybrid orbitals will be available for bonding with
one atom in $\mathcal{B}_{L}^{(0)}$, while the remaining three hybrid orbitals
(designated as auxiliary PHO) are assigned with suitable occupation numbers
and frozen in the SCF calculation. The union of $\mathcal{B}_{L}^{(0)}$ and
its boundary atoms is called cap $\mathcal{C}_{L}^{(0)}$. A subsystem is hence
the union of a pFrag and its cap, i.e.,
$\mathcal{S}_{L}^{(0)}=\mathcal{F}_{L}^{(0)}\cup\mathcal{C}_{L}^{(0)}$.

  5. (5)

Carry out standard SCF calculation for each subsystem (designated as
microiteration), where the auxiliary PHOs are frozen (i.e., projected out of
the variational space; see Appendix B);

  6. (6)

Localize the CMOs of each subsystem using, e.g., the Boys localization109.
Note that the auxiliary PHOs are not included in the localization, since they
are not part of the subsystem but only serve to provide a realistic effective
potential during the SCF procedure;

  7. (7)

Enlarge successively the subsystems until they are just large enough for
expanding the pFLMOs up to a given precision:

    1. (a)

For each subsystem, pick up those occupied and virtual LMOs with the Löwdin
populations on its pFrag larger than $\eta_{\mathcal{F}}=0.1$, leading to
pFLMOs $\\{\psi^{\mathrm{o,FLMO}}_{\mu}\\}_{L}^{(\tilde{m})}$ and
$\\{\psi^{\mathrm{v,FLMO}}_{\mu}\\}_{L}^{(\tilde{m})}$ for the occupied and
virtual subspaces, respectively;

    2. (b)

Merge neighboring pFrags of macroiteration $\tilde{m}$ and cap the results to
yield larger subsystems. First, construct $N_{\mathcal{F}}^{(\tilde{m}+1)}$
pFrags for macroiteration $\tilde{m}+1$ such that the $L$th pFrag
$\mathcal{F}_{L}^{(\tilde{m}+1)}$ is just the union of exactly two proximal
pFrags (denoted as
$\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(1)}^{(\tilde{m})}$ and
$\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(2)}^{(\tilde{m})}$) of
macroiteration $\tilde{m}$, i.e.,
$\mathcal{F}_{L}^{(\tilde{m}+1)}=\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(1)}^{(\tilde{m})}\bigcup\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(2)}^{(\tilde{m})}$.
Note that every pFrag of macroiteration $\tilde{m}$ is used exactly one time.
The pFrags $\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(I)}^{(\tilde{m})}$
($I=1,2$) are called parents of $\mathcal{F}_{L}^{(\tilde{m}+1)}$. For
example, pFrags $\mathcal{F}_{4}^{(0)}$ and $\mathcal{F}_{5}^{(0)}$ in Fig. 1
are merged to give pFrag $\mathcal{F}_{4}^{(1)}$ (i.e.,
$\mathcal{P}_{4}^{(1)}(1)=4$ and $\mathcal{P}_{4}^{(1)}(2)=5$). Note in
passing that it is possible to merge more than two pFrags at a time. However,
the subsystem sizes may increase too rapidly as the macroiteration proceeds.
Second, cap the pFrags $\mathcal{F}_{L}^{(\tilde{m}+1)}$ in the same way as
described in Steps (3) and (4) but with increased buffer radii (see Appendix
C), thereby leading to $N_{\mathcal{F}}^{(\tilde{m}+1)}$ subsystems
$\mathcal{S}_{L}^{(\tilde{m}+1)}$ for macroiteration $\tilde{m}+1$. Among the
cap atoms added in this step, those that do not belong to
$\mathcal{S}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(\cdot)}^{(\tilde{m})}$
constitute the incremental cap of subsystem $\mathcal{S}_{L}^{(\tilde{m}+1)}$,
$\Delta\mathcal{C}_{L}^{(\tilde{m}+1)}=\mathcal{C}_{L}^{(\tilde{m}+1)}\setminus\bigcup_{I}\mathcal{S}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(I)}^{(\tilde{m})}$,
where $\mathcal{C}_{L}^{(\tilde{m}+1)}$ is the cap of
$\mathcal{S}_{L}^{(\tilde{m}+1)}$. Increase the macroiteration count
$\tilde{m}$ by 1.

The above merging scheme improves the quality of the pFLMOs by using more AOs
and improved effective potential $V_{eff}(\bm{r})$. Yet, there exist several
exceptions to these rules:

      1. i.

It is obvious that, if the AOs of subsystem $\mathcal{S}_{L}^{(\tilde{m})}$
are already sufficient for accurately expanding the pFLMOs,
$\mathcal{S}_{L}^{(\tilde{m})}$ need not be further enlarged. This can readily
be checked by the total Löwdin population, $(\Delta_{e})_{L}^{(\tilde{m})}$,
of the occupied pFLMOs $\\{\psi^{\mathrm{o,FLMO}}_{\mu}\\}_{L}^{(\tilde{m})}$
on the incremental cap $\Delta\mathcal{C}_{L}^{(\tilde{m})}$,

$\displaystyle(\Delta_{e})_{L}^{(\tilde{m})}$ $\displaystyle=$
$\displaystyle\sum_{\nu\in\Delta\mathcal{C}_{L}^{(\tilde{m})}}\left((\mathcal{S}_{L}^{(\tilde{m})})^{1/2}\mathbf{\tilde{D}}_{L}^{(\tilde{m})}(\mathcal{S}_{L}^{(\tilde{m})})^{1/2}\right)_{\nu\nu}$
(2) $\displaystyle=$
$\displaystyle\sum_{\nu\in\Delta\mathcal{C}_{L}^{(\tilde{m})}}\sum_{i\in\\{\psi^{\mathrm{o,FLMO}}_{\mu}\\}_{L}^{(\tilde{m})}}\left|\left((\mathcal{S}_{L}^{(\tilde{m})})^{1/2}\mathbf{C}_{L}^{(\tilde{m})}\right)_{\nu
i}\right|^{2},$
$\displaystyle(\mathbf{\tilde{D}}_{L}^{(\tilde{m})})_{\nu\kappa}$
$\displaystyle=$
$\displaystyle\sum_{i\in\\{\psi^{\mathrm{o,FLMO}}_{\mu}\\}_{L}^{(\tilde{m})}}(\mathbf{C}_{L}^{(\tilde{m})})_{\nu
i}(\mathbf{C}_{L}^{(\tilde{m})})_{\kappa i},$ (3)

where $\mathcal{S}_{L}^{(\tilde{m})}$ and $\mathbf{C}_{L}^{(\tilde{m})}$ are
the AO overlap matrix and the LMO coefficients of
$\mathcal{F}_{L}^{(\tilde{m})}$, respectively. If
$(\Delta_{e})_{L}^{(\tilde{m})}<\eta_{\mathrm{tail}}=0.1$, subsystem
$\mathcal{S}_{L}^{(\tilde{m})}$ is considered to be converged and remains
unchanged in subsequent macroiterations. Such converged subsystems have just
enough sizes. For example, subsystems $\mathcal{S}_{1}^{(1)}$,
$\mathcal{S}_{3}^{(1)}$ and $\mathcal{S}_{4}^{(1)}$ in Fig. 1 are converged at
macroiteration 1, but $\mathcal{S}_{2}^{(1)}$ is not yet so and has to be
enlarged once more to give $\mathcal{S}_{2}^{(2)}$.

      2. ii.

When the total number of nonconvergent subsystems is odd or when a pFrag at
macroiteration $\tilde{m}$ is spatially distant from the remaining ones, it is
impossible to merge all pFrags in a pairwise fashion. In these cases, some
pFrags of macroiteration $\tilde{m}+1$ may contain only one pFrag of
macroiteration $\tilde{m}$, i.e.,
$\mathcal{F}_{L}^{(\tilde{m}+1)}=\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(1)}^{(\tilde{m})}$.
However, $\mathcal{S}_{L}^{(\tilde{m}+1)}$ is still larger than
$\mathcal{S}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(1)}^{(\tilde{m})}$ since
$\mathcal{C}_{L}^{(\tilde{m}+1)}$ is larger than
$\mathcal{C}_{\mathcal{P}_{L}^{(\tilde{m}+1)}(1)}^{(\tilde{m})}$. For example,
at macroiteration 0, there are five pFrags for the example shown in Fig. 1.
For load balance, $\mathcal{F}_{2}^{(0)}$ is chosen as the pFrag that will not
undergo merging for $\mathcal{S}_{2}^{(0)}$ is the largest of the five
subsystems. This separates $\mathcal{F}_{1}^{(0)}$ from the other three
pFrags, such that $\mathcal{F}_{1}^{(0)}$ and $\mathcal{F}_{3}^{(0)}$ are also
directly carried over to the next macroiteration. Then, only
$\mathcal{F}_{4}^{(0)}$ and $\mathcal{F}_{5}^{(0)}$ are merged (for details,
see Appendix C).

      3. iii.

When all subsystems have been converged (large enough), merge them together
and stop the macroiterations. If all (nonconvegent) subsystems will be merged
together at a macroiteration, also stop the macroiterations. The iOI procedure
is then finished by performing a small number of global SCF iterations (see
Step (8)).

    3. (c)

Use the pFLMOs of macroiteration $\tilde{m}-1$ as initial guess for
macroiteration $\tilde{m}$. Specifically, find all subsystems
$\mathcal{S}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}$
($i\in[1,\mathcal{M}_{L}^{(\tilde{m})}])$ whose pFrags overlap with subsystem
$\mathcal{S}_{L}^{(\tilde{m})}$, i.e., find $\\{L_{i}^{(\tilde{m}-1)}\\}$ with
$\mathcal{S}_{L}^{(\tilde{m})}\cap\mathcal{F}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}\neq\varnothing$.
For example, in Fig. 1, we have $\mathcal{M}_{1}^{(1)}=2$, $1_{1}^{(0)}=1$ and
$1_{2}^{(0)}=2$, since $\mathcal{S}_{1}^{(1)}$ overlaps with
$\mathcal{F}_{1}^{(0)}$ and $\mathcal{F}_{2}^{(0)}$ but not with
$\mathcal{F}_{3}^{(0)}$, $\mathcal{F}_{4}^{(0)}$ or $\mathcal{F}_{5}^{(0)}$.
Note that
$\mathcal{P}_{L}^{(\tilde{m})}(\cdot)\in\\{L_{i}^{(\tilde{m}-1)}\\}$, since
$\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m})}(\cdot)}^{(\tilde{m}-1)}\subseteq\mathcal{F}_{L}^{(\tilde{m})}$,
which implies necessarily that
$\mathcal{F}_{\mathcal{P}_{L}^{(\tilde{m})}(\cdot)}^{(\tilde{m}-1)}$ overlaps
with $\mathcal{S}_{L}^{(\tilde{m})}$. The pFLMOs of subsystems
$\\{\mathcal{S}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}\\}$ are then projected
onto subsystem $\mathcal{S}_{L}^{(\tilde{m})}$, i.e.,

$(\mathbf{C}_{L}^{(\tilde{m})})_{\kappa
p}=\sum_{\nu\in\mathcal{S}_{L}^{(\tilde{m})}}((\mathcal{S}_{L}^{(\tilde{m})})^{-1})_{\kappa\nu}(\mathcal{S}_{L^{(\tilde{m})},L_{i}^{(\tilde{m}-1)}}\mathbf{C}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)})_{\nu
p},$
$p\in\\{\psi^{\mathrm{o,FLMO}}_{\mu}\\}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}\cup\\{\psi^{\mathrm{v,FLMO}}_{\mu}\\}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)},$
(4)

where $\mathcal{S}_{L}^{(\tilde{m})}$ and
$\mathcal{S}_{L^{(\tilde{m})},L_{i}^{(\tilde{m}-1)}}$ are the AO overlap
matrix of subsystem $\mathcal{S}_{L}^{(\tilde{m})}$ and the rectangular
overlap matrix between the AOs of $\mathcal{S}_{L}^{(\tilde{m})}$ and
$\mathcal{S}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}$, respectively.
Physically, the projection Eq. (4) can be interpreted as finding the pFLMOs
$\\{(\mathbf{C}_{L}^{(\tilde{m})})_{p}\\}$ that resemble the input pFLMOs
$\\{(\mathbf{C}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)})_{p}\\}$ in the least-
squares sense, i.e. the $(\mathbf{C}_{L}^{(\tilde{m})})_{p}$ defined in Eq.
(4) minimizes the differential inner product
$\langle(\mathbf{C}_{L}^{(\tilde{m})})_{p}-(\mathbf{C}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)})_{p}|(\mathbf{C}_{L}^{(\tilde{m})})_{p}-(\mathbf{C}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)})_{p}\rangle$
subject to the constraint that $\\{(\mathbf{C}_{L}^{(\tilde{m})})_{p}\\}$
forms an orthonormal set110. As mentioned above, the parental subsystems of
$\mathcal{S}_{L}^{(\tilde{m})}$,
$\mathcal{S}_{\mathcal{P}_{L}^{(\tilde{m})}(\cdot)}^{(\tilde{m}-1)}$, must be
members of subsystems
$\\{\mathcal{S}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}\\}$ since the pFrags
of the former overlap with $\mathcal{S}_{L}^{(\tilde{m})}$ by construction. Of
all pFLMOs from $\\{\mathcal{S}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}\\}$,
those from
$\mathcal{S}_{\mathcal{P}_{L}^{(\tilde{m})}(\cdot)}^{(\tilde{m}-1)}$ are
denoted as central pFLMOs of $\mathcal{S}_{L}^{(\tilde{m})}$, i.e.,
$\\{\xi_{\mu}^{\mathrm{o}}\\}_{L}^{(\tilde{m})}$ and
$\\{\xi_{\mu}^{\mathrm{v}}\\}_{L}^{(\tilde{m})}$, whereas the remaining ones
are called buffer pFLMOs of $\mathcal{S}_{L}^{(\tilde{m})}$. Note that, in the
presence of PHO boundary atoms, the contributions of the auxiliary PHOs should
be projected out from the central pFLMOs;

    4. (d)

For each subsystem, perform “local elimination of linear dependence” (LELD) in
$\\{\xi_{\mu}^{\mathrm{o}}\\}_{L}^{(\tilde{m})}$ and
$\\{\xi_{\mu}^{\mathrm{v}}\\}_{L}^{(\tilde{m})}$. In the original LELD
algorithm,1 we eliminate the $\xi^{\mathrm{o}}_{\mu}$ that has the largest
overlap with the eigenvector $\\{v_{p}\\}$ of the overlap matrix
$\langle\xi^{\mathrm{o}}_{\mu}|\xi^{\mathrm{o}}_{\nu}\rangle$ with the
smallest eigenvalue. In case of degeneracy, the $\xi_{\mu}^{\mathrm{o}}$ with
the largest orbital spread
$O[\xi_{\mu}^{\mathrm{o}}]=2(\langle\xi_{\mu}^{\mathrm{o}}|\widehat{r^{2}}|\xi_{\mu}^{\mathrm{o}}\rangle-|\langle\xi_{\mu}^{\mathrm{o}}|\hat{\mathbf{r}}|\xi_{\mu}^{\mathrm{o}}\rangle|^{2})$
is first eliminated. The overlap matrix of the survival pFLMOs is then re-
diagonalized, and the process is repeated until the number of remaining
occupied and virtual orbitals match those expected from the numbers of
electrons and basis functions of the subsystem. A naïve implementation of this
algorithm requires diagonalizing the global pFLMO overlap matrix
$N_{\mathrm{red}}$ times (NB: $N_{\mathrm{red}}$ is the number of redundant
pFLMOs), which scales as
$\mathcal{O}(N_{\mathrm{pFLMO}}^{3}N_{\mathrm{red}})$. However, the pFLMO
overlap matrix has a blocked structure, so that its diagonalization can be
approximately replaced by the diagonalization of $O(N_{\mathrm{pFLMO}})$ small
matrices, each of dimension $O(1)$. The revised LELD algorithm thus scales
linearly with respect to system size (the detailed algorithm will be reported
in a separate paper). The linearly independent occupied pFLMOs are then Löwdin
orthonormalized. The virtual pFLMOs are treated in the same way after
projecting out the orthonormal occupied pFLMOs;

    5. (e)

Do subsystem SCF calculations. To retain the locality of the LMOs throughout
the SCF calculation, we here use the bottom-up, least-change block-
diagonalization1 instead of the conventional full diagonalization of the Fock
matrix $\mathbf{F}$ represented in the pFLMO basis and written in block form,

$\mathbf{F}=\begin{pmatrix}\mathbf{F}_{\mathrm{\tilde{o}\tilde{o}}}&\mathbf{F}_{\mathrm{\tilde{o}v}}\\\
\mathbf{F}_{\mathrm{v\tilde{o}}}&\mathbf{F}_{\mathrm{vv}}\end{pmatrix},$ (5)

where the subscript $\mathrm{\tilde{o}}$ emphasizes that the occupied buffer
pFLMOs and the auxiliary PHOs are all frozen. While freezing the auxiliary
PHOs is required by the PHO method itself108, the rationale for freezing also
the occupied buffer pFLMOs of subsystem $\mathcal{S}_{L}^{(\tilde{m})}$ is
that they come from the pFrags
$\mathcal{F}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}$ of neighboring
subsystems, where they experience a very realistic environment provided by the
caps $\mathcal{C}_{L_{i}^{(\tilde{m}-1)}}^{(\tilde{m}-1)}$ of these
subsystems. In contrast, as the cap of $\mathcal{S}_{L}^{(\tilde{m})}$, the
occupied buffer pFLMOs experience an inaccurate effective potential.
Therefore, their relaxation would actually worsen themselves and in turn also
worsen the central pFLMOs of $\mathcal{S}_{L}^{(\tilde{m})}$. Following the
same mathematics as we go from the matrix Dirac equation to the exact two-
component relativistic theory111, 112, we can define the following unitary
transformation113

$\displaystyle\mathbf{U}$ $\displaystyle=$
$\displaystyle\mathbf{D}\mathbf{P},$ (6) $\displaystyle\mathbf{D}$
$\displaystyle=$
$\displaystyle\begin{pmatrix}\mathbf{I}_{\mathrm{\tilde{o}\tilde{o}}}&-\mathbf{X}^{\dagger}\\\
\mathbf{X}&\mathbf{I}_{\mathrm{vv}}\end{pmatrix},$ (7)
$\displaystyle\mathbf{P}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\tilde{\mathcal{S}}_{\mathrm{\tilde{o}\tilde{o}}}^{-1/2}&\mathbf{0}\\\
\mathbf{0}&\tilde{\mathcal{S}}_{\mathrm{vv}}^{-1/2}\end{pmatrix},$ (8)
$\displaystyle\tilde{\mathcal{S}}_{\mathrm{\tilde{o}\tilde{o}}}$
$\displaystyle=$
$\displaystyle\mathbf{I}_{\mathrm{\tilde{o}\tilde{o}}}+\mathbf{X}^{\dagger}\mathbf{X},\quad\tilde{\mathcal{S}}_{\mathrm{vv}}=\mathbf{I}_{\mathrm{vv}}+\mathbf{X}\mathbf{X}^{\dagger},$
(9)

which block-diagonalizes the Fock matrix (5) as

$\displaystyle\bar{\mathbf{F}}$ $\displaystyle=$
$\displaystyle\mathbf{U}^{\dagger}\mathbf{F}\mathbf{U}=\begin{pmatrix}\bar{\mathbf{F}}_{\mathrm{\tilde{o}\tilde{o}}}&\mathbf{0}\\\
\mathbf{0}&\bar{\mathbf{F}}_{\mathrm{vv}}\end{pmatrix},$ (10)
$\displaystyle\bar{\mathbf{F}}_{\mathrm{\tilde{o}\tilde{o}}}$ $\displaystyle=$
$\displaystyle\tilde{\mathcal{S}}_{\mathrm{\tilde{o}\tilde{o}}}^{-1/2}\tilde{\mathbf{F}}_{\mathrm{\tilde{o}\tilde{o}}}\tilde{\mathcal{S}}_{\mathrm{\tilde{o}\tilde{o}}}^{-1/2},$
(11) $\displaystyle\bar{\mathbf{F}}_{\mathrm{vv}}$ $\displaystyle=$
$\displaystyle\tilde{\mathcal{S}}_{\mathrm{vv}}^{-1/2}\tilde{\mathbf{F}}_{\mathrm{vv}}\tilde{\mathcal{S}}_{\mathrm{vv}}^{-1/2},$
(12) $\displaystyle\tilde{\mathbf{F}}_{\mathrm{\tilde{o}\tilde{o}}}$
$\displaystyle=$
$\displaystyle\mathbf{F}_{\mathrm{\tilde{o}\tilde{o}}}+\mathbf{F}_{\mathrm{\tilde{o}v}}\mathbf{X}+\mathbf{X}^{\dagger}\mathbf{F}_{\mathrm{v\tilde{o}}}+\mathbf{X}^{\dagger}\mathbf{F}_{\mathrm{vv}}\mathbf{X},$
(13) $\displaystyle\tilde{\mathbf{F}}_{\mathrm{vv}}$ $\displaystyle=$
$\displaystyle\mathbf{F}_{\mathrm{vv}}-\mathbf{X}\mathbf{F}_{\mathrm{\tilde{o}v}}-\mathbf{F}_{\mathrm{v\tilde{o}}}\mathbf{X}^{\dagger}+\mathbf{X}\mathbf{F}_{\mathrm{\tilde{o}\tilde{o}}}\mathbf{X}^{\dagger},$
(14) $\displaystyle\tilde{\mathbf{F}}_{\mathrm{v\tilde{o}}}$ $\displaystyle=$
$\displaystyle\mathbf{F}_{\mathrm{v\tilde{o}}}-\mathbf{X}\mathbf{F}_{\mathrm{\tilde{o}\tilde{o}}}+\mathbf{F}_{\mathrm{vv}}\mathbf{X}-\mathbf{X}\mathbf{F}_{\mathrm{\tilde{o}v}}\mathbf{X}=\tilde{\mathbf{F}}_{\mathrm{\tilde{o}v}}^{\dagger}=\mathbf{0}.$
(15)

Once the pFLMO Fock matrix (5) is available, the decoupling condition (15) for
$\mathbf{X}$ can be solved iteratively with the Jacobi sweep method. The
$\mathbf{U}$ matrix (6) is just the coefficient matrix of the LMOs in the
pFLMO basis (which is further expanded in the AO basis).

    6. (f)

go to Step (7a).

  8. (8)

Do global SCF. Since the target here is to have fully converged occupied and
virtual LMOs of the whole system, this global SCF or last macroiteration is
necessary.

    1. (a)

Gather the central pFLMOs of all subsystems and generate linearly independent,
orthonormal pFLMOs in the same way as described in Step (7d);

    2. (b)

Like subsystem calculations, the bottom-up, least-change block-diagonalization
of the Fock matrix is used. At each iteration the dimension of the Fock matrix
$\mathbf{F}$ (5) can be reduced by freezing the already converged LMOs: (1)
occupied LMOs $i$ that satisfy $\max_{a}|F_{ia}|<\eta_{\mathrm{Fov}}=10^{-4}$
are frozen, where $a$ runs over all virtual LMOs; (2) virtual LMOs $a$ that
satisfy $\max_{i}|F_{ia}|<\eta_{\mathrm{Fov}}$ are frozen, where $i$ runs over
the active occupied LMOs only. Since the off-diagonal elements $\\{F_{ia}\\}$
tend to diminish along the SCF iterations, usually more and more LMOs can be
frozen. However, some frozen LMOs of an iteration may no longer be frozen in
the subsequent iteration, such that the number of frozen LMOs may occasionally
decrease.

Some further remarks are in order.

  1. (I)

If the global SCF calculation (the last macroiteration of iOI) with a
sufficiently small threshold $\eta_{\mathrm{Fov}}$ for orbital freezing is
iterated till convergence, both the total energy and the wave function (LMOs)
will agree completely with those by the conventional SCF calculation.
Therefore, iOI is an exact method. Even in this case, for sufficiently large
systems, iOI can still achieve considerable speedup over the conventional SCF
algorithm because the gain in efficiency due to the much reduced number of
global SCF iterations overcompensates for the cost of the subsystem
calculations, which need not be fully converged (e.g., to $10^{-3}$ a.u. for
the energy change and $10^{-2}$ a.u. for the maximum change of the density
matrix elements) and can be made embarassingly parallel. Both occupied and
virtual LMOs of the global system are obtained for free at the end of the
global SCF calculation. This is important for the top-down, iterative
localization of the CMOs of large systems is very costly or even impossible.
Moreover, iOI offers an effective means to converge to a user-specified
configuration by assigning proper occupations of the pFrags105 when the global
system has more than one stable SCF solutions (which is very common in, e.g.,
antiferromagnetically coupled systems).

  2. (II)

iOI shares the same spirit as iCI100, 101 for solving the full CI problem: (a)
the orbitals of the small and larger subsystems in iOI correspond,
respectively, to the reference functions and singly and doubly excited
configurations in iCI; (b) the subsystem and global SCF calculations
correspond to the microiterations of iCI for relaxing the contraction
coefficients of the CI vectors; (c) both iOI and iCI can be viewed as a
particular sequential, exact (partial) diagonalization of a large matrix, by
getting first the roots of one portion of the matrix and then those of an
enlarged portion, until the full matrix has been sampled; (d) from the
mathematical point of view, the only distinction between iOI and iCI lies in
that the Fockian in the former depends on the unknowns whereas the Hamiltonian
in the latter does not; (d) if only the active part of an exceedingly large
systems is of interest, some selection of the orbitals localized on
neighboring buffers can be introduced to iOI, just like the selected iCI102,
103.

##  3 Results and discussion

As showcases, we herein report iOI calculations of a series of double-stranded
DNA molecules of various lengths, hereafter termed as
DNA${}_{n}(n=1,2,4,8,16)$, where one strand is composed of $n$ consecutive
adenosines (A) and the other is composed of $n$ thymidines (T). All
calculations were performed with the BDF program package106, 114, 115, 116,
107 on a server equipped with 16 Intel(R) Xeon(R) CPU E5-2640 v3 @ 2.60GHz
cores. The B3LYP functional117, 118 and the def2-SV(P)119 basis set were used
throughout. Unless otherwise noted, all SCF calculations of the global system
(including the exact SCF and the last macroiteration of iOI) are converged to
the following criteria: $<10^{-6}\mathrm{~{}a.u.}$ for the energy change and
$<10^{-4}\mathrm{~{}a.u.}$ for the maximum change of density matrix elements.
The subsystem calculations of DNAn were embarrassingly parallelized over
$N_{\mathrm{procs}}=\max(2n,16)$ processes, each with $16/N_{\mathrm{procs}}$
OpenMP threads, while the global system calculations were performed using 16
OpenMP threads.

The first molecule to be studied is the DNA16 molecule (1052 atoms, 10024
basis functions). The statistics of the iOI macroiterations is documented in
Table 1. It is first seen that iOI takes 4 macroiterations to converge (which,
of course, depends on the particular design of subsystem sizes). The number of
subsystems decreases along the macroiteration as a result of subsystem
merging, which is accompanied by an increased average size of the subsystems.
Since all subsystems are converged at macroiteration 2, all of them are merged
together at macroiteration 3 for final SCF iterations of the whole system.
Note that the subsystems at macroiteration 2 are still much smaller in size
than the whole system, thereby explaining the shorter timing of macroiteration
2 than that of the final macroiteration, needless to say the first two
macroiterations.

Table 1: Statistics of the iOI macroiterations for DNA16.

Macroiteration | 0 | 1 | 2 | 3a
---|---|---|---|---
No. of subsystems | 96 | 48 | 38 | 1
No. of converged subsystems | 0 | 21 | 38 | 1
Minimal No. of atoms in a subsystem | 7 | 36 | 36 | 1052
Maximum No. of atoms in a subsystem | 56 | 114 | 195 | 1052
Average No. of atoms in a subsystem | 26 | 76 | 101 | 1052
Wall time (s) | 1441 | 2476 | 5845 | 75683

a For the whole system.

The details of the final macroiteration of iOI is shown in Fig. 2 to compare
with the exact SCF calculation starting from the superposition of atomic
densities. It can first be seen that the last macroiteration of iOI converges
within 11 iterations, whereas the exact SCF calculation takes 23 iterations to
converge. The first three iOI macroiterations thus generate a much better
guess than the atomic guess and consequently save 12 iterations, at a cost
(9762 s) that corresponds roughly to one global SCF iteration (Fig. 2, bottom
left). The iOI total energy is in error of $4\times 10^{-5}$ a.u. ($4\times
10^{-8}$ a.u. per atom) compared to the exact SCF energy (Fig. 2, top left),
which is a result of freezing the nearly converged LMOs. The number of active
LMOs decreases rapidly as the SCF iterations proceed, approaching eventually
to only one fifth of the total number of LMOs (Fig. 2, top right), which
results in a marked speedup of the Fock matrix block-diagonalization over the
exact diagonalization (Fig. 2, bottom right). It is interesting to note that
the proportion of frozen virtual LMOs is much higher than that of frozen
occupied LMOs (Fig. 2, top right), which can be rationalized by the fact that
there are typically many more Rydberg orbitals than core orbitals with
extended basis sets.

Figure 2: Comparison of the final iOI macroiteration and exact SCF
calculations of DNA16. Top left: error per atom of the iOI energy compared
with the exact SCF energy; top right: number of active occupied and virtual
orbitals; bottom left: wall time of each SCF iteration; bottom right: wall
time of the Fock matrix (block-)diagonalization step of each SCF iteration.

To estimate the scaling behaviors of the energy error and computational time
with respect to system size, we performed iOI and exact SCF calculations on
the series of DNA${}_{n}(n=1,2,4,8,16)$, with the results summarized in Fig.
3. The error of the iOI energy per atom is very small for all systems studied
but tends to increase linearly with respect to system size (Fig. 3, top left).
The latter can be understood as follows: if every LMO is converged to an error
threshold $\epsilon$, the (long-range) Coulomb potential felt by a given LMO
will be in error approximately of $O(\epsilon N_{\mathrm{LMO}})$, thereby
leading to a linear dependence of the energy error per atom on the system
size. Intriguingly, the number of the exact SCF iterations from the atomic
guess increases essentially linearly with respect to the system size, while
the number of iterations of the final iOI macroiteration is almost independent
of the system size (Fig. 3, top right). This implies simply that the atomic
guess becomes worse and worse as the system enlarges but the fragmental guess
inherent in the iOI procedure has a constant quality. As a result, iOI has a
lower scaling than the exact SCF: while iOI is slower than the exact SCF for
DNA1 (62 atoms, 544 basis functions) and DNA2 (128 atoms, 1176 basis
functions), a crossover is observed at DNA4 (260 atoms, 2440 basis functions),
and iOI becomes about twice as fast as the exact SCF for DNA16 (1052 atoms,
10024 basis functions) (Fig. 3, bottom left). The results highlight the
frequently overlooked fact that, for an iterative method, the scaling of a
single iteration does not necessarily reflect the overall scaling of the
method itself, as the number of iterations can depend on the system size. In
this sense, iOI improves the overall scaling of SCF by a much reduced number
of iterations, even though it does not improve the scaling of individual SCF
iterations. Moreover, the block-diagonalization is noticeably faster than the
full diagonalization of the Fock matrix for the largest systems (Fig. 3,
bottom right), which contributes further (albeit marginally) to the speedup of
iOI over the exact SCF. It deserves to be mentioned that the cubic scaling of
the iOI Fock matrix block-diagonalization is the result of using dense matrix
algebra. As the matrices $\mathbf{F}$, $\mathbf{D}$ and $\mathbf{P}$ in Eqs.
(5), (7), and (8), respectively, are all numerically sparse, with only
$O(N_{\mathrm{bas}})$ non-negligible elements, the Fock matrix block-
diagonalization can be made linear scaling when sparse matrix algebra is used.
This is to be contrasted with the full diagonalization, which can at best be
made quadratic scaling even if the sparsity of $\mathbf{F}$ is explored: the
result of the full diagonalization, the CMO coefficient matrix, is a full
matrix with generally $O(N_{\mathrm{bas}}^{2})$ non-negligible elements, which
necessarily requires at least $O(N_{\mathrm{bas}}^{2})$ time to generate.

Figure 3: Comparison of iOI and exact SCF calculations of
DNA${}_{n}(n=1,2,4,8,16)$. Top left: error of the iOI energy compared with the
exact SCF energy; top right: number of SCF iterations of the final iOI
macroiteration and exact SCF; bottom left: total wall time; bottom right: wall
time of the Fock matrix (block-)diagonalization averaged over number of SCF
iterations.

Finally, we examine the quality of the iOI wave function (LMOs) by computing
the lowest 10 excited states of DNA4 with time-dependent density functional
theory. As shown in Fig. 4, the final iOI macroiteration, which is the most
expensive step of an iOI calculation for sufficiently large molecules, need
not be fully converged: 3 SCF iterations are already enough to generate
sufficiently good LMOs, in terms of which the absorption spectrum is
practically indistinguishable from that derived from the exact LMOs. Stopped
at this point, the iOI energy (-9989.825528 a.u.) is 112 mHartree (0.43
mHartree/atom) higher than the exact SCF energy (-9989.937405 a.u.); the total
elapsed time is 3566 s, about half that of the fully converged iOI (6636 s) or
the exact SCF calculation (6769 s). It can therefore be concluded that iOI not
only provides virtually exact SCF energies but also excellent SCF wave
functions (LMOs).

Figure 4: Electronic absorption spectra of DNA4 using LMOs from partially
converged iOI calculations and fully converged exact SCF calculation. All
transitions are Gaussian broadened with 0.2 eV for the full width at half
maximum. The numbers of SCF iterations in the final iOI macroiteration are
given in the legend.

##  4 Conclusions and outlook

A bottom-up, fragment-based SCF solver, iOI, has been developed. Unlike most
fragment-based methods which conquer only the energy approximately but not the
wave function at all, iOI has the following merits: (1) the starting small
fragments are merged automatically to become just large enough for the
description of the LMOs of the whole system; (2) the number of global SCF
iterations is essentially independent of system size, in contrast to the
traditional SCF approach which requires more iterations for large than for
small systems; (3) it can yield not only the exact SCF energy with notable
time savings but also fully converged occupied and virtual LMOs without
additional cost. Some additional merits of iOI, which will be demonstrated in
future, include (a) facilitation of SCF convergence in difficult situations,
(b) directed convergence to a user-specified configuration when there exist
multiple SCF solutions or when an excited configuration is desired, and (c)
automatic determination of the size of active space in multi-scale
calculations of exceedingly large systems.

## Acknowledgement

This work was supported by the National Key R&D Program of China (Grant No.
2017YFB0203402), National Natural Science Foundation of China (Grant Nos.
21833001 and 21973054), Mountain Tai Climb Program of Shandong Province, and
Key-Area Research and Development Program of Guangdong Province (Grant No.
2020B0101350001).

## Data Availability Statement

The data that supports the findings of this study are available within the
article.

##  Appendix A Automatic selection of buffer atoms for each fragment

In Steps (3) and (4) of the iOI algorithm, some atoms that are sufficiently
close to a pFrag are joined with the pFrag to form a subsystem, so as to
provide a more realistic environment for the pFrag on one hand and to saturate
the dangling bonds on the other. The term “sufficiently close” can however
adopt many possible definitions. The most straightforward way is to define
“close” as being spatially close. Thus, for each atom in the pFrag, one can
draw a sphere of a certain radius (hereafter called the buffer radius),
$R_{L}^{(\tilde{m})}$ (where $\tilde{m}$, the macroiteration number, is equal
to 0 here), centered at that atom. Any atom external to the fragment that
falls into at least one of the spheres

$A\in\mathcal{B}_{L}^{(\tilde{m})}\quad\mathrm{iff}\quad\exists
B\in\mathcal{F}_{L}^{(\tilde{m})}\quad\mathrm{s.t.}\quad
R_{AB}<R_{L}^{(\tilde{m})}$ (16)

is then chosen as a buffer atom. The buffering process is then followed by
saturating the remaining dangling bonds by PHO boundary atoms to give the cap
$\mathcal{C}_{L}^{(\tilde{m})}$, which is the union of the buffer atoms and
the PHO boundary atoms. This approach however has a few drawbacks:

  1. (1)

The importance of an atom as a buffer atom is not always correlated with its
spatial distance from the pFrag. Rather, orbital-based metrics are usually
more reliable than real space-based metrics. Since one of the purposes of
introducing the buffer atoms is to provide enough AOs for accurately expanding
the pFLMOs, the most rigorous metric would be the population of the pFLMOs on
the AOs of the buffer atoms, i.e., a buffer atom is only included if the
pFLMOs have a large enough population on that atom. However, the pFLMOs are
not known at the time the buffer atoms are chosen. Thus, we instead measure
the importance of a buffer atom by the overlap matrix elements of the AOs of
the buffer atom with the AOs of the pFrag. Specifically, we define the
following “effective distance” between two atoms:

$R_{AB}^{\mathrm{eff}}=\kappa\sqrt{-\ln\max_{\mu\in A,\nu\in
B}|S_{\mu\nu}|},\quad\kappa=2.0\mbox{ \AA}.$ (17)

Herein $\kappa$ is chosen such that $R_{AB}^{\mathrm{eff}}$ is similar to the
real-space distance between the atoms $A$ and $B$ (in Angstrom) when typical
double-zeta basis sets are used. We then choose the buffer atoms via Eq. (16),
but using $R_{AB}^{\mathrm{eff}}$ in Eq. (17) in place of $R_{AB}$. The
threshold radius at macroiteration 0, $R_{L}^{(0)}$, is chosen to be 2.0
Åirrespective of the subsystem $L$. However, $R_{L}^{(\tilde{m})}$ will become
subsystem-dependent when $\tilde{m}>0$ (vide infra).

  2. (2)

In case that the use of PHO boundary atoms to saturate the dangling bonds is
inappropriate, the corresponding PHO boundary atoms are redesignated as buffer
atoms, and the search for dangling bonds is performed again. This happens
under four circumstances:

    1. (a)

The PHO boundary atom is not a $p$-block element atom (see Fig. 5(a)). This is
because the PHO algorithm involves the hybridization of the valence $s$ and
$p$ orbitals to yield 4 hybrid orbitals, which makes sense only if the PHO
boundary atoms are $p$-block element atoms.

    2. (b)

The dangling bond is a multiple bond (see Fig. 5(b)). This is because the PHO
algorithm works only for the $\sigma$ type of active hybrid orbitals.

    3. (c)

The same PHO boundary atom is used to saturate more than one dangling bond
(Fig. 5(c)). While the original PHO method108 supports PHO boundary atoms with
more than one active hybrid orbitals, the present implementation is confined
to PHO boundary atoms with a single active hybrid orbital.

    4. (d)

All the environment atoms of a PHO boundary atom are hydrogen atoms (Fig.
5(d)). This usually happens when the PHO boundary atom is a methyl carbon
atom. In this case, it is a better strategy to convert the carbon atom to a
buffer atom and include the three methyl hydrogen atoms in the buffer, since
this hardly increases the computational cost but makes the description of the
methyl group much more reliable. Note that this includes PHO boundary atoms
that have zero environment atoms (e.g. halogens) as a special case.

Figure 5: Special treatment of dangling bonds: (a) the PHO boundary atom is
not a $p$-block element; (b) the dangling bond is not a single bond; (c) the
PHO boundary atom is connected to more than one buffer atoms; (d) all
environment atoms bonded to a PHO boundary atom are hydrogens. The central and
buffer atoms are indicated in black, while the PHO boundary atoms are
indicated in blue.

##  Appendix B Modified PHO

The PHO method developed by Gao et al.108 is a convenient way for saturating
dangling bonds in QM/MM calculations. Compared with the use of link atoms, the
PHO method suffers less from artifacts such as over-polarization of the QM
region by the MM region. Herein we introduce some modifications to the PHO
approach to describe more complex chemical environments.

Suppose that the global system is fragmented at a C-C bond (Fig. 6).
Traditionally, the free valence of the CQ atom is saturated by a hydrogen
atom; the atoms on the left side of the bond are thus assigned to subsystem
$L$, while those on the right side are termed environment. In the PHO method,
however, atom CB is included in subsystem $L$ and is called PHO boundary atom.
The four valence orbitals (one $2s$ and three $2p$ orbitals) of CB are
hybridized to form four hybrid orbitals (four $sp^{3}$ orbitals in Fig. 6(a);
three $sp^{2}$ orbitals and one $p$ orbital in Fig. 6(b)). The one that is
oriented towards CQ is termed active, whereas the other three are termed
auxiliary. The former is allowed to mix with the AOs of the remaining parts of
the subsystem, whereas the latter represent free valences of CB. Since they
are actually saturated in the whole system by bonding with the Rx groups, the
auxiliary hybrid orbitals must be frozen, i.e., projected out of the
variational space. Otherwise one is effectively calculating a carbon
triradical R1R2R3C-C$\vdots$, instead of the closed-shell system
R1R2R3C-CR4R5R6 or R1R2R3C-CR4R5. Nevertheless, the auxiliary hybrid orbitals
can accommodate electrons, so as to provide an effective external potential
for the rest of the subsystem.

Figure 6: The boundary of two subsystems capped by PHO boundary atoms (CB):
(a) the dangling bond is C(sp3)-C(sp3); (b) the dangling bond is
C(sp3)-C(sp2). The active and auxiliary hybrid orbitals are colored in red and
grey, respectively.

To uniquely define our modified PHO method, several details need to be
addressed:

  1. (1)

Generation of valence atomic orbitals
While in minimal basis set calculations some AOs can be unambiguously assigned
as the atomic valence orbitals, this is generally not true for extended basis
sets. In the original PHO method108, the STO-3G $2s$ and $2p$ are taken to be
the atomic valence orbitals, even when the calculation is performed with
extended basis sets. Although the deficiency of the STO-3G basis is partly
remedied by transforming the density matrix elements of the PHO boundary atom
from STO-3G to the actually employed basis, this algorithm can still be
potentially unreliable when the STO-3G orbitals are quite different from the
exact atomic valence orbitals (e.g., when the PHO boundary atom is a heavy
element). To resolve this issue, we perform a spherically averaged and
unpolarized atomic ROHF/ROKS calculation for each PHO boundary atom with the
actual basis, so that the valence orbitals of the atomic calculation form a
minimal valence basis which is of much better quality than STO-3G. We choose
only $p$-block atoms as PHO boundary atoms.

  2. (2)

Generation of hybrid orbitals
Hybridizing the valence $s$ and $p$ orbitals is trivial if the boundary atom
adopts a perfect tetrahedral ($sp^{3}$), trigonal ($sp^{2}$), or linear ($sp$)
geometry. However, real-world molecules seldom possess such exact symmetries.
In particular, if the PHO boundary atom has a lone pair, the deviations from
exact symmetries can be quite significant. Before discussing the generation of
hybrid orbitals under nonideal symmetries, we first introduce the concept of a
_hybrid orbital pointing towards_ an atom. A hybrid orbital of atom $B$ of the
following form,

$\psi^{B}=c_{s}^{B}\chi_{s}^{B}+c_{p_{x}}^{B}\chi_{p_{x}}^{B}+c_{p_{y}}^{B}\chi_{p_{y}}^{B}+c_{p_{z}}^{B}\chi_{p_{z}}^{B},$
(18)

is said to be pointing towards atom $A$, if its coefficients satisfy
$(c_{p_{x}}^{B},c_{p_{y}}^{B},c_{p_{z}}^{B})=k\vec{r}_{BA}$, with
$\vec{r}_{BA}$ being the vector from atom B to atom A and $k$ being a positive
constant of proportionality. In order for the definition to make sense, the
phases of the atomic valence orbials
$\\{\chi_{s}^{B},\chi_{p_{x}}^{B},\chi_{p_{y}}^{B},\chi_{p_{z}}^{B}\\}$ must
be such that when a hybrid orbital of atom $B$ is pointing towards atom $A$,
the lobe of that hybrid orbital to the direction of $A$ is larger than the
lobe pointing away from $A$. Intuitively, when a hybrid orbital of atom $B$
points towards atom $A$, its $\sigma$-overlap with the AOs of atom $A$ is
maximized, which is the most energetically favorable scenario for bonding.
Thus, for a $p$-block element atom $B$ bonded to $n$ atoms, we expect that $n$
hybrid orbitals of the atom will point (at least approximately) towards the
atoms directly bonded to $B$, and the remaining $4-n$ hybrid orbitals are
either non-bonding orbitals or form $\pi$ bonds with the neighboring atoms.
While for $n\leq 3$ the former hybrid orbitals can indeed be made exactly
pointing towards the neighboring atoms, for $n=4$ this is generally not
possible because the problem is overdetermined: requiring all four hybrid
orbitals to point towards specified atoms introduces 8 constraints on the AO
coefficients of the hybrid orbitals (2 for each hybrid orbital), while the
orthonormality conditions of the four hybrid orbitals introduce another
$4\times(4+1)/2=10$ constraints, giving a total of 18 constraints, but there
are only 16 variables to be determined (4 for each hybrid orbital). This is to
be contrasted with the $n=3$ case, where there are exactly 16 constraints,
giving a well-determined problem. The $n=2$ case is accordingly
underdetermined, which can also be problematic as the solution to the hybrid
orbitals is not unique.

Gao’s approach108, 120 for $n=4$ amounts to requiring three of the hybrid
orbitals to point towards the three environment atoms R4, R5, and R6, so that
the fourth hybrid orbital does not necessarily point exactly towards the QM
atom CQ (Fig. 6). Thus, one of the “exactly pointing” constraints is relaxed,
making the problem well-determined. It was found by Exner and coworkers121
that this will deteriorate the description of the CQ-CB bond, especially if
the local geometry of the boundary atom CB is significantly distorted from
ideal tetrahedral geometry. Therefore, they proposed to force three of the
hybrid orbitals to point towards CQ, R4 and R5, at the expense that the fourth
hybrid orbital will not point exactly towards R6. This approach improves the
accuracy for some molecules, but the price to pay is that the result is no
longer invariant with respect to permutations of the environment atoms R4, R5,
and R6, which is clearly undesirable. They then proposed to fix the latter
problem by generating a set of hybrid orbitals for each of these permutations,
and average the resulting hybrid orbitals, followed by a symmetric
orthogonalization, which solves the aforementioned problem at the expense of
increased complexity of the method. Herein, we propose a simple and elegant
remedy to this problem by generating four hybrid orbitals
$\\{\tilde{\psi}^{B}\\}$ that are normalized but not necessarily mutually
orthogonal:

$\tilde{\psi}^{B}=\tilde{c}_{s}^{B}\chi_{s}^{B}+\tilde{c}_{p_{x}}^{B}\chi_{p_{x}}^{B}+\tilde{c}_{p_{y}}^{B}\chi_{p_{y}}^{B}+\tilde{c}_{p_{z}}^{B}\chi_{p_{z}}^{B},$
$\tilde{c}_{s}^{2}=1/4,\quad\tilde{c}_{p_{x}}^{2}+\tilde{c}_{p_{y}}^{2}+\tilde{c}_{p_{z}}^{2}=3/4.$
(19)

To be unique, $\\{\tilde{\psi}^{B}\\}$ are required to be $sp^{3}$-hybridized,
as implied by the latter conditions in Eq. (19). Being nonorthogonal,
$\\{\tilde{\psi}^{B}\\}$ can be made to point exactly towards the respective
atoms. They are then symmetrically orthogonalized to yield the final
orthonormal hybrid orbitals. Since the symmetric orthogonalization procedure
produces orbitals that are closest to the original ones, the present hybrid
orbitals are more realistic than the previous ones108, 120, 121, especially
for heavily distorted structures (which will be demonstrated in a separate
paper).

When $n=3$ or $n=2$, we still generate the $n$ hybrid orbitals
$\\{\psi^{B}\\}$ by the above procedure (i.e., the initial hybrid orbitals
$\\{\tilde{\psi}^{B}\\}$ are still taken as $sp^{3}$-hybridized), and find the
remaining $4-n$ hybrid orbitals by Schmidt orthogonalization of the atomic
valence $s$ and $p$ orbitals against the $n$ hybrid orbitals $\\{\psi^{B}\\}$.
Note that this works even when the real hybrid orbitals are not
$sp^{3}$-hybridized, because the symmetric orthogonalization step
automatically adjusts the $s$ components of the hybrid orbitals. For example,
for a trigonal planar PHO boundary atom like the CB atom in Fig. 6(b), we
generate three nonorthonormal $sp^{3}$ hybrid orbitals
$\\{\tilde{\psi}^{B}\\}$ that point towards the atoms CQ, R4, and R5,
respectively. If the geometry of CB is exactly threefold symmetric, each
hybrid orbital must have exactly $1/3$ of $s$ character to be mutually
orthogonal, thereby giving the desired three $sp^{2}$ hybrid orbitals. The
Schmidt orthogonalization step then finds the $p$ orbital perpendicular to the
trigonal plane as the fourth hybrid orbital, since this is the orthogonal
complement of the aforementioned three $sp^{2}$ hybrid orbitals in the four-
dimensional space of the atomic valence orbitals of CB.

  3. (3)

Projection of auxiliary hybrid orbitals out of variational space
This is required under two circumstances, one is the projection of the
components of the auxiliary hybrid orbitals from the pFLMOs before the
elimination of linear dependence, to ensure that the linearly independent
pFLMOs are free of contributions from the auxiliary hybrid orbitals, and the
other is the projection of the contributions of the auxiliary hybrid orbitals
from the Fock matrix before (block-)diagonalization (Steps (5) and (7e)), to
ensure that the new canonical (for macroiteration 0) or localized (for
macroiteration $\tilde{m}>0$) MOs from the (block-)diagonalization of the Fock
matrix are free of contributions from the auxiliary hybrid orbitals. The
projection is trivial when there is only one PHO boundary atom in the system,
where the auxiliary hybrid orbitals are mutually orthogonal, as described by
Gao108. When there are more than one PHO boundary atoms, the auxiliary hybrid
orbitals that belong to different PHO boundary atoms are generally not exactly
orthogonal with each other (although usually very nearly so, since the PHO
boundary atoms of a subsystem are usually relatively far from each other). In
this case, they must be first symmetrically orthogonalized to give the
orthonormal auxiliary hybrid orbitals $\\{\psi_{\mu}^{B,\mathrm{PHO}}\\}$.
Then, we project out the contributions of $\\{\psi_{\mu}^{B,\mathrm{PHO}}\\}$
from the pFLMOs and/or the Fock matrix, as usually done122.

  4. (4)

Occupation assignment of auxiliary hybrid orbitals
The original PHO method108, which was developed in the context of QM/MM,
requires that the atomic charge of the PHO boundary atom $q_{B}$ be known
before the SCF iterations, which is generally available from the definition of
the MM force field. Then, if $n=4$ (i.e., the number of auxiliary hybrid
orbitals is 3), the occupation numbers are determined as

$f(\psi_{\mu}^{B,\mathrm{PHO}})=1-q_{B}/3,$ (20)

so that if the bond CQ-CB is exactly non-polar, the QM charge of the atom CB
reproduces the MM charge $q_{B}$. The description of the charge of the PHO
boundary atom by fractional occupations of the auxiliary hybrid orbitals,
rather than by point charges, is the key to avoid over-polarization effects at
the boundary of the QM region108. However, in macroiteration 0 of the iOI
algorithm, the charge $q_{B}$ is not known in advance and is set simply to
zero. For macroiteration $\tilde{m}>0$, we identify the subsystem at
macroiteration $\tilde{m}-1$ whose pFrag contains the PHO boundary atom in
question. The Löwdin charge of CB in that subsystem is then assigned to
$q_{B}$.

For $n=3$ or $n=2$, it is generally advocated to partition the occupation
numbers unevenly between the hybrid orbitals when some of the auxiliary hybrid
orbitals are lone pairs, e.g., the occupation numbers are partitioned in a
2:1:1 ratio among the three auxiliary hybrid orbitals for amine nitrogen atom,
with the first hybrid orbital being the lone pair orbital. The ratio can be
even more delicately tuned by considering the electronegativities of the
neighboring atoms121, 122. However there are many situations where the ratio
of occupation numbers among the auxiliary hybrid orbitals cannot be
unambiguously assigned, such as when the PHO boundary atom is an amide
nitrogen atom, where the occupation of the non-bonding auxiliary hybrid
orbital is difficult to assign accurately without inventing a rule
specifically for amides. Therefore, for simplicity, we always partition the
fractional occupations of a PHO boundary orbital evenly across all auxiliary
hybrid orbitals regardless of the nature of the PHO boundary atom. As will be
demonstrated in a separate paper, the results are not noticeably inferior to
those schemes where the occupation numbers are partitioned unevenly, at least
when only the quality of the LMOs are concerned.

  5. (5)

Coulomb interaction between the subsystem and environment
In the original PHO implementation108 one generally describes all MM atoms as
point charges, which can incur a non-trivial cost when the MM region is very
large. In the iOI algorithm, however, it is only necessary to accurately
describe the polarization effect of the point charges on the subsystem, but
not the electrostatic energy between the point charges and the subsystem,
since the electronic energy is ultimately obtained by a global SCF
calculation. As the polarization effect is much more short-ranged ($R^{-4}$)
than the electrostatic energy ($R^{-1}$), we find it only necessary to include
the point charges of the environment atoms that are directly bonded to the PHO
boundary atoms. As before, the point charges are zero for macroiteration 0 and
are taken as the Löwdin charges from the previous macroiteration for
macroiteration $\tilde{m}>0$.

##  Appendix C Merging subsystems

As described in Step (7b) of the iOI algorithm, in the end of each
macroiteration, the pFrags will be merged together to give new pFrags, in a
way that every new pFrag consists of either one or two old pFrags, which is
followed by a capping step. The merging process only makes sense if the two
old pFrags comprising a new one are spatially adjacent. This can be a non-
trivial condition to fulfill for molecules with complex topologies. Consider,
e.g., a star-shaped molecule where a pFrag $\mathcal{F}_{L}$ is connected to
three mutually disconnected pFrags $\mathcal{F}_{M_{1}}$,
$\mathcal{F}_{M_{2}}$ and $\mathcal{F}_{M_{3}}$. If the three pFrags are
sufficiently similar but not identical, then it is not easy to decide which of
them should be merged with $\mathcal{F}_{L}$. Choosing the spatially closest
one to merge with $\mathcal{F}_{L}$ is not always the best solution for the
remaining two pFrags may be too far from each other. Keeping such situations
in mind, we propose the following merging process:

  1. (1)

Make a list of the pFrags of all not yet converged subsystems, the ordering of
which is immaterial;

  2. (2)

If the number of nonconvergent subsystems is odd, there will be at least one
pFrag that cannot be merged with others. To make the sizes of the new pFrags
as uniform as possible, the pFrag with the largest number of basis functions
is removed from the list, such that it is directly carried over to the next
macroiteration as a new pFrag. The remaining pFrags can then be merged in a
pairwise fashion;

  3. (3)

Calculate the distance
$d(\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{M}^{(\tilde{m})})$ between ever
pair of pFrags, $\mathcal{F}_{L}^{(\tilde{m})}$ and
$\mathcal{F}_{M}^{(\tilde{m})}$, which is defined as the minimum of the
effective distances, Eq. (17), between atoms in
$\mathcal{F}_{L}^{(\tilde{m})}$ and in $\mathcal{F}_{M}^{(\tilde{m})}$;

  4. (4)

Starting from the first pFrag in the list, form a pFrag pair with pFrag
$\mathcal{F}_{L^{\prime}}^{(\tilde{m})}$ and its spatially closest pFrag
$\mathcal{F}_{L}^{(\tilde{m})}$. Remove $\mathcal{F}_{L}^{(\tilde{m})}$ and
$\mathcal{F}_{L^{\prime}}^{(\tilde{m})}$ from the list and repeat the process
for the next surviving entry of the list;

  5. (5)

For each _pair of pFrag pairs_
$\\{\\{\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{L^{\prime}}^{(\tilde{m})}\\},\\{\mathcal{F}_{M}^{(\tilde{m})},\mathcal{F}_{M^{\prime}}^{(\tilde{m})}\\}\\}$,
permute the four pFrags so as to reduce the maximum intrapair distance
$\max(d(\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{L^{\prime}}^{(\tilde{m})}),d(\mathcal{F}_{M}^{(\tilde{m})},\mathcal{F}_{M^{\prime}}^{(\tilde{m})}))$.
This is repeated for every possible pair of pFrag pairs, until the maximum
intrapair distance cannot be reduced by permutations. There is however an
exception: if the minimum value for the maximum intrapair distance is larger
than a threshold, $d_{\mathrm{max}}=4.0$ Å, no more than one of the two pFrag
pairs will actually give rise to a new pFrag (see next step). In this case we
minimize
$\min(d(\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{L^{\prime}}^{(\tilde{m})}),d(\mathcal{F}_{M}^{(\tilde{m})},\mathcal{F}_{M^{\prime}}^{(\tilde{m})}))$
instead of
$\max(d(\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{L^{\prime}}^{(\tilde{m})}),d(\mathcal{F}_{M}^{(\tilde{m})},\mathcal{F}_{M^{\prime}}^{(\tilde{m})}))$.
The pFrag pairs are then as mutually proximal as possible. Another benefit of
this step is that it minimizes (though not necessarily eliminates) the
dependence of the result on the ordering of the pFrags in the list;

  6. (6)

For each pFrag pair
$\\{\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{M}^{(\tilde{m})}\\}$, if
$d(\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{M}^{(\tilde{m})})\leq
d_{\mathrm{max}}$, merge the two pFrags to form a new pFrag. If
$d(\mathcal{F}_{L}^{(\tilde{m})},\mathcal{F}_{M}^{(\tilde{m})})>d_{\mathrm{max}}$,
the two pFrags are not considered to be adjacent and are therefore carried
over directly to the next macroiteration as two new pFrags;

  7. (7)

Generate the caps $\mathcal{C}_{L}^{(\tilde{m}+1)}$ for all the new pFrags
$\mathcal{F}_{L}^{(\tilde{m}+1)}$, according to Eq. (16). The buffer radii,
$R_{L}^{(\tilde{m}+1)}$, are calculated as follows: (a) set
$R_{L}^{(\tilde{m}+1)}=\max_{I}R_{\mathcal{P}_{L}^{(\tilde{m}+1)}(I)}^{(\tilde{m})}$,
i.e., the maximum of the buffer radii of the parental subsystems of
$\mathcal{S}_{L}^{(\tilde{m}+1)}$; (b) increase $R_{L}^{(\tilde{m}+1)}$ until
exactly one additional atom is included in the buffer; (c) increase
$R_{L}^{(\tilde{m}+1)}$ by 1.0 Å. Step (b) guarantees that the new subsystem
$\mathcal{S}_{L}^{(\tilde{m}+1)}$ is strictly larger than its parental
subsystems, while Step (c) ensures that the incremental cap
$\Delta\mathcal{C}_{L}^{(\tilde{m}+1)}$ is not too small.

## References

  * Li et al. 2014 Li, Z.; Li, H.; Suo, B.; Liu, W. Localization of molecular orbitals: from fragments to molecule. _Acc. Chem. Res._ 2014, _47_ , 2758–2767
  * Høyvik and Jørgensen 2016 Høyvik, I.-M.; Jørgensen, P. Characterization and generation of local occupied and virtual Hartree-Fock orbitals. _Chem. Rev._ 2016, _116_ , 3306–3327
  * Prodan and Kohn 2005 Prodan, E.; Kohn, W. Nearsightedness of electronic matter. _Proc. Nat. Acad. Sci._ 2005, _102_ , 11635–11638
  * Warshel and Karplus 1972 Warshel, A.; Karplus, M. Calculation of ground and excited state potential surfaces of conjugated molecules. I. Formulation and parametrization. _J. Am. Chem. Soc._ 1972, _94_ , 5612–5625
  * Day et al. 1996 Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An effective fragment method for modeling solvent effects in quantum mechanical calculations. _J. Chem. Phys._ 1996, _105_ , 1968–1986
  * Gordon et al. 2001 Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. The effective fragment potential method: A QM-based MM approach to modeling environmental effects in chemistry. _J. Phys. Chem. A_ 2001, _105_ , 293–307
  * Xie and Gao 2007 Xie, W.; Gao, J. Design of a next generation force field: the X-POL potential. _J. Chem. Theory Comput._ 2007, _3_ , 1890–1900
  * Dapprich et al. 1999 Dapprich, S.; Komáromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. A new ONIOM implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives. _J. Mol. Struct. (THEOCHEM)_ 1999, _461_ , 1–21
  * Vreven et al. 2006 Vreven, T.; Byun, K. S.; Komáromi, I.; Dapprich, S.; Montgomery Jr, J. A.; Morokuma, K.; Frisch, M. J. Combining quantum mechanics methods with molecular mechanics methods in ONIOM. _J. Chem. Theory Comput._ 2006, _2_ , 815–826
  * Guo et al. 2010 Guo, W.; Wu, A.; Xu, X. XO: An extended ONIOM method for accurate and efficient geometry optimization of large molecules. _Chem. Phys. Lett._ 2010, _498_ , 203–208
  * Wang et al. 2012 Wang, Y.; Sosa, C. P.; Cembran, A.; Truhlar, D. G.; Gao, J. Multilevel X-Pol: A fragment-based method with mixed quantum mechanical representations of different fragments. _J. Phys. Chem. B_ 2012, _116_ , 6781–6788
  * Yang 1991 Yang, W. Direct calculation of electron density in density-functional theory. _Phys. Rev. Lett._ 1991, _66_ , 1438
  * Yang and Lee 1995 Yang, W.; Lee, T.-S. A density-matrix divide-and-conquer approach for electronic structure calculations of large molecules. _J. Chem. Phys._ 1995, _103_ , 5674–5678
  * Akama et al. 2007 Akama, T.; Kobayashi, M.; Nakai, H. Implementation of divide-and-conquer method including Hartree-Fock exchange interaction. _J. Comput. Chem._ 2007, _28_ , 2003–2012
  * He and Merz Jr 2010 He, X.; Merz Jr, K. M. Divide and conquer Hartree- Fock calculations on proteins. _J. Chem. Theory Comput._ 2010, _6_ , 405–411
  * Wesolowski and Warshel 1993 Wesolowski, T. A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. _J. Phys. Chem._ 1993, _97_ , 8050–8053
  * Assfeld and Rivail 1996 Assfeld, X.; Rivail, J.-L. Quantum chemical computations on parts of large molecules: the ab initio local self consistent field method. _Chem. Phys. Lett._ 1996, _263_ , 100–106
  * Iannuzzi et al. 2006 Iannuzzi, M.; Kirchner, B.; Hutter, J. Density functional embedding for molecular systems. _Chem. Phys. Lett._ 2006, _421_ , 16–20
  * Fux et al. 2010 Fux, S.; Jacob, C. R.; Neugebauer, J.; Visscher, L.; Reiher, M. Accurate frozen-density embedding potentials as a first step towards a subsystem description of covalent bonds. _J. Chem. Phys._ 2010, _132_ , 164101
  * Chulhai and Jensen 2015 Chulhai, D. V.; Jensen, L. Frozen Density Embedding with External Orthogonality in Delocalized Covalent Systems. _J. Chem. Theory Comput._ 2015, _11_ , 3080–3088
  * Huang et al. 2011 Huang, C.; Pavone, M.; Carter, E. A. Quantum mechanical embedding theory based on a unique embedding potential. _J. Chem. Phys._ 2011, _134_ , 154110
  * Manby et al. 2012 Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller III, T. F. A simple, exact density-functional-theory embedding scheme. _J. Chem. Theory Comput._ 2012, _8_ , 2564–2568
  * Knizia and Chan 2012 Knizia, G.; Chan, G. K.-L. Density matrix embedding: A simple alternative to dynamical mean-field theory. _Phys. Rev. Lett._ 2012, _109_ , 186404
  * Knizia and Chan 2013 Knizia, G.; Chan, G. K.-L. Density matrix embedding: A strong-coupling quantum embedding theory. _J. Chem. Theory Comput._ 2013, _9_ , 1428–1432
  * Tamukong et al. 2014 Tamukong, P. K.; Khait, Y. G.; Hoffmann, M. R. Density differences in embedding theory with external orbital orthogonality. _J. Phys. Chem. A_ 2014, _118_ , 9182–9200
  * Welborn et al. 2016 Welborn, M.; Tsuchimochi, T.; Van Voorhis, T. Bootstrap embedding: An internally consistent fragment-based method. _J. Chem. Phys._ 2016, _145_ , 074102
  * Hégely et al. 2016 Hégely, B.; Nagy, P. R.; Ferenczy, G. G.; Kállay, M. Exact density functional and wave function embedding schemes based on orbital localization. _J. Chem. Phys._ 2016, _145_ , 064107
  * Culpitt et al. 2017 Culpitt, T.; Brorsen, K. R.; Hammes-Schiffer, S. Communication: Density functional theory embedding with the orthogonality constrained basis set expansion procedure. _J. Chem. Phys._ 2017, _146_ , 211101
  * Chen et al. 2019 Chen, M.; Baer, R.; Neuhauser, D.; Rabani, E. Overlapped embedded fragment stochastic density functional theory for covalently-bonded materials. _J. Chem. Phys._ 2019, _150_ , 034106
  * Niffenegger et al. 2019 Niffenegger, K.; Oueis, Y.; Nafziger, J.; Wasserman, A. Density embedding with constrained chemical potential. _Mol. Phys._ 2019, _117_ , 2188–2194
  * Macetti et al. 2021 Macetti, G.; Wieduwilt, E. K.; Genoni, A. QM/ELMO: A Multi-Purpose Fully Quantum Mechanical Embedding Scheme Based on Extremely Localized Molecular Orbitals. _J. Phys. Chem. A_ 2021, _125_ , 2709–2726
  * Stoll et al. 1980 Stoll, H.; Wagenblast, G.; Preu$\beta$, H. On the use of local basis sets for localized molecular orbitals. _Theor. Chim. Acta_ 1980, _57_ , 169–178
  * Stoll 1992 Stoll, H. The correlation energy of crystalline silicon. _Chem. Phys. Lett._ 1992, _191_ , 548–552
  * Friedrich et al. 2007 Friedrich, J.; Hanrath, M.; Dolg, M. Fully automated implementation of the incremental scheme: Application to CCSD energies for hydrocarbons and transition metal compounds. _J. Chem. Phys._ 2007, _126_ , 154110
  * Li et al. 2002 Li, S.; Ma, J.; Jiang, Y. Linear scaling local correlation approach for solving the coupled cluster equations of large systems. _J. Comput. Chem._ 2002, _23_ , 237–244
  * Li and Piecuch 2010 Li, W.; Piecuch, P. Improved design of orbital domains within the cluster-in-molecule local correlation framework: Single-environment cluster-in-molecule ansatz and its application to local coupled-cluster approach with singles and doubles. _J. Phys. Chem. A_ 2010, _114_ , 8644–8657
  * Flocke and Bartlett 2004 Flocke, N.; Bartlett, R. J. A natural linear scaling coupled-cluster method. _J. Chem. Phys._ 2004, _121_ , 10935–10944
  * Hughes et al. 2008 Hughes, T. F.; Flocke, N.; Bartlett, R. J. Natural linear-scaled coupled-cluster theory with local transferable triple excitations: Applications to peptides. _J. Phys. Chem. A_ 2008, _112_ , 5994–6003
  * Ziółkowski et al. 2010 Ziółkowski, M.; Jansik, B.; Kjærgaard, T.; Jørgensen, P. Linear scaling coupled cluster method with correlation energy based error control. _J. Chem. Phys._ 2010, _133_ , 014107
  * Jacobson and Herbert 2011 Jacobson, L. D.; Herbert, J. M. An efficient, fragment-based electronic structure method for molecular systems: Self-consistent polarization with perturbative two-body exchange and dispersion. _J. Chem. Phys._ 2011, _134_ , 094118
  * Aoki and Gu 2012 Aoki, Y.; Gu, F. L. An elongation method for large systems toward bio-systems. _Phys. Chem. Chem. Phys._ 2012, _14_ , 7640–7668
  * Gadre et al. 1994 Gadre, S. R.; Shirsat, R. N.; Limaye, A. C. Molecular tailoring approach for simulation of electrostatic properties. _J. Phys. Chem._ 1994, _98_ , 9165–9169
  * Ganesh et al. 2006 Ganesh, V.; Dongare, R. K.; Balanarayan, P.; Gadre, S. R. Molecular tailoring approach for geometry optimization of large molecules: Energy evaluation and parallelization strategies. _J. Chem. Phys._ 2006, _125_ , 104109
  * Kitaura et al. 1999 Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment molecular orbital method: an approximate computational method for large molecules. _Chem. Phys. Lett._ 1999, _313_ , 701–706
  * Huang et al. 2005 Huang, L.; Massa, L.; Karle, J. Kernel energy method illustrated with peptides. _Int. J. Quantum Chem._ 2005, _103_ , 808–817
  * Zhang and Zhang 2003 Zhang, D. W.; Zhang, J. Z. H. Molecular fractionation with conjugate caps for full quantum mechanical calculation of protein–molecule interaction energy. _J. Chem. Phys._ 2003, _119_ , 3599–3605
  * Xu et al. 2019 Xu, M.; He, X.; Zhu, T.; Zhang, J. Z. H. A fragment quantum mechanical method for metalloproteins. _J. Chem. Theory Comput._ 2019, _15_ , 1430–1439
  * Deev and Collins 2005 Deev, V.; Collins, M. A. Approximate ab initio energies by systematic molecular fragmentation. _J. Chem. Phys._ 2005, _122_ , 154102
  * Collins 2012 Collins, M. A. Systematic fragmentation of large molecules by annihilation. _Phys. Chem. Chem. Phys._ 2012, _14_ , 7744–7751
  * Jiang et al. 2006 Jiang, N.; Ma, J.; Jiang, Y. Electrostatic field-adapted molecular fractionation with conjugated caps for energy calculations of charged biomolecules. _J. Chem. Phys._ 2006, _124_ , 114112
  * Li et al. 2007 Li, W.; Li, S.; Jiang, Y. Generalized energy-based fragmentation approach for computing the ground-state energies and properties of large molecules. _J. Phys. Chem. A_ 2007, _111_ , 2193–2199
  * Hua et al. 2010 Hua, S.; Hua, W.; Li, S. An efficient implementation of the generalized energy-based fragmentation approach for general large molecules. _J. Phys. Chem. A_ 2010, _114_ , 8126–8134
  * Dahlke and Truhlar 2007 Dahlke, E. E.; Truhlar, D. G. Electrostatically embedded many-body expansion for large systems, with applications to water clusters. _J. Chem. Theory Comput._ 2007, _3_ , 46–53
  * Dahlke and Truhlar 2007 Dahlke, E. E.; Truhlar, D. G. Electrostatically Embedded Many-Body Correlation Energy, with Applications to the Calculation of Accurate Second-Order Møller- Plesset Perturbation Theory Energies for Large Water Clusters. _J. Chem. Theory Comput._ 2007, _3_ , 1342–1348
  * Řezáč and Salahub 2010 Řezáč, J.; Salahub, D. R. Multilevel fragment-based approach (MFBA): A novel hybrid computational method for the study of large molecules. _J. Chem. Theory Comput._ 2010, _6_ , 91–99
  * Beran and Nanda 2010 Beran, G. J.; Nanda, K. Predicting organic crystal lattice energies with chemical accuracy. _J. Phys. Chem. Lett._ 2010, _1_ , 3480–3487
  * Mayhall and Raghavachari 2011 Mayhall, N. J.; Raghavachari, K. Molecules-in-molecules: An extrapolated fragment-based approach for accurate calculations on large molecules and materials. _J. Chem. Theory Comput._ 2011, _7_ , 1336–1343
  * Mayhall and Raghavachari 2012 Mayhall, N. J.; Raghavachari, K. Many-overlapping-body (MOB) expansion: A generalized many body expansion for nondisjoint monomers in molecular fragmentation calculations of covalent molecules. _J. Chem. Theory Comput._ 2012, _8_ , 2669–2675
  * Richard and Herbert 2012 Richard, R. M.; Herbert, J. M. A generalized many-body expansion and a unified view of fragment-based methods in electronic structure theory. _J. Chem. Phys._ 2012, _137_ , 064113
  * Le et al. 2012 Le, H.-A.; Tan, H.-J.; Ouyang, J. F.; Bettens, R. P. Combined fragmentation method: A simple method for fragmentation of large molecules. _J. Chem. Theory Comput._ 2012, _8_ , 469–478
  * Saha and Raghavachari 2014 Saha, A.; Raghavachari, K. Dimers of dimers (DOD): A new fragment-based method applied to large water clusters. _J. Chem. Theory Comput._ 2014, _10_ , 58–67
  * Gordon et al. 2012 Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation methods: A route to accurate calculations on large systems. _Chem. Rev._ 2012, _112_ , 632–672
  * He et al. 2014 He, X.; Zhu, T.; Wang, X.; Liu, J.; Zhang, J. Z. Fragment quantum mechanical calculation of proteins and its applications. _Acc. Chem. Res._ 2014, _47_ , 2748–2757
  * Li et al. 2014 Li, S.; Li, W.; Ma, J. Generalized energy-based fragmentation approach and its applications to macromolecules and molecular aggregates. _Acc. Chem. Res._ 2014, _47_ , 2712–2720
  * Collins et al. 2014 Collins, M. A.; Cvitkovic, M. W.; Bettens, R. P. The combined fragmentation and systematic molecular fragmentation methods. _Acc. Chem. Res._ 2014, _47_ , 2776–2785
  * Tanaka et al. 2014 Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. Electron-correlated fragment-molecular-orbital calculations for biomolecular and nano systems. _Phys. Chem. Chem. Phys._ 2014, _16_ , 10310–10344
  * Gadre et al. 2014 Gadre, S. R.; Yeole, S. D.; Sahu, N. Quantum chemical investigations on molecular clusters. _Chem. Rev._ 2014, _114_ , 12132–12173
  * Collins and Bettens 2015 Collins, M. A.; Bettens, R. P. Energy-based molecular fragmentation methods. _Chem. Rev._ 2015, _115_ , 5607–5642
  * Raghavachari and Saha 2015 Raghavachari, K.; Saha, A. Accurate composite and fragment-based quantum chemical models for large molecules. _Chem. Rev._ 2015, _115_ , 5643–5677
  * Gordon 2017 Gordon, M. S. _Fragmentation: Toward Accurate Calculations on Complex Molecular Systems_ ; John Wiley & Sons, 2017
  * Saha and Raghavachari 2015 Saha, A.; Raghavachari, K. Analysis of different fragmentation strategies on a variety of large peptides: Implementation of a low level of theory in fragment-based methods can be a crucial factor. _J. Chem. Theory Comput._ 2015, _11_ , 2012–2023
  * Kumar and Iyengar 2019 Kumar, A.; Iyengar, S. S. Fragment-based electronic structure for potential energy surfaces using a superposition of fragmentation topologies. _J. Chem. Theory Comput._ 2019, _15_ , 5769–5786
  * Shang et al. 2010 Shang, H.; Xiang, H.; Li, Z.; Yang, J. Linear scaling electronic structure calculations with numerical atomic basis set. _Int. Rev. Phys. Chem._ 2010, _29_ , 665–691
  * Goedecker 1999 Goedecker, S. Linear scaling electronic structure methods. _Rev. Mod. Phys._ 1999, _71_ , 1085–1123
  * Bowler and Miyazaki 2012 Bowler, D. R.; Miyazaki, T. \mathcal$\\{$O$\\}$(N) methods in electronic structure calculations. _Reports on Progress in Physics_ 2012, _75_ , 036503
  * Kussmann et al. 2013 Kussmann, J.; Beer, M.; Ochsenfeld, C. Linear-scaling self-consistent field methods for large molecules. _WIREs Comput. Mol. Sci._ 2013, _3_ , 614–636
  * Chiba et al. 2007 Chiba, M.; Fedorov, D. G.; Kitaura, K. Time-dependent density functional theory based upon the fragment molecular orbital method. _J. Chem. Phys._ 2007, _127_ , 104108
  * Khait and Hoffmann 2010 Khait, Y. G.; Hoffmann, M. R. Embedding theory for excited states. _J. Chem. Phys._ 2010, _133_ , 044107
  * Daday et al. 2013 Daday, C.; König, C.; Valsson, O.; Neugebauer, J.; Filippi, C. State-specific embedding potentials for excitation-energy calculations. _J. Chem. Theory Comput._ 2013, _9_ , 2355–2367
  * Hedegård et al. 2013 Hedegård, E. D.; List, N. H.; Jensen, H. J. A.; Kongsted, J. The multi-configuration self-consistent field method within a polarizable embedded framework. _J. Chem. Phys._ 2013, _139_ , 044101
  * Yoshikawa et al. 2013 Yoshikawa, T.; Kobayashi, M.; Fujii, A.; Nakai, H. Novel approach to excited-state calculations of large molecules based on divide-and-conquer method: application to photoactive yellow protein. _J. Phys. Chem. B_ 2013, _117_ , 5565–5573
  * Daday et al. 2014 Daday, C.; König, C.; Neugebauer, J.; Filippi, C. Wavefunction in density functional theory embedding for excited states: Which wavefunctions, which densities? _ChemPhysChem_ 2014, _15_ , 3205–3217
  * Artiukhin et al. 2015 Artiukhin, D. G.; Jacob, C. R.; Neugebauer, J. Excitation energies from frozen-density embedding with accurate embedding potentials. _J. Chem. Phys._ 2015, _142_ , 234101
  * Gurunathan et al. 2016 Gurunathan, P. K.; Acharya, A.; Ghosh, D.; Kosenkov, D.; Kaliman, I.; Shao, Y.; Krylov, A. I.; Slipchenko, L. V. Extension of the effective fragment potential method to macromolecules. _J. Phys. Chem. B_ 2016, _120_ , 6562–6574
  * Li et al. 2016 Li, W.; Li, Y.; Lin, R.; Li, S. Generalized Energy-Based Fragmentation Approach for Localized Excited States of Large Systems. _J. Phys. Chem. A_ 2016, _120_ , 9667–9677
  * Liu et al. 2019 Liu, J.; Sun, H.; Glover, W. J.; He, X. Prediction of excited-state properties of oligoacene crystals using fragment-based quantum mechanical method. _J. Phys. Chem. A_ 2019, _123_ , 5407–5417
  * Chen et al. 2019 Chen, W.-K.; Fang, W.-H.; Cui, G. A multi-layer energy-based fragment method for excited states and nonadiabatic dynamics. _Phys. Chem. Chem. Phys._ 2019, _21_ , 22695–22699
  * Chen et al. 2019 Chen, W.-K.; Fang, W.-H.; Cui, G. Integrating Machine Learning with Multi-Layer Energy-Based Fragment Method for Excited States of Large Systems. _J. Phys. Chem. Lett._ 2019,
  * Macetti and Genoni 2020 Macetti, G.; Genoni, A. Quantum Mechanics/Extremely Localized Molecular Orbital Embedding Strategy for Excited States: Coupling to Time-Dependent Density Functional Theory and Equation-of-Motion Coupled Cluster. _J. Chem. Theory Comput._ 2020, _16_ , 7490–7506
  * Al Hajj et al. 2005 Al Hajj, M.; Malrieu, J.-P.; Guihéry, N. Renormalized excitonic method in terms of block excitations: Application to spin lattices. _Phys. Rev. B_ 2005, _72_ , 224412
  * Zhang et al. 2012 Zhang, H.; Malrieu, J.-P.; Ma, H.; Ma, J. Implementation of renormalized excitonic method at ab initio level. _J. Comput. Chem._ 2012, _33_ , 34–43
  * Ma et al. 2012 Ma, Y.; Liu, Y.; Ma, H. A new fragment-based approach for calculating electronic excitation energies of large systems. _J. Chem. Phys._ 2012, _136_ , 024113
  * Liu and Herbert 2015 Liu, J.; Herbert, J. M. An efficient and accurate approximation to time-dependent density functional theory for systems of weakly coupled monomers. _J. Chem. Phys._ 2015, _143_ , 034106
  * Liu and Herbert 2016 Liu, J.; Herbert, J. M. Local excitation approximations to time-dependent density functional theory for excitation energies in solution. _J. Chem. Theory Comput._ 2016, _12_ , 157–166
  * Closser et al. 2015 Closser, K. D.; Ge, Q.; Mao, Y.; Shao, Y.; Head-Gordon, M. Superposition of fragment excitations for excited states of large clusters with application to helium clusters. _J. Chem. Theory Comput._ 2015, _11_ , 5791–5803
  * Chulhai and Jensen 2016 Chulhai, D. V.; Jensen, L. External orthogonality in subsystem time-dependent density functional theory. _Phys. Chem. Chem. Phys._ 2016, _18_ , 21032–21039
  * Nakai and Yoshikawa 2017 Nakai, H.; Yoshikawa, T. Development of an excited-state calculation method for large systems using dynamical polarizability: A divide-and-conquer approach at the time-dependent density functional level. _J. Chem. Phys._ 2017, _146_ , 124123
  * Ge et al. 2017 Ge, Q.; Mao, Y.; White, A. F.; Epifanovsky, E.; Closser, K. D.; Head-Gordon, M. Simulating the absorption spectra of helium clusters (N= 70, 150, 231, 300) using a charge transfer correction to superposition of fragment single excitations. _J. Chem. Phys._ 2017, _146_ , 044111
  * Liu et al. 2014 Liu, J.; Zhang, Y.; Liu, W. Photoexcitation of Light-Harvesting C–P–C60 Triads: A FLMO-TD-DFT Study. _J. Chem. Theory Comput._ 2014, _10_ , 2436–2448
  * Liu and Hoffmann 2014 Liu, W.; Hoffmann, M. R. SDS: the ‘static-dynamic-static’ framework for strongly correlated electrons. _Theor. Chem. Acc._ 2014, _133_ , 1481
  * Liu and Hoffmann 2016 Liu, W.; Hoffmann, M. R. iCI: Iterative CI toward full CI. _J. Chem. Theory Comput._ 2016, _12_ , 1169–1178, (E)12, 3000 (2016).
  * Zhang et al. 2020 Zhang, N.; Liu, W.; Hoffmann, M. R. Iterative Configuration Interaction with Selection. _J. Chem. Theory Comput._ 2020,
  * Zhang et al. 2021 Zhang, N.; Liu, W.; Hoffmann, M. R. Further Development of iCIPT2 for Strongly Correlated Electrons. _J. Chem. Theory Comput._ 2021, _17_ , 949–964
  * Wu et al. 2011 Wu, F.; Liu, W.; Zhang, Y.; Li, Z. Linear-scaling time-dependent density functional theory based on the idea of “from fragments to molecule”. _J. Chem. Theory Comput._ 2011, _7_ , 3643–3660
  * Li et al. 2017 Li, H.; Liu, W.; Suo, B. Localization of open-shell molecular orbitals via least change from fragments to molecule. _J. Chem. Phys._ 2017, _146_ , 104104
  * Liu et al. 1997 Liu, W.; Hong, G.; Dai, D.; Li, L.; Dolg, M. The Beijing 4-component density functional theory program package (BDF) and its application to EuO, EuS, YbO and YbS. _Theor. Chem. Acc._ 1997, _96_ , 75–83
  * Zhang et al. 2020 Zhang, Y.; Suo, B.; Wang, Z.; Zhang, N.; Li, Z.; Lei, Y.; Zou, W.; Gao, J.; Peng, D.; Pu, Z.; Xiao, Y.; Sun, Q.; Wang, F.; Ma, Y.; Wang, X.; Guo, Y.; Liu, W. BDF: A relativistic electronic structure program package. _J. Chem. Phys._ 2020, _152_ , 064113
  * Wang and Gao 2015 Wang, Y.; Gao, J. Projected Hybrid Orbitals: A General QM/MM Method. _J. Phys. Chem. B_ 2015, _119_ , 1213–1224
  * Boys 1960 Boys, S. F. Construction of some molecular orbitals to be approximately invariant for changes from one molecule to another. _Rev. Mod. Phys._ 1960, _32_ , 296–299
  * Boughton and Pulay 1993 Boughton, J. W.; Pulay, P. Comparison of the Boys and Pipek-Mezey localizations in the local correlation approach and automatic virtual basis selection. _J. Comput. Chem._ 1993, _14_ , 736–740
  * Kutzelnigg and Liu 2005 Kutzelnigg, W.; Liu, W. Quasirelativistic theory equivalent to fully relativistic theory. _J. Chem. Phys._ 2005, _123_ , 241102
  * Liu and Peng 2009 Liu, W.; Peng, D. Exact two-component Hamiltonians revisited. _J. Chem. Phys._ 2009, _131_ , 031104
  * Liu 2010 Liu, W. Ideas of relativistic quantum chemistry. _Mol. Phys._ 2010, _108_ , 1679–1706
  * Liu et al. 2003 Liu, W.; Wang, F.; Li, L. The Beijing density functional (BDF) program package: Methodologies and applications. _J. Theor. Comput. Chem._ 2003, _2_ , 257–272
  * Liu et al. 2004 Liu, W.; Wang, F.; Li, L. In _Recent Advances in Relativistic Molecular Theory_ ; Hirao, K., Ishikawa, Y., Eds.; World Scientific: Singapore, 2004; pp 257–282
  * Liu et al. 2004 Liu, W.; Wang, F.; Li, L. In _Encyclopedia of Computational Chemistry_ ; von Ragué Schleyer, P., Allinger, N. L., Clark, T., Gasteiger, J., Kollman, P. A., Schaefer III, H. F., Eds.; Wiley: Chichester, UK, 2004
  * 117 Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange.
  * Stephens et al. 1994 Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. _J. Phys. Chem._ 1994, _98_ , 11623–11627
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. _Phys. Chem. Chem. Phys._ 2005, _7_ , 3297–3305
  * Amara et al. 2000 Amara, P.; Field, M. J.; Alhambra, C.; Gao, J. L. The generalized hybrid orbital method for combined quantum mechanical/molecular mechanical calculations: formulation and tests of the analytical derivatives. _Theor. Chem. Acc._ 2000, _104_ , 336–343
  * Eckard and Exner 2009 Eckard, S.; Exner, T. E. Improvements in the generalized hybrid orbital method. _Int. J. Quantum Chem._ 2009, _109_ , 1451–1463
  * Jung et al. 2007 Jung, J.; Choi, C. H.; Sugita, Y.; Ten-no, S. New implementation of a combined quantum mechanical and molecular mechanical method using modified generalized hybrid orbitals. _J. Chem. Phys._ 2007, _127_ , 204102