Fragmentation-Based Linear-Scaling Method for Strongly Correlated Systems: Divide-and-Conquer Hartree–Fock–Bogoliubov Method, Its Energy Gradient, and Applications to Graphene Nano-Ribbon Systems

1. Introduction

A simple yet accurate description of the static electron correlation, especially in large systems, is still a challenging problem in the electronic structure theory. The static electron correlation is often closely related to the chemical structure. For example, the bond-alternating polyene chain has a weaker static correlation than the chain with a uniform bond length. The transition-state structures often suffer from stronger static correlation than the stable structures. Polyacenes [1] and graphene nanoribbons (GNRs) [2] exhibit different strengths of static correlation depending on their size, shape, and topology (sheet, moebius, or twisted) [3]. In these systems, the strength of the static correlation is related to the diradical character [4]. The diradical character was originally defined in the context of the valence configuration interaction wavefunction [5]. The definition based on the natural orbital occupation numbers obtained from the unrestricted Hartree–Fock (UHF) density matrix was also proposed [6,7]. In particular, Nakano has conducted detailed research into the relationship between the diradical character and excited state properties such as the optical response [8,9,10].

The standard method for describing the static electron correlation is the active-space self-consistent field (SCF) method, typified by the complete active-space SCF (CASSCF) method [11]. These methods require the specification of the active space that defines the strongly correlated orbitals and electrons. Although the recent advent of the density matrix renormalization group method [12,13,14,15] allows much larger active space to be handled than before, its application is still limited to medium-sized systems.

Another way to account for the static correlation is based on the two-electron wavefunction, called geminal [16,17,18]. These schemes, including generalized valence bond methods [19]. only treat two opposite-spin electrons in the variational way, and the other dynamical correlation may be incorporated in the perturbational manner [20,21]. Although the formal computational scaling for these methods can be lower than the CASSCF method with large active space, the simultaneous optimization of orbitals and CI coefficients usually causes serious convergence problems. Note that the geminal method is also related to the pair coupled cluster theories [22,23].

An alternative approach to static correlation problems inspired by the geminal-based method is the Hartree–Fock–Bogoliubov (HFB) method reformulated for quantum chemical problems [24]. Since the HFB wavefunction is related to the antisymmetrized geminal power wavefunction by the particle-number projection [25], the HFB method can effectively describe the static electron correlation, although the obtained as-is HFB wavefunction is contaminated by the different particle-number wavefunctions. The formal computational scaling for the HFB method is similar to that of the Hartree–Fock (HF) method with doubled basis functions. However, we have proposed a method to apply the linear-scaling divide-and-conquer (DC) scheme [26,27,28] to the HFB method [29].

To predict molecular geometry from the quantum chemical calculations, energy gradient must be formulated and implemented in the computational code [30]. One of the authors (MK) has formulated the HFB energy gradient [31]. Note that Chai and coworkers have proposed thermally-assisted occupation density functional theory (TAO-DFT) [32,33,34]. This method can consider the static electron correlation by introducing fictitious electronic temperature and its energy gradient is available. On the other hand, in the DC method, the energy gradient was proposed in the earliest Yang and Lee paper in HF and DFT frameworks [26]. Since this DC energy gradient expression includes an approximation, MK and coworkers delived the exact DC-HF/DFT energy gradient and proposed an improved approximate energy gradient expression [35].

In this paper, we derived two energy gradient expressions for the DC-HFB method in accordance with Refs. [YangDM-DC,KKASNDC-SCFgrad]. The accuracy of these two gradient expressions is investigated in calculations of GNRs including polyacenes. In addition, the strength of the static electron correlation is discussed from the natural orbitals of the HFB wavefunction and their occupation numbers, which are closely related to the diradical, tetraradical, and hexaradical characters. Although the simultaneous consideration of the dynamical electron correlation as well as the static correlation is necessary for the quantitative discussion, we focus here only on the static correlation for the basic evaluation of the method. The organization of this paper is the following. The theoretical perspective of the present method is given in Section 2. The numerical assessment in calculations of GNRs are presented in Section 3. Finally, Section 4 gives the concluding remarks.

2. Theory

2.1. Restricted DC-HFB Method

We here focus on the restricted HFB method for simplicity [31]. Through this paper, Greek-letter indices $\{\mu ,\nu ,\cdots \}$ refer to the basis functions [i.e., atomic orbitals (AOs) that can be non-orthogonal] and $\{i,j,\cdots \}$ refer to the quasiparticle orbitals with negative eigenvalues. The restricted HFB energy can be expressed as the functional of the density matrix $\overline{\mathbf{D}}\equiv {\mathbf{D}}^{\uparrow \uparrow }$ and the paring matrix $\overline{\mathbf{K}}\equiv {\mathbf{K}}^{\uparrow \downarrow }=-{\mathbf{K}}^{\downarrow \uparrow }$:

$$\begin{array}{c}\hfill E_{\mathrm{HFB}}=2\sum _{\mu \nu }{h}_{\mu \nu }{\overline{D}}_{\nu \mu }+\sum _{\mu \nu \lambda \rho }(2\langle \mu \nu |\lambda \rho \rangle -\langle \mu \nu |\rho \lambda \rangle ){\overline{D}}_{\lambda \mu }{\overline{D}}_{\rho \nu }-\zeta \sum _{\mu \nu \lambda \rho }\langle \mu \nu |\lambda \rho \rangle {\overline{K}}_{\mu \nu }{\overline{K}}_{\lambda \rho }+2\Lambda \left(N-\sum _{\mu \nu }{S}_{\mu \nu }{\overline{D}}_{\nu \mu }\right),\end{array}$$

(1)

where $\langle \mu \nu |\lambda \rho \rangle =\int \int \mu \left({\mathbf{r}}_1\right)\nu \left({\mathbf{r}}_2\right)r_{12}^{-1}\lambda \left({\mathbf{r}}_1\right)\rho \left({\mathbf{r}}_2\right)d{\mathbf{r}}_{1}d{\mathbf{r}}_2$, $\mathbf{h}$ is the one-electron integral matrix, $\mathbf{S}$ represents the overlap matrix of the non-orthogonal basis functions, and $\zeta$ is a scaling factor introduced by Staroverov and Scuseria [24]. $\Lambda$ is the chemical potential that is tuned to keep the number of electrons in the system to be $2N$. The density and pairing matrices constitute the generalized density matrix $\overline{\mathbf{R}}$:

$$\begin{array}{c}\hfill \overline{\mathbf{R}}=\left(\begin{array}{cc}\overline{\mathbf{D}}& \overline{\mathbf{K}}\\ \overline{\mathbf{K}}& {\mathbf{S}}^{-1}-\overline{\mathbf{D}}\end{array}\right).\end{array}$$

(2)

Variation of the energy (1) with respect to the generalized density matrix (2) leads to the restricted HFB Hamiltonian

$$\begin{array}{c}\hfill \mathbf{H}=\left(\begin{array}{cc}\tilde{\mathbf{F}}& \Delta \\ \Delta & -\tilde{\mathbf{F}}\end{array}\right)\end{array}$$

(3)

where $\tilde{\mathbf{F}}$ is the shifted Fock matrix,

$$\begin{array}{ccc}\hfill {\tilde{F}}_{\mu \lambda }& =& h_{\mu \lambda }-\Lambda S_{\mu \lambda }+\sum _{\nu \rho }(2\langle \mu \nu |\lambda \rho \rangle -\langle \mu \nu |\rho \lambda \rangle ){\overline{D}}_{\rho \nu }\hfill \\ & =& {\tilde{h}}_{\mu \lambda }+\sum _{\nu \rho }(2\langle \mu \nu |\lambda \rho \rangle -\langle \mu \nu |\rho \lambda \rangle ){\overline{D}}_{\rho \nu },\hfill \end{array}$$

(4)

with $\tilde{\mathbf{h}}=\mathbf{h}-\Lambda \mathbf{S}$, and $\Delta$ is the pairing-field matrix,

$$\begin{array}{c}\hfill {\Delta }_{\mu \nu }=-\zeta \sum _{\lambda \rho }\langle \mu \nu |\lambda \rho \rangle {\overline{K}}_{\lambda \rho }.\end{array}$$

(5)

The eigenvectors of the HFB Hamiltonian (3), ${\left({\mathbf{Y}}_i^{\mathrm{T}}{\mathbf{X}}_i^{\mathrm{T}}\right)}^{\mathrm{T}}$,

$$\begin{array}{c}\hfill \left(\begin{array}{cc}\tilde{\mathbf{F}}& \Delta \\ \Delta & -\tilde{\mathbf{F}}\end{array}\right)\left(\begin{array}{c}{\mathbf{Y}}_i\\ {\mathbf{X}}_i\end{array}\right)=-\left(\begin{array}{cc}\mathbf{S}& \mathbf{0}\\ \mathbf{0}& \mathbf{S}\end{array}\right)\left(\begin{array}{c}{\mathbf{Y}}_i\\ {\mathbf{X}}_i\end{array}\right){\epsilon }_i,\end{array}$$

(6)

represent the quasiparticle orbitals with negative eigenvalues (${\epsilon }_i\ge 0$), which are connected to the density and pairing matrices as follows:

$$\begin{array}{c}\hfill {\overline{D}}_{\mu \nu }=\sum _i^M{Y}_{\mu i}{Y}_{\nu i},\end{array}$$

(7)

$$\begin{array}{c}\hfill {\overline{K}}_{\mu \nu }=\sum _i^M{Y}_{\mu i}{X}_{\nu i}=\sum _i^M{X}_{\mu i}{Y}_{\nu i}.\end{array}$$

(8)

Here, M represents the number of quasiparticle orbitals. Note that the generalized density matrix $\overline{\mathbf{R}}$ of the converged HFB solution is idempotent, namely,

$$\begin{array}{c}\hfill \overline{\mathbf{D}}\mathbf{S}\overline{\mathbf{D}}-\overline{\mathbf{D}}=-\overline{\mathbf{K}}\mathbf{S}\overline{\mathbf{K}},\end{array}$$

(9)

$$\begin{array}{c}\hfill \overline{\mathbf{D}}\mathbf{S}\overline{\mathbf{K}}=\overline{\mathbf{K}}\mathbf{S}\overline{\mathbf{D}}.\end{array}$$

(10)

In the DC method, the entire system is spatially divided into disjoint subsystems called the central region, where the set of AOs in subsystem $\alpha$ is denoted by $\mathcal{S}\left(\alpha \right)$ and the union of $\mathcal{S}\left(\alpha \right)$ becomes the set of basis functions over the entire system denoted by $\mathcal{T}$ as follows:

$$\begin{array}{c}\hfill \mathcal{S}\left(\alpha \right)\cap \mathcal{S}\left(\beta \right)=\varnothing \forall \alpha \ne \beta ,\end{array}$$

(11)

$$\begin{array}{c}\hfill \bigcup _{\alpha }\mathcal{S}\left(\alpha \right)=\mathcal{T}.\end{array}$$

(12)

When constructing MOs for subsystem $\alpha$, the AOs of atoms adjacent to the central region, called the buffer region denoted by $\mathcal{B}\left(\alpha \right)$, are taken into consideration in addition to $\mathcal{S}\left(\alpha \right)$. The union of the central and buffer regions of subsystem $\alpha$ is called the localization region, in which the set of AOs is denoted by $\mathcal{L}\left(\alpha \right)$:

$$\begin{array}{c}\hfill \mathcal{S}\left(\alpha \right)\cup \mathcal{B}\left(\alpha \right)\equiv \mathcal{L}\left(\alpha \right).\end{array}$$

(13)

In the DC-HFB method [29], the generalized density matrix of Equation (2) is approximated as sum of subsystem contributions:

$$\begin{array}{c}\hfill {\overline{D}}_{\mu \nu }\approx {\overline{D}}_{\mu \nu }^{\mathrm{DC}}=\sum _{\alpha }{P}_{\mu \nu }^{\alpha }{\overline{D}}_{\mu \nu }^{\alpha },\end{array}$$

(14)

$$\begin{array}{c}\hfill {\overline{K}}_{\mu \nu }\approx {\overline{K}}_{\mu \nu }^{\mathrm{DC}}=\sum _{\alpha }{P}_{\mu \nu }^{\alpha }{\overline{K}}_{\mu \nu }^{\alpha },\end{array}$$

(15)

where ${\mathbf{P}}^{\alpha }$ is the partition matrix defined by

$$\begin{array}{c}\hfill P_{\mu \nu }^{\alpha }=\left\{\begin{array}{cc}1\hfill & (\mu \in \mathcal{S}(\alpha )\wedge \nu \in \mathcal{S}(\alpha \left)\right)\hfill \\ 1/2\hfill & (\mu \in \mathcal{S}(\alpha )\wedge \nu \in \mathcal{B}(\alpha ),\mathrm{or}\mathit{vice}\mathit{versa})\hfill \\ 0\hfill & \left(\mathrm{otherwise}\right).\hfill \end{array}\right.\end{array}$$

(16)

Using the submatrices of $\tilde{\mathbf{F}}$ and $\Delta$ for $\mathcal{L}\left(\alpha \right)$, denoted by ${\tilde{\mathbf{F}}}^{\alpha }$ and ${\Delta }^{\alpha }$, subsystem density and pairing matrices are obtained with the subsystem quasiparticle orbitals obtained as the eigenvectors of the following subsystem HFB equation

$$\begin{array}{c}\hfill \left(\begin{array}{cc}{\tilde{\mathbf{F}}}^{\alpha }& {\Delta }^{\alpha }\\ {\Delta }^{\alpha }& -{\tilde{\mathbf{F}}}^{\alpha }\end{array}\right)\left(\begin{array}{c}{\mathbf{Y}}_i^{\alpha }\\ {\mathbf{X}}_i^{\alpha }\end{array}\right)=-\left(\begin{array}{cc}{\mathbf{S}}^{\alpha }& \mathbf{0}\\ \mathbf{0}& {\mathbf{S}}^{\alpha }\end{array}\right)\left(\begin{array}{c}{\mathbf{Y}}_i^{\alpha }\\ {\mathbf{X}}_i^{\alpha }\end{array}\right){\epsilon }_i^{\alpha },\end{array}$$

(17)

as

$$\begin{array}{c}\hfill {\overline{D}}_{\mu \nu }^{\alpha }=\sum _i^{M^{\alpha }}{Y}_{\mu i}^{\alpha }{Y}_{\nu i}^{\alpha },\end{array}$$

(18)

$$\begin{array}{c}\hfill {\overline{K}}_{\mu \nu }^{\alpha }=\sum _i^{M^{\alpha }}{Y}_{\mu i}^{\alpha }{X}_{\nu i}^{\alpha }=\sum _i^{M^{\alpha }}{X}_{\mu i}^{\alpha }{Y}_{\nu i}^{\alpha },\end{array}$$

(19)

with the number of quasiparticle orbitals in subysystem $\alpha$, $M^{\alpha }$. Another expression that is consistent to the DC-HF method with the finite electronic temperature of $1/\left(k_{\mathrm{B}}\beta \right)$ is the following:

$$\begin{array}{c}\hfill {\overline{D}}_{\mu \nu }^{\alpha }=\sum _p{f}_{\beta }\left({\epsilon }_p^{\alpha }\right)Y_{\mu p}^{\alpha }{Y}_{\nu p}^{\alpha },\end{array}$$

(20)

$$\begin{array}{c}\hfill {\overline{K}}_{\mu \nu }^{\alpha }=\sum _p{f}_{\beta }\left({\epsilon }_p^{\alpha }\right)Y_{\mu p}^{\alpha }{X}_{\nu p}^{\alpha }=\sum _p{f}_{\beta }\left({\epsilon }_p^{\alpha }\right)X_{\mu p}^{\alpha }{Y}_{\nu p}^{\alpha },\end{array}$$

(21)

where $f_{\beta }\left(x\right)={\left[1+\exp (-\beta x)\right]}^{-1}$ is the Fermi distribution function. Note that the summation over p runs all the quasiparticle orbitals regardless of the sign of the eigenvalues.

2.2. DC-HFB Energy Gradient with Respect to Atomic Coordinate

According to Ref. [KHFBgrad], the differentiation of Equation (1) with respect to atomic coordinate Q gives

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial E_{\mathrm{HFB}}}{\partial Q}}& =& 2\sum _{\mu \nu }{\displaystyle \frac{\partial (h_{\mu \nu }-\Lambda S_{\mu \nu })}{\partial Q}}{\overline{D}}_{\nu \mu }+\sum _{\mu \nu \lambda \rho }{\displaystyle \frac{\partial (2\langle \mu \nu |\lambda \rho \rangle -\langle \mu \nu |\rho \lambda \rangle )}{\partial Q}}{\overline{D}}_{\lambda \mu }{\overline{D}}_{\rho \nu }-\zeta \sum _{\mu \nu \lambda \rho }{\displaystyle \frac{\partial \langle \mu \nu |\lambda \rho \rangle }{\partial Q}}{\overline{K}}_{\mu \nu }{\overline{K}}_{\lambda \rho }\hfill \\ & & +2\sum _{\mu \nu }(h_{\mu \nu }-\Lambda S_{\mu \nu }){\displaystyle \frac{\partial {\overline{D}}_{\nu \mu }}{\partial Q}}+2\sum _{\mu \nu \lambda \rho }(2\langle \mu \nu |\lambda \rho \rangle -\langle \mu \nu |\rho \lambda \rangle ){\displaystyle \frac{\partial {\overline{D}}_{\lambda \mu }}{\partial Q}}{\overline{D}}_{\rho \nu }-2\zeta \sum _{\mu \nu \lambda \rho }\langle \mu \nu |\lambda \rho \rangle {\displaystyle \frac{\partial {\overline{K}}_{\mu \nu }}{\partial Q}}{\overline{K}}_{\lambda \rho }\hfill \\ & & +2{\displaystyle \frac{\partial \Lambda }{\partial Q}}\left(N-\sum _{\mu \nu }{S}_{\mu \nu }{\overline{D}}_{\nu \mu }\right)\hfill \\ & =& G_{\mathrm{H}--\mathrm{F}}^Q+2\mathrm{tr}\left(\tilde{\mathbf{F}}{\displaystyle \frac{\partial \overline{\mathbf{D}}}{\partial Q}}+\Delta {\displaystyle \frac{\partial \overline{\mathbf{K}}}{\partial Q}}\right).\hfill \end{array}$$

(22)

The first term, $G_{\mathrm{H}--\mathrm{F}}^Q$, represents the Hellmann–Feynman (H–F) force that relates to the derivatives of the Hamiltonian matrix elements, and the second term is the so-called Pulay force. In case of the variational solution, the evaluation of $\overline{\mathbf{D}}$ and $\overline{\mathbf{K}}$ derivatives with respect to Q can be avoided by using the idempotency condition of the generalized density matrix, $\overline{\mathbf{R}}$, as [31]

$$\begin{array}{c}\hfill 2\mathrm{tr}\left(\tilde{\mathbf{F}}{\displaystyle \frac{\partial \overline{\mathbf{D}}}{\partial Q}}+\Delta {\displaystyle \frac{\partial \overline{\mathbf{K}}}{\partial Q}}\right)=-2\mathrm{tr}\left(\overline{\mathbf{W}}{\displaystyle \frac{\partial \mathbf{S}}{\partial Q}}\right),\end{array}$$

(23)

where $\overline{\mathbf{W}}$ is the energy-weighted density matrix defined by

$$\begin{array}{c}\hfill {\overline{W}}_{\mu \nu }=-\sum _i^M{\epsilon }_i{Y}_{\mu i}{Y}_{\nu i}.\end{array}$$

(24)

Note that the eigenvalues of Equation (6) are $\{-{\epsilon }_i\}$.

As similar to the DC-HF energy gradient [35], the Pulay force of the DC-HFB energy gradient cannot be simplified to the form of Equation (23) because the DC-HFB density and pairing matrices are not variationally determined. However, as Yang and Lee [26] proposed, the Pulay force can be approximated in the similar manner to the DC density and pairing matrices of Equations (14) and (15) as

$$\begin{array}{c}\hfill {\overline{W}}_{\mu \nu }\approx {\overline{W}}_{\mu \nu }^{\mathrm{DC}}=\sum _{\alpha }{P}_{\mu \nu }^{\alpha }{\overline{W}}_{\mu \nu }^{\alpha },\end{array}$$

(25)

$$\begin{array}{c}\hfill {\overline{W}}_{\mu \nu }^{\alpha }=-\sum _i^{M^{\alpha }}{\epsilon }_i^{\alpha }{Y}_{\mu i}^{\alpha }{Y}_{\nu i}^{\alpha },\end{array}$$

(26)

Therefore, the DC-HFB energy gradient can be approximated as the following:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial E_{\mathrm{DC}-\mathrm{HFB}}}{\partial Q}}\approx G_{\mathrm{H}-\mathrm{F}}^Q-2\mathrm{tr}\left({\overline{\mathbf{W}}}^{\mathrm{DC}}{\displaystyle \frac{\partial \mathbf{S}}{\partial Q}}\right).\end{array}$$

(27)

According to Ref. [KKASNDC-SCFgrad], this formula is referred to as the YL formula throughout this paper.

Alternatively, the Pulay force in the DC-HFB energy can be reformulated as the following:

$$\begin{array}{c}\hfill 2\mathrm{tr}\left(\tilde{\mathbf{F}}{\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\mathrm{DC}}}{\partial Q}}+\Delta {\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\mathrm{DC}}}{\partial Q}}\right)=2\sum _{\alpha }\sum _{\mu \in \mathcal{S}\left(\alpha \right)}{\left({\tilde{\mathbf{F}}}^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}}+{\Delta }^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\alpha }}{\partial Q}}\right)}_{\mu \mu }.\end{array}$$

(28)

On the analogy to the formulation in Ref. [KKASNDC-SCFgrad], we here assume that the subsystem generalized density matrix,

$$\begin{array}{c}\hfill {\overline{\mathbf{R}}}^{\alpha }=\left(\begin{array}{cc}{\overline{\mathbf{D}}}^{\alpha }& {\overline{\mathbf{K}}}^{\alpha }\\ {\overline{\mathbf{K}}}^{\alpha }& {\left({\mathbf{S}}^{\alpha }\right)}^{-1}-{\overline{\mathbf{D}}}^{\alpha }\end{array}\right),\end{array}$$

(29)

also satisfies the idempotency condition, leading

$$\begin{array}{c}\hfill {\overline{\mathbf{D}}}^{\alpha }{\mathbf{S}}^{\alpha }{\overline{\mathbf{D}}}^{\alpha }-{\overline{\mathbf{D}}}^{\alpha }=-{\overline{\mathbf{K}}}^{\alpha }{\mathbf{S}}^{\alpha }{\overline{\mathbf{K}}}^{\alpha },\end{array}$$

(30)

$$\begin{array}{c}\hfill {\overline{\mathbf{D}}}^{\alpha }{\mathbf{S}}^{\alpha }{\overline{\mathbf{K}}}^{\alpha }={\overline{\mathbf{K}}}^{\alpha }{\mathbf{S}}^{\alpha }{\overline{\mathbf{D}}}^{\alpha }.\end{array}$$

(31)

Furthermore, we regard that the subsystem generalized density matrices, $\left\{{\overline{\mathbf{R}}}^{\alpha }\right\}$, are variationally determined. Then, any simultaneous solution of the derivatives of Equations (30) and (31),

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}}{\mathbf{S}}^{\alpha }{\overline{\mathbf{D}}}^{\alpha }+{\overline{\mathbf{D}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{D}}}^{\alpha }+{\overline{\mathbf{D}}}^{\alpha }{\mathbf{S}}^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}}-{\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}}=-{\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\alpha }}{\partial Q}}{\mathbf{S}}^{\alpha }{\overline{\mathbf{K}}}^{\alpha }-{\overline{\mathbf{K}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{K}}}^{\alpha }-{\overline{\mathbf{K}}}^{\alpha }{\mathbf{S}}^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\alpha }}{\partial Q}},\end{array}$$

(32)

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}}{\mathbf{S}}^{\alpha }{\overline{\mathbf{K}}}^{\alpha }+{\overline{\mathbf{D}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{K}}}^{\alpha }+{\overline{\mathbf{D}}}^{\alpha }{\mathbf{S}}^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\alpha }}{\partial Q}}={\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\alpha }}{\partial Q}}{\mathbf{S}}^{\alpha }{\overline{\mathbf{D}}}^{\alpha }+{\overline{\mathbf{K}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{D}}}^{\alpha }+{\overline{\mathbf{K}}}^{\alpha }{\mathbf{S}}^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}},\end{array}$$

(33)

gives the same force [36]. A simple set of solutions is the following [31]:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}}={\overline{\mathbf{K}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{K}}}^{\alpha }-{\overline{\mathbf{D}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{D}}}^{\alpha }.\end{array}$$

(34)

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\alpha }}{\partial Q}}=-{\overline{\mathbf{K}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{D}}}^{\alpha }-{\overline{\mathbf{D}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{K}}}^{\alpha }.\end{array}$$

(35)

Applying Equations (34) and (35) to the Pulay force of Equation (28) and considering the subsystem HFB equation of Equation (17) leads to

$$\begin{array}{c}\hfill 2\sum _{\alpha }\sum _{\mu \in \mathcal{S}\left(\alpha \right)}{\left({\tilde{\mathbf{F}}}^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{D}}}^{\alpha }}{\partial Q}}+{\Delta }^{\alpha }{\displaystyle \frac{\partial {\overline{\mathbf{K}}}^{\alpha }}{\partial Q}}\right)}_{\mu \mu }=-2\sum _{\alpha }\sum _{\mu \in \mathcal{S}\left(\alpha \right)}{\left({\mathbf{S}}^{\alpha }{\overline{\mathbf{W}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{D}}}^{\alpha }+{\mathbf{S}}^{\alpha }{\overline{\mathbf{V}}}^{\alpha }{\displaystyle \frac{\partial {\mathbf{S}}^{\alpha }}{\partial Q}}{\overline{\mathbf{K}}}^{\alpha }\right)}_{\mu \mu },\end{array}$$

(36)

where ${\overline{\mathbf{V}}}^{\alpha }$ is the subsystem energy-weighted pairing matrix defined by

$$\begin{array}{c}\hfill {\overline{V}}_{\mu \nu }^{\alpha }=-\sum _i^{M^{\alpha }}{\epsilon }_i^{\alpha }{X}_{\mu i}^{\alpha }{Y}_{\nu i}^{\alpha }.\end{array}$$

(37)

Finally, the following approximated DC-HFB energy gradient can be obtained:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial E_{\mathrm{DC}-\mathrm{HFB}}}{\partial Q}}\approx G_{\mathrm{H}--\mathrm{F}}^Q-2\mathrm{tr}\left[{\displaystyle \frac{\partial \mathbf{S}}{\partial Q}}\sum _{\alpha }\left({\mathbf{A}}^{\alpha }+{\mathbf{B}}^{\alpha }\right)\right],\end{array}$$

(38)

where

$$\begin{array}{c}\hfill {\mathbf{A}}^{\alpha }={\overline{\mathbf{D}}}^{\alpha \left[\mathcal{L}\right(\alpha )\times \mathcal{S}(\alpha \left)\right]}{\mathbf{S}}^{\alpha \left[\mathcal{S}\right(\alpha )\times \mathcal{L}(\alpha \left)\right]}{\overline{\mathbf{W}}}^{\alpha },\end{array}$$

(39)

$$\begin{array}{c}\hfill {\mathbf{B}}^{\alpha }={\overline{\mathbf{K}}}^{\alpha \left[\mathcal{L}\right(\alpha )\times \mathcal{S}(\alpha \left)\right]}{\mathbf{S}}^{\alpha \left[\mathcal{S}\right(\alpha )\times \mathcal{L}(\alpha \left)\right]}{\overline{\mathbf{V}}}^{\alpha }.\end{array}$$

(40)

Here, the superscript $\left[\mathcal{L}\right(\alpha )\times \mathcal{S}(\alpha \left)\right]$ represents the submatrix whose row and column indices are $\mathcal{L}\left(\alpha \right)$ and $\mathcal{S}\left(\alpha \right)$, respectively. This formula is referred to as the KKASN formula [35] throughout this paper. As similar to the DC-HF energy gradient [35], ${\mathbf{A}}^{\alpha }$ and ${\mathbf{B}}^{\alpha }$ matrices in Equation (38) can be Hermitized.

3. Numerical Application to Graphene Nanoribbons

Single-layer graphene sheet shows characteristic electronic properties such as fractional quantum Hall effect [37] and ultrafast carrier mobility [38] due to its special band structure derived from the electronic structure confined in a plane. Graphene-based spintronics devices are also expected to be developed, [38,39] although the band gap has to be generated for realizing these devices. GNRs were theoretically expected to generate band gap [40,41] and have also been experimentally reported. Since the electronic properties of GNRs depend on the geometrical factors such as width, length, and edge shape, it is essential to understand their relation with the electronic structure in order to design GNR-based nanodevices. To obtain the correct electronic structures of GNRs, known as strongly correlated systems, high-level ab initio multireference methods such as the density matrix renormalization group (DMRG) [4,42] and the second-order reduced density matrix-based variational method [43] are required. However, due to their large computational demands, their applications to larger GNRs are difficult. In addition, at the current stage, the geometry optimization based on these high-level ab initio methods is not tructable.

Here, we performed the geometry optimizations of GNRs with the DC-HFB energy gradient and investigated the polyradical character of GNRs from the HFB natural orbital occupation. After a (DC-)HFB calculation, the natural orbitals can be obtained by diagonalizing the density matrix, $\overline{\mathbf{D}}$. The eigenvalues correspond to the natural orbital occupation numbers (NOONs), which has been discussed as the indicator of diradical character in UHF calculations [44]. For single-reference systems, NOON is either 0 (unoccupied) or 2 (two-electron occupied), while for singlet diradical systems with $N_{\mathrm{e}}$ electrons, the highest occupied (i.e., $N_{\mathrm{e}}⁄2$-th) and lowest unoccupied (i.e., $(N_{\mathrm{e}}⁄2+1)$-th) natural orbitals (HONO and LUNO) are degenerate, resulting in NOON being 1 [45]. Throughout this paper, the modified version of the Gamess program package [46,47] was used for computations of the standard and DC-HFB energy and gradient.

3.1. Determination of Parameter $\zeta$

Before applying the HFB method to GNRs, the scaling factor, $\zeta$, has to be fixed. Here, this parameter was determined by comparing the optimized geometries of phenalenyl cation, anthracene, and phenanthrene (Figure 1) with those by reference all-$\pi$ CASSCF calculations, which were performed with the OpenMolcas 24.10 program [48]. No dynamical correlation correction was introduced to the CASSCF reference since the HFB method basically cannot account for the dynamical correlation. Throughout this paper, the 6-31G** basis set [49,50] was adopted. The obtained CASSCF geometries for phenalenyl cation, anthracene, and phenanthrene have the $D_{3h}$, $D_{2h}$, and $C_{2v}$ symmetries, respectively, the stabilities of which are confirmed by the frequency analysis.

Table 1 and Table 2 summarize the optimized C–C and C–H bond lengths, respectively, for these molecules obtained with CASSCF and HFB methods. For HFB calculations, the parameter $\zeta$ is varied between 0.75 and 1.0. The mean absolute deviations (MADs) of the HFB results from the CASSCF reference results are also listed. The deviations of the C–C bond lengths are one-order (or more) larger than those of the C–H lengths. The HFB results with $\zeta =1.0$ include obviously larger deviation than the others. For both the C–C and C–H lengths, the HFB results with $\zeta =0.85$ gives the smallest MADs from the CASSCF ones. Hereafter, $\zeta =0.85$ was adopted.

3.2. Polyradicality of Graphene Nanoribbons

The HFB calculations were performed for GNRs of different widths (W) and lengths (L) to investigate the effect of edge structure on polyradicality. The width W and length L are defined as shown in Figure 2. For this sample molecule, we denote $\{L,W\}=\{4,2\}$. NOONs for (HONO$-2$) to (LUNO$+2$) were examined as a measure of polyradicality. The geometries were also optimized at the HFB/6-31G** level of theory.

First, we fixed the width $W=2$ and varied the zigzag edge length $L=2-8$. NOONs for six frontier natural orbitals are summarized in Figure 3. The occupation number difference between HONO and LUNO rapidly decreases as L increases, indicating that the diradical character of the system is increasing along L. In addition, the differences between (HONO$-1$)/(HONO$-2$) and (LUNO$+1$)/(LUNO$+2$) also decreases gradually, indicating the increase of tetra- and hexa-radical characters in these systems.

Next, we fixed the length $L=2$ and varied the armchair edge width $W=1-7$, of which the results are depicted in Figure 4. Again, the occupation number difference between HONO and LUNO decreases as W increases, but the rate of decrease is slower than the $W=2$ result in Figure 3. The reason can be deduced from the shapes of HONO and LUNO. The HONOs and LUNOs of the GNRs with $\{L,W\}=\{8,2\}$ and $\{L,W\}=\{3,4\}$ are illustrated in Figure 5 and Figure 6, respectively. Since the occupation numbers for these NOs are $\sim 1$, the electron distribution represented by these NOs corresponds to the distribution of the singlet diradical. It can be confirmed that the radical electrons are mainly localized on zigzag edges. Therefore, diradical character is enhanced when the zigzag edge is elongated than when the armchair edge is elongated.

Figure 7 (a) shows the HONO-LUNO occupation number difference of GNRs plotted against the number of carbon atoms. For the same number of carbon atoms, the $W=1$ series, in which the zigzag edge increases with the smallest armchair edge fraction, has the highest diradical character. Conversely, the $L=2$ series, in which the armchair edge increases with the smallest zigzag edge fraction, has the lowest diradical character. The difference in the occupation numbers between HONO$-1$ and LUNO$+1$, which is an indicator of tetraradicality, is shown in Figure 7 (b), and a similar trend is seen especially in W-constant series. In the L-constant series, if the number of carbon atoms is taken as the horizontal axis, the $L=2$ and $L=3$ series are reversed in the middle. This is probably due to the fact that radicals tend to localize at the armchair edge compared to the center of the GNR.

3.3. DC-HFB Optimization of GNR Structures

In this subsection, the geometry optimization of the GNR with $\{L,W\}=\{15,1\}$ was performed by the standard and DC-HFB method. In the DC-HFB calculations, two different gradient formalisms, namely KKASN and YL formulas, as well as two different fragmentation schemes, namely type A and B as shown in Figure 8 were investigated.

First the fragmentation scheme was set to type A. Figure 9 shows the optimized C–C bond lengths at the zigzag edge with respect to the sequential position. In the DC calculations, up to $n_{\mathrm{b}}$-th nearest fragments were included in the buffer region. Although alternating single and double bond structure can be found around the end points, the bond alternation rapidly decreases around the center. For DC with $n_{\mathrm{b}}=3$, the structure obtained with the YL formula significantly differs from the standard one. The difference significantly diminishes when adopting KKASN formula. The mean absolute deviations (MADs) from the standard HFB result are $0.0246$ Å and $0.0043$ Å for YL and KKASN formulas, respectively. As the buffer size increases to $n_{\mathrm{b}}=4$, the structure deviation decreases by about one order of magnitude both for KKASN and YL formulas. The MADs are $0.0040$ Å and $0.0004$ Å, respectively.

Figure 10 shows the energy convergence process in the geometry optimization of GNR with $\{L,W\}=\{15,1\}$. The initial structure was set with the uniform bond lengths and angles, namely, $1.4$ Å, $1.1$ Å, and 120 for C–C bond lengths, C–H bond lengths, and bond angles, respectively. For the DC-HFB results with $n_{\mathrm{b}}=4$, the energy convergence process is almost equivalent to the standard HFB one. The deviations of the energy at the optimized structure are $-1.12$ $\mathrm{m}$ and $-0.07$ $\mathrm{m}$ for KKASN and YL formulas, respectively. At first glance, the YL formula seems to produce better result than the KKASN one, but the energy difference between the DC and standard HFB results for the initial structure is $-1.14$ $\mathrm{m}$, so this is due to the error cancellation caused by the slightly poorer structure by the YL formula. For the combination of smaller buffer size of $n_{\mathrm{b}}=3$ and the YL formula, the energy gradually increased from the third iteration and converged to higher energy structure after more than 20 iterations. However, when combined with the KKASN formula, the energy converged as smoothly as the $n_{\mathrm{b}}=4$ case, while the energy is $28.77$ $\mathrm{m}$ lower than the standard HFB result.

Finally, the dependence on the fragmentation scheme was investigated. Figure 11 shows the optimized C–C bond lengths at the zigzag edge with respect to the sequential position. The type A results are equivalent to those given in Figure 9. Especially for the small buffer size of $n_{\mathrm{b}}=3$, the trend is different between different fragmentation schemes, namely, the type A fragmentation emphasizes the bond alternation than the type B one. However, for the large buffer size of $n_{\mathrm{b}}=4$, the geometry is almost identical. For the type B fragmentation, the MADs from the standard HFB result are $0.0026$ Å and $0.0001$ Å for $n_{\mathrm{b}}=3$ and 4, respectively, which are slightly smaller than for the type A fragmentation.

4. Concluding Remarks

This study presented the divide-and-conquer Hartree–-Fock-–Bogoliubov (DC-HFB) method as an effective approach for the treatment of static electron correlation in large molecular systems. The analysis of natural orbitals obtained from the HFB denstiy matrix revealed that GNRs exhibit strong polyradical character, particularly in structures with extended zigzag edges, which is consistent with previous computational reports. The proposed two approximate energy gradient formulations, specifically the YL and KKASN methods, demonstrated their utility in the geometry optimization of graphene nanoribbons (GNRs). In general, the KKASN method gives more accurate gradient and geometry than the YL one. These results highlight the DC-HFB method’s scalability and accuracy for large-scale systems, making it a promising tool for quantum chemical calculations in strongly correlated systems. As mentioned in the introduction, incorporating dynamical electron correlation is required to further enhance the quantitative accuracy of the method. A simple scheme to incorporate dynamical electron correlation is to combine the HFB method with the Kohn-Sham density functional theory (DFT), which has already been accomplished by Tsuchimochi et al. [51]. Its DC-based linear-scaling implementation as well as its energy gradient formulation will be straightforward, since the DC method was originally proposed in DFT.

Recently, we have proposed a time-dependent (TD) HFB method for the excited-state calculation of molecular systems [52]. Since the real-time propagation scheme is adopted in this excited-state calculation method, the DC-TD-HFB formulation will be possible, opening up the excited-state calculation of large-scale statically correlated systems. In addition, the HFB wavefunction (and of course the DC-HFB method) suffers from the contamination of the different particle number wavefunctions, which is related to the spin contamination in the UHF wavefunction. The projected HFB or UHF method, which extracts the desired particle-number or spin states, respectively, from the contaminated wavefunction, has been proposed [25,53]. The application of the DC scheme to these projected HF theories is currently underway and will be presented in the near future.