- Email: [email protected]
Fractional-occupation-number based divide-and-conquer coupled-cluster theory
Accepted Manuscript Research paper Fractional-occupation-number based divide-and-conquer coupled-cluster theory Takeshi Yoshikawa, Hiromi Nakai PII: DOI: Reference:
S0009-2614(18)30787-5 https://doi.org/10.1016/j.cplett.2018.09.056 CPLETT 35968
To appear in:
Chemical Physics Letters
Received Date: Accepted Date:
4 August 2018 24 September 2018
Please cite this article as: T. Yoshikawa, H. Nakai, Fractional-occupation-number based divide-and-conquer coupled-cluster theory, Chemical Physics Letters (2018), doi: https://doi.org/10.1016/j.cplett.2018.09.056
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Fractional-occupation-number based divide-and-conquer coupled-cluster theory
Takeshi Yoshikawaa, Hiromi Nakaia-c *
a
Department of Chemistry and Biochemistry, School of Advanced Science and
Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555, Japan b
Research Institute for Science and Engineering, Waseda University, 3-4-1 Okubo,
Shinjuku-ku, Tokyo 169-8555, Japan c
Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University,
Katsura, Kyoto 615-8520, Japan _____________________________________________________________________ *
Corresponding author. FAX: +81-3-3205-2504, E-mail address: [email protected] (H. Nakai). URL: http://www.chem.waseda.ac.jp/nakai/ (H. Nakai).
1
Abstract We have extended the divide-and-conquer (DC) coupled-cluster with singles and doubles (CCSD) to a fractional occupation number (FON) formalism, denoted as FON-DC-CCSD, using the thermal Wick theorem. The motivation is to address the inconsistency in the treatment of orbital occupations between the DC-based Hartree– Fock and the DC-CCSD methods, which adopt the Fermi distribution function and the step function for orbital occupation, respectively. Numerical applications involving polyene chains and single-walled carbon nanotubes confirm that the proposed FON-DC-CCSD method reduces both energy errors and computational costs compared with the conventional DC-CCSD method.
2
1. Introduction The coupled-cluster (CC) theory is a robust theoretical approach that systematically improves the accuracy of calculation of the total energy and the wavefunction of atoms and molecules by considering the increase in the order of electron excitation. Based on the correlated method, the energetics of various molecules can be predicted with chemical accuracy. The CC method with single, double, and perturbative triple excitations (CCSD(T)) has been widely used in practical applications and benchmark studies over the years [1,2]. However, its high scaling with respect to the system size N remains a severe impediment for large systems. For example, the computational costs of CCSD and CCSD(T) are asymptotically O(N6) and O(N7), respectively. Many researchers have reported a wide variety of linear-scaling methods based on fragmentation approaches [3-19], which partition the system under consideration into several subsystems. Li et al. [9] categorized the fragmentation approaches into energy-based approaches [3-14] and density matrix (DM)-based approaches [15-29]. Energy-based approaches directly evaluate the total energy of the system from subsystem contributions. The approach is generally applicable to various theory levels such as Hartree–Fock (HF), Møller–Plesset perturbation, and CC theories. In the
3
DM-based approaches, the total energy is calculated from the assembled DM, which is constructed from the DMs of a series of subsystems. One of the features of the DM-based approaches is that a non-integer number of electrons can be treated for each subsystem. Thus, it is possible to reasonably treat charge- and/or spin-delocalized systems using the DM-based approach. Gordon et al. [30] recognized one-step and two-step techniques for the energy-based and DM-based approaches, respectively. The divide-and-conquer (DC) method proposed by Yang and Lee [15,16] can be directly applied to the density functional theory, and is an example of a DM-based approach. These authors investigated the performance of this approach for calculations including HF exchange interactions [20-22]. Furthermore, DC-based correlation methods [23-29] were developed by exploiting energy density analysis (EDA) techniques [31]. The correlation energy of the total system is evaluated by summing up subsystem contributions. Therefore, the DC-based correlation methods are categorized as an energy-based approach. Note that the subsystem orbitals are obtained by diagonalizing the subsystem Fock matrices, which are constructed using the assembled DM. However, the occupation numbers of the subsystem orbitals need to be approximated from non-integer to integer, to solve the subsystem correlation problem. Recently, Kobayashi and Taketsugu proposed the DC-based second-order Møller–
4
Plesset perturbation (DC-MP2) method based on the thermal Wick theorem [32-34] and improved its accuracy [35]. In this study, we have extended the fractional occupation number (FON) technique to the DC-CCSD theory, denoted as FON-DC-CCSD. The rest of this report is organized as follows: in Section 2, we present the theoretical aspects; Section 3 describes the numerical applications; Section 4 presents the conclusion.
5
2. Theory 2.1. DC-HF method This subsection provides a brief explanation of the DC-HF method, the results of which are the basis for the DC-CCSD calculations described in the following subsection. For simplicity, we consider the case of closed-shell systems. The extension to open-shell systems is straightforward. In the DC-HF method, the system under consideration is spatially divided into disjoint subsystems, which are called the central regions. A set of atomic orbitals (AOs) corresponding to a central region α is denoted by S ( . To improve the description of the subsystem, the area adjacent to the central region, which is called the buffer region, is considered during the expansion of the subsystem’s molecular orbitals (MOs). A set of AOs corresponding to the buffer region of the subsystem α, which is denoted by B( , is added to S ( to construct a set of AOs in the localization region of the
subsystem α, L( , as follows: L( S( B .
(1)
The DC-HF method evaluates the MOs C and the orbital energies of the subsystem α by solving the following local Roothaan–Hall equation for the localization region:
6
F C S C ,
(2)
where Sα and Fα represent the local overlap and Fock matrices, respectively, which are extracted from the total matrices. The DM of the entire system, D, is constructed from the local DM for the subsystem α, Dα, as follows: . D P D
(3)
Pα represents the partition matrix with elements given by
1 P 1 2 0
S S B S S B
(4)
otherwise.
The local DM for the subsystem α, Dα, is obtained using the subsystem MOs Cα, the subsystem orbital energies ɛα, and the common Fermi level F D 2 f F p Cp*Cp ,
(5)
f x 1 exp x .
(6)
p
1
Here, indices {μ, ν} and {p, q} refer to AOs and MOs, respectively. Unlike the standard HF method, the DM is constructed on the basis of the Fermi distribution function fβ(x) for the inverse temperature kBT
1
instead of the step function θ(x). Note that the
Fermi distribution function is identical to the step function in the limit of β → ∞, i.e., T → 0 K. As the occupation numbers of the subsystem orbitals are non-integer as
7
shown in Fig. 1(a), it is possible to reasonably treat charge- and/or spin-delocalized systems. The common Fermi level F is determined by a constraint on the total number of electrons in the entire system, N e , as follows:
Ne D S .
(7)
,
Similarly, the number of electrons for the subsystem α, N eL , is determined by
NeL =
D
, L
S ,
(8)
whose value is a non-integer because of fractional occupation. Note that N e is not consistent with the sum N eL because of double counting of the number of electrons in the buffer regions. The partition matrix Pα in Eq. (4) avoids such double counting. From Eqs. (3) and (7), Ne P D S N eS .
,
(9)
Here, N eS is defined and rewritten as follows:
NeS =
S
Thus, N e
P
, L
D S
D
S L
S .
(10)
is regarded as the number of electrons in the central region. Notably, the S
summation of N e
is consistent with the total number of electrons in the entire S
system, i.e., integer, although N e
is a non-integer.
8
2.2. DC-CCSD method This subsection briefly summarizes the conventional DC-CCSD method. As the DC-HF method yields the MOs C and the orbital energies of the subsystem α defined by L( , the local CCSD calculations can be performed for the localization region. The conventional DC-CCSD method approximately separates the occupied orbitals { i , j , …} and virtual orbitals { a , b , …} by replacing the Fermi distribution function f F p with the step function F p . The DC-based correlation energy is estimated by summing the correlation energies corresponding to the individual subsystems. To avoid double counting of the correlation energy owing to the overlap of the buffer regions, EDA is adopted to extract the correlation energy for the central region. DC Ecorr Ecorr .
S
(11)
S Ecorr
where
C S
i , j a ,b
* i
j a b 2ti j ,a b ti j ,b a ,
(12)
pq rs represents a two-electron integral notation. For the DC-CCSD method,
the two-electron excitation coefficient for the subsystem, t , becomes i j ,a b
ti j ,a b ti j ,a b ti ,a t j ,b ,
(13)
where t and t are the single and double excitation amplitudes of the local i j ,a b i ,a
9
CCSD, respectively. The local DM in this treatment is given as follows: D 2 F p Cp*Cp 2 Ci*Ci . occ
p
(14)
i
Thus, the number of electrons in the localization region L( is evaluated as
NeL =
L
Note that N e
D
, L
S .
(15)
is an integer. Furthermore, the number of electrons corresponding to
the central region S ( is given in a similar way to Eq. (10)
NeS =
P
, L
D S
S
Notably, the summation of N e
D
S L
S .
(16)
over the subsystems is not consistent with the total
number of electrons in the entire system
Ne NeS .
(17)
2.3. FON-DC-CCSD method This subsection presents a novel DC-CCSD method by adopting the FON formalism, i.e., FON-DC-CCSD. This formulation aims to address the inconsistency of the orbital occupations of the subsystem between the DC-HF and the conventional DC-CCSD methods—namely, the former treats the FON by utilizing the Fermi distribution function, whereas the latter approximately uses the integer occupation
10
number with the step function. Table 1 shows a comparison between the conventional and FON-type CCSD methods. In the FON formalism, arbitrary orbitals {p, q, r, s, …} possess the FONs of electrons and holes given by
f p f F p ,
(18)
f p 1 f F p .
(19)
Thus, the occupied and virtual orbital spaces are no longer separable. The excitation operators are defined as
ˆp ,q aq†a p ,
(20)
ˆpq ,rs ar†as†aq a p ,
(21)
where a †p and a p represent the creation and annihilation operators, respectively. The non-zero commutation and anti-commutation relations become
a†p , aq f p pq ,
(22)
a p , aq† f p pq .
(23)
Consequently, the thermal Wick theorem is given by the following non-zero contractions:
a†p aq f p pq ,
11
(24)
a p aq† f p pq .
(25)
With this theorem, the energy and amplitude equations for the FON-CCSD are derived as shown in Refs. 36 and 37. As in the case of the conventional DC-CCSD method, the correlation energy of the FON-DC-CCSD approach is estimated by summing the correlation energies for the central region using EDA. Instead of Eq. (12), the correlation energy in the case of the FON-DC-CCSD is evaluated as S Ecorr
S
p
Cp* q r s f p f q f r f s
q r s
2t
f p f F p ,
p q ,r s
(26)
t p q ,s r ,
f p 1 f F p ,
(27)
where indices { p , q , r , s } refer to the arbitrary orbitals of the subsystem α. As the numbers of electrons in L( and S ( are evaluated using Eqs. (8) and (10), respectively, the FON-DC-CCSD method satisfies the constraint on the total number of electrons and the DC-HF method. A comparison between Eqs. (12) and (26) reveals that the computational cost
L 2 L increases significantly from O Nocc N vir 2
L 4 to O N orb , where N occ , N vir , and
represent the numbers of occupied, virtual, and arbitrary orbitals of the subsystem N orb
α, respectively. If the FONs of electrons and holes are close to zero, the contribution to
12
the correlation energy in Eq. (24) is negligible. Therefore, by introducing a threshold λ, we define the quasi-occupied and quasi-virtual orbitals that satisfy f p and f p , respectively. The present study adopts 1015 . The numbers of
quasi-occupied and quasi-virtual orbitals are represented as
and N qvir , N qocc
respectively. Although the relations Nocc and N vir are Nqocc Norb Nqvir Norb
satisfied, Nocc and N vir are established in practical cases, as shown in Nqocc Nqvir
the following section. Consequently, the correlation energy for the central region in the FON-DC-CCSD case becomes
Ecorr
Nqocc Nqvir
C S
i
* i
j a b fi f j f a fb 2ti j ,a b ti j ,b a .
, j a ,b
(28)
The aforementioned FON-DC-CCSD method was implemented in the house-code of the GAMESS program [38].
13
3. Illustrative applications 3.1. Polyene chain In this section, the results of a numerical examination of the performance of the FON-DC-CCSD method in comparison with the conventional DC-CCSD calculations are presented. The examined systems were bond-alternating (BA) and the uniform (U) polyene chains, C40H42, are as shown in Fig. 2. All the C-C-H and C-C-C bond angles were fixed at 120°, and the C-H bond length was set at 1.07 Å. For the BA chain, the single-bond C-C and double-bond C=C lengths were set as 1.46 and 1.36 Å, respectively, whereas for the U chain, a constant C-C length of 1.40 Å was adopted. The utilized basis set was 6-31G** for C and H [40,41]. For all the CCSD calculations, MOs corresponding to the C 1s orbitals were treated as frozen cores. For the DC calculations, the temperature parameter of the Fermi function was fixed to 1578 K, which was a default value in the GAMESS program [38]. Individual C2H2(3) units were used as central regions, which were labeled from the left edge as shown in Fig. 2. The six units on the left and right sides were treated as the HF buffer region for the DC-HF calculations. By applying the dual-buffer scheme [39], the electron correlation was treated with different buffer sizes, namely, from one to six units on the left and right sides. In the dual-buffer scheme, we redetermined the Fermi level to conserve the total
14
number of electrons. Although the DM obtained using DC-CCSD is different from that obtained using DC-HF, the number of electrons within the localization region is conserved through the DC-correlation calculations. Table 2 shows the numbers of electrons in the localization and central regions obtained using the conventional and FON-type DC-CCSD calculations with a buffer size nb = 2, namely, { N eL , N eS } and { N eL , N eS }, respectively. Table 2 numerically confirms that N eL becomes an integer as theoretically predicted, although N eL is a non-integer. However, in the case of the BA chain, the differences between N eL and N eL are small, i.e., less than 0.218. The U chain leads to large differences between N eL and N eL , i.e., a maximum of 1.965. This is due to the delocalization nature of the electronic structure of the U chain. The formal numbers of electrons in the central regions, C2H2 for #2-#19 and C2H3 for #1 and #20, are 14 and 15, respectively. Although both N eS and N eS are non-integers, the corresponding values for the BA chain are considerably closer to the formal number of electrons, i.e., the differences are 0.027 for N eS and 0.014 for N eS at the maximum. In the case of the U chain, the discrepancies are in the range
-0.054 to 0.282 for N eS and -0.155 to 0.095 for N eS . The summation of N eS is consistent with the total number of electrons in the entire system for both the BA and U
15
chains, which satisfies the theoretical condition. However, the summation of N eS is not consistent with these results. In particular, in the case of U chain, the difference is 3.054. Figure 3 shows the buffer-size dependence of the conventional and FON-type DC-CCSD energy errors of C40H42 from the standard CCSD energy. In all cases, the errors decrease for larger buffer sizes. For the BA polyene chain, the errors of the conventional and FON-type DC-CCSD exhibit a fast convergence to zero: for nb ≥ 2, 1 kcal/mol or less, i.e., the so-called chemical accuracy. The difference in the energy errors between the conventional and FON-type DC-CCSD is small except for nb = 1. For the U chain, the errors of conventional DC-CCSD calculations slowly converge to zero, whereas the errors of the FON-DC-CCSD results show a fast convergence. Table 3 shows the buffer-size dependence of the numbers of (quasi-)occupied and (quasi-)virtual orbitals at subsystem #10, which is located in the middle of the polyene chain. Although N occ and N vir in DC-CCSD are well defined, N qocc and N qvir in
FON-DC-CCSD change owing to the FONs. The increase from N occ to N qocc as well
as from N vir to N qvir results in an increase in the computational cost. However, in the
polyene chains, the increases are extremely small: namely, {+1, +1} and {+3, +3} for { Nqocc , Nqvir Nocc N vir } of the BA and U chains, respectively. In Table 3, the number
16
of CCSD iterations, Niter , which is also related to the computational cost, is shown for
each case. For the BA chain, Niter in FON-DC-CCSD is close to that in DC-CCSD.
However, Niter in FON-DC-CCSD is considerably smaller than that in DC-CCSD. This
indicates that the FON treatment improves the poor convergence owing to the inconsistency of the orbital occupations of the subsystem between the DC-HF and DC-CCSD calculations, in particular, for delocalized systems. Finally, we examine the computational cost of the conventional and FON-type DC-CCSD calculations. Figure 4 shows the energy errors from the standard CCSD results versus the central processing unit (CPU) time for the conventional and FON-type DC-CCSD calculations of the BA and U polyene chains, using the CPU times over eight cores of two Intel Xeon X5680 (3.33 GHz) processors. Each plot corresponds to the case of nb = 1, 2, 3, 4, and 5. At the same buffer size, it is determined that the CPU times of the FON-DC-CCSD calculations are typically smaller than those of the conventional DC-CCSD. This is because the effect of the decrease in Niter is greater than that of the
increases in N qocc and N qvir in terms of the computational cost. Furthermore, at the
same accuracy, i.e., the same energy error, the effectiveness of FON-DC-CCSD is apparent, especially, for the delocalized U chain.
17
3.2. Single-walled carbon nanotubes This subsection examines the applicability of the FON-DC-CCSD method to single-walled carbon nanotubes (SWCNTs). SWCNTs have different characteristics depending on the chiral vectors Ch (n,m). SWCNTs are metallic if mod(n – m, 3) = 0, and semiconducting if mod(n – m, 3) = 1 or 2. The examined systems were metallic and semiconducting SWCNTs with chiral vectors (3,3) and (3,2), whose lengths corresponded to 13 benzene units, terminated by hydrogen atoms. It is difficult to accurately treat SWCNTs using standard fragmentation methods, because the aromaticity results in complicated subsystems. Generally, fragmentation with chemical intuition might achieve higher accuracy. However, in the present study, we adopted the automatic fragmentation technique, which defines the subsystem, i.e., the central region, as a cubic box with a certain length without any chemical intuition. We used 3.0 Å for the cubic box as illustrated in Fig. S1 in the supporting information. Based on the dual-buffer scheme, only the electron correlation was treated with the DC approach after the conventional HF calculations. The atoms within rb of each subsystem are treated as the correlation buffer region. The numbers of electrons in the localization and central regions calculated using the conventional and FON-type DC-CCSD methods are summarized in Table S1 in the supporting information. A similar tendency to that of the
18
polyene chains was observed, as shown in Table 2. Figure 5 shows the buffer-radius dependence rb [Å] of the conventional and FON-type DC-CCSD energy errors of the metallic and semiconducting SWCNTs from the standard CCSD energy. Although the subsystems are heterogeneous, the exhibited errors converge to zero for larger buffer sizes. For the semiconducting SWCNT, the difference in the energy errors between the conventional and FON-type DC-CCSD calculations is small. For the metallic SWCNTs, the FON-DC-CCSD calculations show a rapid decrease in the energy error compared with the conventional DC-CCSD calculations.
19
Conclusion In this study, we applied the FON formalism to address the inconsistency in the treatment of orbital occupations between the DC-HF and DC-CCSD methods, where the FONs of the subsystem orbitals are kept from the DC-HF with FON-DC-CCSD. Numerical applications to polyene chains and SWCNTs confirmed that the proposed FON-DC-CCSD method reduces the energy errors from the standard CCSD results in comparison with the conventional DC-CCSD method. Consequently, smaller buffer sizes can be adopted in the case of the FON-DC-CCSD method. Although the occupied and virtual spaces slightly increase, a difference in the maximum of three orbitals is observed through an examination of the polyene chains. The computational times for the proposed FON-DC-CCSD calculations are less than those for the conventional DC-CCSD calculations, because the numbers of CCSD iteration are considerably reduced. The aforementioned characteristics are more apparent for systems with delocalized electronic structures, for which difficulties are generally encountered when applying fragmentation methods. Consequently, the proposed FON-DC-CCSD method can expand the applicability of fragmentation methods.
20
Acknowledgments Some of the presented calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, Institutes of Natural Sciences (NINS). This study was supported by the program “Elements Strategy Initiative to Form Core Research Center” of the Ministry of Education, Culture, Sports, Science and Technology (MEXT).
21
References [1]
K. Raghavachari, G. W. Trucks, J. A. Pople, M. Head-Gordon, Chem. Phys. Lett. 157 (1989) 479.
[2]
T. J. Lee, G. E. Scuseria, Quantum Mechanical Electronic Structure Calculations with Chemical Accuracy, Kluwer Academic Publishers: Dordrecht (1995).
[3]
H. Stoll, Chem. Phys. Lett. 191 (1992) 548.
[4]
K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayashi, Chem. Phys. Lett. 313 (1999) 701.
[5]
S. Li, J. Ma, Y. Jiang, J. Comput. Chem. 23 (2002) 237.
[6]
D. W. Zhang, J. Z. H. Zhang, J. Chem. Phys. 119 (2003) 3599.
[7]
K. Babu, S. R. Gadre, J. Comput. Chem. 24 (2003) 484.
[8]
V. Deev, M. A. Collins, J. Chem. Phys. 112 (2005) 154102.
[9]
W. Li, S. Li, Y. S. Jiang, J. Phys. Chem. A 111 (2007) 2193.
[10] X. He, B. Wang, K. M. Merz, J. Phys. Chem. B 113 (2009) 10380. [11] N. J. Mayhall, K. Raghavachari, J. Chem. Theor. Comput. 7 (2011) 1336. [12] H. A. Le, H. J. Tan, J. F. Ouyang, R. P. A. Bettens, J. Chem. Theor. Comput. 8 (2012) 469. [13] R. M. Richard, J. M. Herbert, J. Chem. Phys. 137 (2012) 064113.
22
[14] J. J. Eriksen, P. Baudin, P. Ettenhuber, K. Kristensen, T. Kjærgaard, P. Jørgensen, J. Chem. Theor. Comput. 11 (2015) 2984. [15] W. Yang, Phys. Rev. Lett. 66 (1991) 1438. [16] W. Yang, T. S. Lee, J. Chem. Phys. 103 (1995) 5674. [17] T. E. Exner, P. G. Mezey, J. Phys. Chem. A 108 (2004) 4301. [18] X. He, J. Z. H. Zhang, J. Chem. Phys. 122 (2005) 031103. [19] W. Li, S. Li, J. Chem. Phys. 122 (2005) 194109. [20] T. Akama, M. Kobayashi, H. Nakai, J. Comput. Chem. 28 (2007) 2003. [21] T. Akama, M. Kobayashi, H. Nakai, Int. J. Quantum Chem. 109 (2009) 2706. [22] T. Akama, A. Fujii, M. Kobayashi, H. Nakai, Mol. Phys. 105 (2010) 2799. [23] M. Kobayashi, Y. Imamura, H. Nakai, J. Chem. Phys. 127 (2007) 074103. [24] M. Kobayashi, H. Nakai, J. Chem. Phys. 129 (2008) 044103. [25] M. Kobayashi, H. Nakai, J. Chem. Phys. 131 (2009) 114108. [26] T. Yoshikawa, M. Kobayashi, H. Nakai, Theor. Chem. Acc. 130 (2011) 411. [27] T. Yoshikawa, M. Kobayashi, H. Nakai, Int. J. Quant. Chem. 113 (2013) 218. [28] J. Seino, H. Nakai, J. Chem. Phys. 139 (2013) 034109. [29] M. Nakano, T. Yoshikawa, S. Hirata, J. Seino, H. Nakai, J. Comput. Chem. 38 (2017) 2520.
23
[30] M. S. Gordon, D. G. Fedorov, S. R. Pruitt, L. V. Slipchenko, Chem. Rev. 112 (2012) 632. [31] H. Nakai, Chem. Phys. Lett. 363 (2002) 73. [32] T. Matsubara, Theor. Phys. 14 (1955) 351. [33] D. J. Thouless, Phys. Rev. 107 (1957) 1162. [34] G. Sanyal, S. H. Mandai, S. Guha, D. Mukherjee, Phys. Rev. E 48 (1993) 3373. [35] M. Kobayashi, T. Taketsugu, Theor. Chem. Acc. 134 (2015) 107. [36] X. He, S. Ryu, S. Hirata, J. Chem. Phys. 140 (2014) 024702. [37] S. Hirata, X. He, J. Chem Phys. 138 (2013) 204112. [38] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery Jr., J. Comput. Chem. 14 (1993) 1347. [39] M. Kobayashi, H. Nakai, Int. Quant. Chem. 109 (2009) 2227. [40] W. J. Hehre, R. Ditchfield, J. A. Pople, J. Chem. Phys. 56 (1972) 2257. [41] J. D. Dill, J. A. Pople, J. Chem. Phys. 62 (1975) 2921.
24
Figure Captions Figure 1. Schematics of (a) non-integer-type and (b) integer-type occupations in each subsystem
Figure 2. Labeling of central regions of the BA and U polyene chains C40H42
Figure 3. Buffer-size (nb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD results of BA and U polyene chains C40H42
Figure 4. Conventional and FON-type DC-CCSD energies errors from standard CCSD energy of BA and U polyene chains C40H42 versus CPU times
Figure 5. Buffer-radius (rb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD calculations of semiconducting (3,2) and metallic (3,3) SWCNTs
25
Table 1. Comparison between conventional and FON-type DC-CCSD methods based on the canonical HF reference
DC-CCSD Occupied orbital {i,j,…} Virtual orbital {a,b,…} 1 0
FON-DC-CCSD Arbitrary orbital {p,q,…}
Excitation operator
ˆi ,a aa†ai ˆij ,ab aa†ab†a j ai
ˆp ,q aq†a p ˆpq ,rs ar†as†aq a p
Anti-commutation relation
ai† , a j a j , ai† ij † † aa , ab ab , aa ab
a†p , aq f p pq , a p , aq† f p pq
ai† a j ij
a †p aq f p pq
Orbital Occupation
Contraction CCSD energy of subsystem
S Ecor r
f p 1 f F p
aa ab† ab
C S
f p f F p
i
, j a ,b
* i
j a b 2ti j ,a b ti j ,b a
26
EcSorr
a p aq† f p pq
S
p
,q ,r , s
Cp* q r s
2t
p q ,r s
t p q ,s r
Table 2. Numbers of electrons of the localization and central regions calculated based on conventional and FON-type DC-CCSD calculations of BA and U polyene chains C40H42
DC-CCSD Central
BA
FON-DC-CCSD U
region
N eL
N eS
N eL
#1, #20
44.000
15.001
44.000
#2, #19
58.000
13.976
#3, #18
72.000
#4, #17
BA
N eS
U
N eL
N eS
N eL
N eS
15.282
44.000
15.000
44.000
14.888
58.000
13.946
58.000
13.977
58.000
13.845
13.983
72.000
14.026
72.001
13.983
72.000
13.882
72.000
13.975
72.000
14.104
72.064
13.987
72.041
13.891
#5, #16
72.000
13.974
72.000
14.161
72.129
13.998
73.169
14.025
#6, #15
72.000
13.973
72.000
14.188
72.171
14.005
73.848
14.090
#7, #14
72.000
13.973
72.000
14.200
72.195
14.010
73.935
14.095
#8, #13
72.000
13.973
72.000
14.205
72.208
14.012
73.956
14.095
#9, #12
72.000
13.973
72.000
14.207
72.215
14.013
73.963
14.095
#10, #11
72.000
13.973
72.000
14.208
72.218
14.014
73.965
14.095
Sum
-
281.548
-
285.054
-
282.000
27
-
282.000
Table 3. Buffer-size (nb) dependence of the numbers of (quasi-)occupied orbitals and (quasi-)virtual orbitals for conventional and FON-type DC-CCSD calculations of BA and U polyene chains C40H42. The numbers of CCSD iterations for the determination of excitation amplitudes, Niter, in conventional and FT-type DC-CCSD calculations are also listed.
nb
DC-CCSD
FON-DC-CCSD
BA N occ
N vir
U Niter
N occ
N vir
BA Niter
N qocc
N qvir
U Niter
N qocc
N qvir
Niter
1
16
98
53
16
98
183
17
99
52
19
101
92
2
26
164
42
26
164
163
27
165
42
29
167
41
3
36
230
41
36
230
112
36
230
42
39
233
42
4
46
296
41
46
296
75
46
296
41
49
299
41
5
56
362
41
56
362
74
56
362
41
59
365
42
28
Figure 1. Schematics of (a) non-integer-type and (b) integer-type occupations in each subsystem
29
Figure 2. Labeling of central regions of the BA and U polyene chains C40H42
30
Figure 3. Buffer-size (nb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD results of BA and U polyene chains C40H42
31
Figure 4. Conventional and FON-type DC-CCSD energies errors from standard CCSD energy of BA and U polyene chains C40H42 versus CPU times
32
Figure 5. Buffer-radius (rb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD calculations of semiconducting (3,2) and metallic (3,3) SWCNTs
33
34
Highlights
The fractional-occupation-number based divide-and-conquer coupled-cluster (FT-DC-CC) is proposed.
The FON formalism solves the inconsistency of occupations for the DC-CC method.
The FON-DC-CC method reduces the energy errors from the DC-CC method.
FON-DC-CC reduces the computational costs by adopting smaller buffer sizes.
35
S0009-2614(18)30787-5 https://doi.org/10.1016/j.cplett.2018.09.056 CPLETT 35968
To appear in:
Chemical Physics Letters
Received Date: Accepted Date:
4 August 2018 24 September 2018
Please cite this article as: T. Yoshikawa, H. Nakai, Fractional-occupation-number based divide-and-conquer coupled-cluster theory, Chemical Physics Letters (2018), doi: https://doi.org/10.1016/j.cplett.2018.09.056
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Fractional-occupation-number based divide-and-conquer coupled-cluster theory
Takeshi Yoshikawaa, Hiromi Nakaia-c *
a
Department of Chemistry and Biochemistry, School of Advanced Science and
Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555, Japan b
Research Institute for Science and Engineering, Waseda University, 3-4-1 Okubo,
Shinjuku-ku, Tokyo 169-8555, Japan c
Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University,
Katsura, Kyoto 615-8520, Japan _____________________________________________________________________ *
Corresponding author. FAX: +81-3-3205-2504, E-mail address: [email protected] (H. Nakai). URL: http://www.chem.waseda.ac.jp/nakai/ (H. Nakai).
1
Abstract We have extended the divide-and-conquer (DC) coupled-cluster with singles and doubles (CCSD) to a fractional occupation number (FON) formalism, denoted as FON-DC-CCSD, using the thermal Wick theorem. The motivation is to address the inconsistency in the treatment of orbital occupations between the DC-based Hartree– Fock and the DC-CCSD methods, which adopt the Fermi distribution function and the step function for orbital occupation, respectively. Numerical applications involving polyene chains and single-walled carbon nanotubes confirm that the proposed FON-DC-CCSD method reduces both energy errors and computational costs compared with the conventional DC-CCSD method.
2
1. Introduction The coupled-cluster (CC) theory is a robust theoretical approach that systematically improves the accuracy of calculation of the total energy and the wavefunction of atoms and molecules by considering the increase in the order of electron excitation. Based on the correlated method, the energetics of various molecules can be predicted with chemical accuracy. The CC method with single, double, and perturbative triple excitations (CCSD(T)) has been widely used in practical applications and benchmark studies over the years [1,2]. However, its high scaling with respect to the system size N remains a severe impediment for large systems. For example, the computational costs of CCSD and CCSD(T) are asymptotically O(N6) and O(N7), respectively. Many researchers have reported a wide variety of linear-scaling methods based on fragmentation approaches [3-19], which partition the system under consideration into several subsystems. Li et al. [9] categorized the fragmentation approaches into energy-based approaches [3-14] and density matrix (DM)-based approaches [15-29]. Energy-based approaches directly evaluate the total energy of the system from subsystem contributions. The approach is generally applicable to various theory levels such as Hartree–Fock (HF), Møller–Plesset perturbation, and CC theories. In the
3
DM-based approaches, the total energy is calculated from the assembled DM, which is constructed from the DMs of a series of subsystems. One of the features of the DM-based approaches is that a non-integer number of electrons can be treated for each subsystem. Thus, it is possible to reasonably treat charge- and/or spin-delocalized systems using the DM-based approach. Gordon et al. [30] recognized one-step and two-step techniques for the energy-based and DM-based approaches, respectively. The divide-and-conquer (DC) method proposed by Yang and Lee [15,16] can be directly applied to the density functional theory, and is an example of a DM-based approach. These authors investigated the performance of this approach for calculations including HF exchange interactions [20-22]. Furthermore, DC-based correlation methods [23-29] were developed by exploiting energy density analysis (EDA) techniques [31]. The correlation energy of the total system is evaluated by summing up subsystem contributions. Therefore, the DC-based correlation methods are categorized as an energy-based approach. Note that the subsystem orbitals are obtained by diagonalizing the subsystem Fock matrices, which are constructed using the assembled DM. However, the occupation numbers of the subsystem orbitals need to be approximated from non-integer to integer, to solve the subsystem correlation problem. Recently, Kobayashi and Taketsugu proposed the DC-based second-order Møller–
4
Plesset perturbation (DC-MP2) method based on the thermal Wick theorem [32-34] and improved its accuracy [35]. In this study, we have extended the fractional occupation number (FON) technique to the DC-CCSD theory, denoted as FON-DC-CCSD. The rest of this report is organized as follows: in Section 2, we present the theoretical aspects; Section 3 describes the numerical applications; Section 4 presents the conclusion.
5
2. Theory 2.1. DC-HF method This subsection provides a brief explanation of the DC-HF method, the results of which are the basis for the DC-CCSD calculations described in the following subsection. For simplicity, we consider the case of closed-shell systems. The extension to open-shell systems is straightforward. In the DC-HF method, the system under consideration is spatially divided into disjoint subsystems, which are called the central regions. A set of atomic orbitals (AOs) corresponding to a central region α is denoted by S ( . To improve the description of the subsystem, the area adjacent to the central region, which is called the buffer region, is considered during the expansion of the subsystem’s molecular orbitals (MOs). A set of AOs corresponding to the buffer region of the subsystem α, which is denoted by B( , is added to S ( to construct a set of AOs in the localization region of the
subsystem α, L( , as follows: L( S( B .
(1)
The DC-HF method evaluates the MOs C and the orbital energies of the subsystem α by solving the following local Roothaan–Hall equation for the localization region:
6
F C S C ,
(2)
where Sα and Fα represent the local overlap and Fock matrices, respectively, which are extracted from the total matrices. The DM of the entire system, D, is constructed from the local DM for the subsystem α, Dα, as follows: . D P D
(3)
Pα represents the partition matrix with elements given by
1 P 1 2 0
S S B S S B
(4)
otherwise.
The local DM for the subsystem α, Dα, is obtained using the subsystem MOs Cα, the subsystem orbital energies ɛα, and the common Fermi level F D 2 f F p Cp*Cp ,
(5)
f x 1 exp x .
(6)
p
1
Here, indices {μ, ν} and {p, q} refer to AOs and MOs, respectively. Unlike the standard HF method, the DM is constructed on the basis of the Fermi distribution function fβ(x) for the inverse temperature kBT
1
instead of the step function θ(x). Note that the
Fermi distribution function is identical to the step function in the limit of β → ∞, i.e., T → 0 K. As the occupation numbers of the subsystem orbitals are non-integer as
7
shown in Fig. 1(a), it is possible to reasonably treat charge- and/or spin-delocalized systems. The common Fermi level F is determined by a constraint on the total number of electrons in the entire system, N e , as follows:
Ne D S .
(7)
,
Similarly, the number of electrons for the subsystem α, N eL , is determined by
NeL =
D
, L
S ,
(8)
whose value is a non-integer because of fractional occupation. Note that N e is not consistent with the sum N eL because of double counting of the number of electrons in the buffer regions. The partition matrix Pα in Eq. (4) avoids such double counting. From Eqs. (3) and (7), Ne P D S N eS .
,
(9)
Here, N eS is defined and rewritten as follows:
NeS =
S
Thus, N e
P
, L
D S
D
S L
S .
(10)
is regarded as the number of electrons in the central region. Notably, the S
summation of N e
is consistent with the total number of electrons in the entire S
system, i.e., integer, although N e
is a non-integer.
8
2.2. DC-CCSD method This subsection briefly summarizes the conventional DC-CCSD method. As the DC-HF method yields the MOs C and the orbital energies of the subsystem α defined by L( , the local CCSD calculations can be performed for the localization region. The conventional DC-CCSD method approximately separates the occupied orbitals { i , j , …} and virtual orbitals { a , b , …} by replacing the Fermi distribution function f F p with the step function F p . The DC-based correlation energy is estimated by summing the correlation energies corresponding to the individual subsystems. To avoid double counting of the correlation energy owing to the overlap of the buffer regions, EDA is adopted to extract the correlation energy for the central region. DC Ecorr Ecorr .
S
(11)
S Ecorr
where
C S
i , j a ,b
* i
j a b 2ti j ,a b ti j ,b a ,
(12)
pq rs represents a two-electron integral notation. For the DC-CCSD method,
the two-electron excitation coefficient for the subsystem, t , becomes i j ,a b
ti j ,a b ti j ,a b ti ,a t j ,b ,
(13)
where t and t are the single and double excitation amplitudes of the local i j ,a b i ,a
9
CCSD, respectively. The local DM in this treatment is given as follows: D 2 F p Cp*Cp 2 Ci*Ci . occ
p
(14)
i
Thus, the number of electrons in the localization region L( is evaluated as
NeL =
L
Note that N e
D
, L
S .
(15)
is an integer. Furthermore, the number of electrons corresponding to
the central region S ( is given in a similar way to Eq. (10)
NeS =
P
, L
D S
S
Notably, the summation of N e
D
S L
S .
(16)
over the subsystems is not consistent with the total
number of electrons in the entire system
Ne NeS .
(17)
2.3. FON-DC-CCSD method This subsection presents a novel DC-CCSD method by adopting the FON formalism, i.e., FON-DC-CCSD. This formulation aims to address the inconsistency of the orbital occupations of the subsystem between the DC-HF and the conventional DC-CCSD methods—namely, the former treats the FON by utilizing the Fermi distribution function, whereas the latter approximately uses the integer occupation
10
number with the step function. Table 1 shows a comparison between the conventional and FON-type CCSD methods. In the FON formalism, arbitrary orbitals {p, q, r, s, …} possess the FONs of electrons and holes given by
f p f F p ,
(18)
f p 1 f F p .
(19)
Thus, the occupied and virtual orbital spaces are no longer separable. The excitation operators are defined as
ˆp ,q aq†a p ,
(20)
ˆpq ,rs ar†as†aq a p ,
(21)
where a †p and a p represent the creation and annihilation operators, respectively. The non-zero commutation and anti-commutation relations become
a†p , aq f p pq ,
(22)
a p , aq† f p pq .
(23)
Consequently, the thermal Wick theorem is given by the following non-zero contractions:
a†p aq f p pq ,
11
(24)
a p aq† f p pq .
(25)
With this theorem, the energy and amplitude equations for the FON-CCSD are derived as shown in Refs. 36 and 37. As in the case of the conventional DC-CCSD method, the correlation energy of the FON-DC-CCSD approach is estimated by summing the correlation energies for the central region using EDA. Instead of Eq. (12), the correlation energy in the case of the FON-DC-CCSD is evaluated as S Ecorr
S
p
Cp* q r s f p f q f r f s
q r s
2t
f p f F p ,
p q ,r s
(26)
t p q ,s r ,
f p 1 f F p ,
(27)
where indices { p , q , r , s } refer to the arbitrary orbitals of the subsystem α. As the numbers of electrons in L( and S ( are evaluated using Eqs. (8) and (10), respectively, the FON-DC-CCSD method satisfies the constraint on the total number of electrons and the DC-HF method. A comparison between Eqs. (12) and (26) reveals that the computational cost
L 2 L increases significantly from O Nocc N vir 2
L 4 to O N orb , where N occ , N vir , and
represent the numbers of occupied, virtual, and arbitrary orbitals of the subsystem N orb
α, respectively. If the FONs of electrons and holes are close to zero, the contribution to
12
the correlation energy in Eq. (24) is negligible. Therefore, by introducing a threshold λ, we define the quasi-occupied and quasi-virtual orbitals that satisfy f p and f p , respectively. The present study adopts 1015 . The numbers of
quasi-occupied and quasi-virtual orbitals are represented as
and N qvir , N qocc
respectively. Although the relations Nocc and N vir are Nqocc Norb Nqvir Norb
satisfied, Nocc and N vir are established in practical cases, as shown in Nqocc Nqvir
the following section. Consequently, the correlation energy for the central region in the FON-DC-CCSD case becomes
Ecorr
Nqocc Nqvir
C S
i
* i
j a b fi f j f a fb 2ti j ,a b ti j ,b a .
, j a ,b
(28)
The aforementioned FON-DC-CCSD method was implemented in the house-code of the GAMESS program [38].
13
3. Illustrative applications 3.1. Polyene chain In this section, the results of a numerical examination of the performance of the FON-DC-CCSD method in comparison with the conventional DC-CCSD calculations are presented. The examined systems were bond-alternating (BA) and the uniform (U) polyene chains, C40H42, are as shown in Fig. 2. All the C-C-H and C-C-C bond angles were fixed at 120°, and the C-H bond length was set at 1.07 Å. For the BA chain, the single-bond C-C and double-bond C=C lengths were set as 1.46 and 1.36 Å, respectively, whereas for the U chain, a constant C-C length of 1.40 Å was adopted. The utilized basis set was 6-31G** for C and H [40,41]. For all the CCSD calculations, MOs corresponding to the C 1s orbitals were treated as frozen cores. For the DC calculations, the temperature parameter of the Fermi function was fixed to 1578 K, which was a default value in the GAMESS program [38]. Individual C2H2(3) units were used as central regions, which were labeled from the left edge as shown in Fig. 2. The six units on the left and right sides were treated as the HF buffer region for the DC-HF calculations. By applying the dual-buffer scheme [39], the electron correlation was treated with different buffer sizes, namely, from one to six units on the left and right sides. In the dual-buffer scheme, we redetermined the Fermi level to conserve the total
14
number of electrons. Although the DM obtained using DC-CCSD is different from that obtained using DC-HF, the number of electrons within the localization region is conserved through the DC-correlation calculations. Table 2 shows the numbers of electrons in the localization and central regions obtained using the conventional and FON-type DC-CCSD calculations with a buffer size nb = 2, namely, { N eL , N eS } and { N eL , N eS }, respectively. Table 2 numerically confirms that N eL becomes an integer as theoretically predicted, although N eL is a non-integer. However, in the case of the BA chain, the differences between N eL and N eL are small, i.e., less than 0.218. The U chain leads to large differences between N eL and N eL , i.e., a maximum of 1.965. This is due to the delocalization nature of the electronic structure of the U chain. The formal numbers of electrons in the central regions, C2H2 for #2-#19 and C2H3 for #1 and #20, are 14 and 15, respectively. Although both N eS and N eS are non-integers, the corresponding values for the BA chain are considerably closer to the formal number of electrons, i.e., the differences are 0.027 for N eS and 0.014 for N eS at the maximum. In the case of the U chain, the discrepancies are in the range
-0.054 to 0.282 for N eS and -0.155 to 0.095 for N eS . The summation of N eS is consistent with the total number of electrons in the entire system for both the BA and U
15
chains, which satisfies the theoretical condition. However, the summation of N eS is not consistent with these results. In particular, in the case of U chain, the difference is 3.054. Figure 3 shows the buffer-size dependence of the conventional and FON-type DC-CCSD energy errors of C40H42 from the standard CCSD energy. In all cases, the errors decrease for larger buffer sizes. For the BA polyene chain, the errors of the conventional and FON-type DC-CCSD exhibit a fast convergence to zero: for nb ≥ 2, 1 kcal/mol or less, i.e., the so-called chemical accuracy. The difference in the energy errors between the conventional and FON-type DC-CCSD is small except for nb = 1. For the U chain, the errors of conventional DC-CCSD calculations slowly converge to zero, whereas the errors of the FON-DC-CCSD results show a fast convergence. Table 3 shows the buffer-size dependence of the numbers of (quasi-)occupied and (quasi-)virtual orbitals at subsystem #10, which is located in the middle of the polyene chain. Although N occ and N vir in DC-CCSD are well defined, N qocc and N qvir in
FON-DC-CCSD change owing to the FONs. The increase from N occ to N qocc as well
as from N vir to N qvir results in an increase in the computational cost. However, in the
polyene chains, the increases are extremely small: namely, {+1, +1} and {+3, +3} for { Nqocc , Nqvir Nocc N vir } of the BA and U chains, respectively. In Table 3, the number
16
of CCSD iterations, Niter , which is also related to the computational cost, is shown for
each case. For the BA chain, Niter in FON-DC-CCSD is close to that in DC-CCSD.
However, Niter in FON-DC-CCSD is considerably smaller than that in DC-CCSD. This
indicates that the FON treatment improves the poor convergence owing to the inconsistency of the orbital occupations of the subsystem between the DC-HF and DC-CCSD calculations, in particular, for delocalized systems. Finally, we examine the computational cost of the conventional and FON-type DC-CCSD calculations. Figure 4 shows the energy errors from the standard CCSD results versus the central processing unit (CPU) time for the conventional and FON-type DC-CCSD calculations of the BA and U polyene chains, using the CPU times over eight cores of two Intel Xeon X5680 (3.33 GHz) processors. Each plot corresponds to the case of nb = 1, 2, 3, 4, and 5. At the same buffer size, it is determined that the CPU times of the FON-DC-CCSD calculations are typically smaller than those of the conventional DC-CCSD. This is because the effect of the decrease in Niter is greater than that of the
increases in N qocc and N qvir in terms of the computational cost. Furthermore, at the
same accuracy, i.e., the same energy error, the effectiveness of FON-DC-CCSD is apparent, especially, for the delocalized U chain.
17
3.2. Single-walled carbon nanotubes This subsection examines the applicability of the FON-DC-CCSD method to single-walled carbon nanotubes (SWCNTs). SWCNTs have different characteristics depending on the chiral vectors Ch (n,m). SWCNTs are metallic if mod(n – m, 3) = 0, and semiconducting if mod(n – m, 3) = 1 or 2. The examined systems were metallic and semiconducting SWCNTs with chiral vectors (3,3) and (3,2), whose lengths corresponded to 13 benzene units, terminated by hydrogen atoms. It is difficult to accurately treat SWCNTs using standard fragmentation methods, because the aromaticity results in complicated subsystems. Generally, fragmentation with chemical intuition might achieve higher accuracy. However, in the present study, we adopted the automatic fragmentation technique, which defines the subsystem, i.e., the central region, as a cubic box with a certain length without any chemical intuition. We used 3.0 Å for the cubic box as illustrated in Fig. S1 in the supporting information. Based on the dual-buffer scheme, only the electron correlation was treated with the DC approach after the conventional HF calculations. The atoms within rb of each subsystem are treated as the correlation buffer region. The numbers of electrons in the localization and central regions calculated using the conventional and FON-type DC-CCSD methods are summarized in Table S1 in the supporting information. A similar tendency to that of the
18
polyene chains was observed, as shown in Table 2. Figure 5 shows the buffer-radius dependence rb [Å] of the conventional and FON-type DC-CCSD energy errors of the metallic and semiconducting SWCNTs from the standard CCSD energy. Although the subsystems are heterogeneous, the exhibited errors converge to zero for larger buffer sizes. For the semiconducting SWCNT, the difference in the energy errors between the conventional and FON-type DC-CCSD calculations is small. For the metallic SWCNTs, the FON-DC-CCSD calculations show a rapid decrease in the energy error compared with the conventional DC-CCSD calculations.
19
Conclusion In this study, we applied the FON formalism to address the inconsistency in the treatment of orbital occupations between the DC-HF and DC-CCSD methods, where the FONs of the subsystem orbitals are kept from the DC-HF with FON-DC-CCSD. Numerical applications to polyene chains and SWCNTs confirmed that the proposed FON-DC-CCSD method reduces the energy errors from the standard CCSD results in comparison with the conventional DC-CCSD method. Consequently, smaller buffer sizes can be adopted in the case of the FON-DC-CCSD method. Although the occupied and virtual spaces slightly increase, a difference in the maximum of three orbitals is observed through an examination of the polyene chains. The computational times for the proposed FON-DC-CCSD calculations are less than those for the conventional DC-CCSD calculations, because the numbers of CCSD iteration are considerably reduced. The aforementioned characteristics are more apparent for systems with delocalized electronic structures, for which difficulties are generally encountered when applying fragmentation methods. Consequently, the proposed FON-DC-CCSD method can expand the applicability of fragmentation methods.
20
Acknowledgments Some of the presented calculations were performed at the Research Center for Computational Science (RCCS), Okazaki Research Facilities, Institutes of Natural Sciences (NINS). This study was supported by the program “Elements Strategy Initiative to Form Core Research Center” of the Ministry of Education, Culture, Sports, Science and Technology (MEXT).
21
References [1]
K. Raghavachari, G. W. Trucks, J. A. Pople, M. Head-Gordon, Chem. Phys. Lett. 157 (1989) 479.
[2]
T. J. Lee, G. E. Scuseria, Quantum Mechanical Electronic Structure Calculations with Chemical Accuracy, Kluwer Academic Publishers: Dordrecht (1995).
[3]
H. Stoll, Chem. Phys. Lett. 191 (1992) 548.
[4]
K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayashi, Chem. Phys. Lett. 313 (1999) 701.
[5]
S. Li, J. Ma, Y. Jiang, J. Comput. Chem. 23 (2002) 237.
[6]
D. W. Zhang, J. Z. H. Zhang, J. Chem. Phys. 119 (2003) 3599.
[7]
K. Babu, S. R. Gadre, J. Comput. Chem. 24 (2003) 484.
[8]
V. Deev, M. A. Collins, J. Chem. Phys. 112 (2005) 154102.
[9]
W. Li, S. Li, Y. S. Jiang, J. Phys. Chem. A 111 (2007) 2193.
[10] X. He, B. Wang, K. M. Merz, J. Phys. Chem. B 113 (2009) 10380. [11] N. J. Mayhall, K. Raghavachari, J. Chem. Theor. Comput. 7 (2011) 1336. [12] H. A. Le, H. J. Tan, J. F. Ouyang, R. P. A. Bettens, J. Chem. Theor. Comput. 8 (2012) 469. [13] R. M. Richard, J. M. Herbert, J. Chem. Phys. 137 (2012) 064113.
22
[14] J. J. Eriksen, P. Baudin, P. Ettenhuber, K. Kristensen, T. Kjærgaard, P. Jørgensen, J. Chem. Theor. Comput. 11 (2015) 2984. [15] W. Yang, Phys. Rev. Lett. 66 (1991) 1438. [16] W. Yang, T. S. Lee, J. Chem. Phys. 103 (1995) 5674. [17] T. E. Exner, P. G. Mezey, J. Phys. Chem. A 108 (2004) 4301. [18] X. He, J. Z. H. Zhang, J. Chem. Phys. 122 (2005) 031103. [19] W. Li, S. Li, J. Chem. Phys. 122 (2005) 194109. [20] T. Akama, M. Kobayashi, H. Nakai, J. Comput. Chem. 28 (2007) 2003. [21] T. Akama, M. Kobayashi, H. Nakai, Int. J. Quantum Chem. 109 (2009) 2706. [22] T. Akama, A. Fujii, M. Kobayashi, H. Nakai, Mol. Phys. 105 (2010) 2799. [23] M. Kobayashi, Y. Imamura, H. Nakai, J. Chem. Phys. 127 (2007) 074103. [24] M. Kobayashi, H. Nakai, J. Chem. Phys. 129 (2008) 044103. [25] M. Kobayashi, H. Nakai, J. Chem. Phys. 131 (2009) 114108. [26] T. Yoshikawa, M. Kobayashi, H. Nakai, Theor. Chem. Acc. 130 (2011) 411. [27] T. Yoshikawa, M. Kobayashi, H. Nakai, Int. J. Quant. Chem. 113 (2013) 218. [28] J. Seino, H. Nakai, J. Chem. Phys. 139 (2013) 034109. [29] M. Nakano, T. Yoshikawa, S. Hirata, J. Seino, H. Nakai, J. Comput. Chem. 38 (2017) 2520.
23
[30] M. S. Gordon, D. G. Fedorov, S. R. Pruitt, L. V. Slipchenko, Chem. Rev. 112 (2012) 632. [31] H. Nakai, Chem. Phys. Lett. 363 (2002) 73. [32] T. Matsubara, Theor. Phys. 14 (1955) 351. [33] D. J. Thouless, Phys. Rev. 107 (1957) 1162. [34] G. Sanyal, S. H. Mandai, S. Guha, D. Mukherjee, Phys. Rev. E 48 (1993) 3373. [35] M. Kobayashi, T. Taketsugu, Theor. Chem. Acc. 134 (2015) 107. [36] X. He, S. Ryu, S. Hirata, J. Chem. Phys. 140 (2014) 024702. [37] S. Hirata, X. He, J. Chem Phys. 138 (2013) 204112. [38] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery Jr., J. Comput. Chem. 14 (1993) 1347. [39] M. Kobayashi, H. Nakai, Int. Quant. Chem. 109 (2009) 2227. [40] W. J. Hehre, R. Ditchfield, J. A. Pople, J. Chem. Phys. 56 (1972) 2257. [41] J. D. Dill, J. A. Pople, J. Chem. Phys. 62 (1975) 2921.
24
Figure Captions Figure 1. Schematics of (a) non-integer-type and (b) integer-type occupations in each subsystem
Figure 2. Labeling of central regions of the BA and U polyene chains C40H42
Figure 3. Buffer-size (nb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD results of BA and U polyene chains C40H42
Figure 4. Conventional and FON-type DC-CCSD energies errors from standard CCSD energy of BA and U polyene chains C40H42 versus CPU times
Figure 5. Buffer-radius (rb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD calculations of semiconducting (3,2) and metallic (3,3) SWCNTs
25
Table 1. Comparison between conventional and FON-type DC-CCSD methods based on the canonical HF reference
DC-CCSD Occupied orbital {i,j,…} Virtual orbital {a,b,…} 1 0
FON-DC-CCSD Arbitrary orbital {p,q,…}
Excitation operator
ˆi ,a aa†ai ˆij ,ab aa†ab†a j ai
ˆp ,q aq†a p ˆpq ,rs ar†as†aq a p
Anti-commutation relation
ai† , a j a j , ai† ij † † aa , ab ab , aa ab
a†p , aq f p pq , a p , aq† f p pq
ai† a j ij
a †p aq f p pq
Orbital Occupation
Contraction CCSD energy of subsystem
S Ecor r
f p 1 f F p
aa ab† ab
C S
f p f F p
i
, j a ,b
* i
j a b 2ti j ,a b ti j ,b a
26
EcSorr
a p aq† f p pq
S
p
,q ,r , s
Cp* q r s
2t
p q ,r s
t p q ,s r
Table 2. Numbers of electrons of the localization and central regions calculated based on conventional and FON-type DC-CCSD calculations of BA and U polyene chains C40H42
DC-CCSD Central
BA
FON-DC-CCSD U
region
N eL
N eS
N eL
#1, #20
44.000
15.001
44.000
#2, #19
58.000
13.976
#3, #18
72.000
#4, #17
BA
N eS
U
N eL
N eS
N eL
N eS
15.282
44.000
15.000
44.000
14.888
58.000
13.946
58.000
13.977
58.000
13.845
13.983
72.000
14.026
72.001
13.983
72.000
13.882
72.000
13.975
72.000
14.104
72.064
13.987
72.041
13.891
#5, #16
72.000
13.974
72.000
14.161
72.129
13.998
73.169
14.025
#6, #15
72.000
13.973
72.000
14.188
72.171
14.005
73.848
14.090
#7, #14
72.000
13.973
72.000
14.200
72.195
14.010
73.935
14.095
#8, #13
72.000
13.973
72.000
14.205
72.208
14.012
73.956
14.095
#9, #12
72.000
13.973
72.000
14.207
72.215
14.013
73.963
14.095
#10, #11
72.000
13.973
72.000
14.208
72.218
14.014
73.965
14.095
Sum
-
281.548
-
285.054
-
282.000
27
-
282.000
Table 3. Buffer-size (nb) dependence of the numbers of (quasi-)occupied orbitals and (quasi-)virtual orbitals for conventional and FON-type DC-CCSD calculations of BA and U polyene chains C40H42. The numbers of CCSD iterations for the determination of excitation amplitudes, Niter, in conventional and FT-type DC-CCSD calculations are also listed.
nb
DC-CCSD
FON-DC-CCSD
BA N occ
N vir
U Niter
N occ
N vir
BA Niter
N qocc
N qvir
U Niter
N qocc
N qvir
Niter
1
16
98
53
16
98
183
17
99
52
19
101
92
2
26
164
42
26
164
163
27
165
42
29
167
41
3
36
230
41
36
230
112
36
230
42
39
233
42
4
46
296
41
46
296
75
46
296
41
49
299
41
5
56
362
41
56
362
74
56
362
41
59
365
42
28
Figure 1. Schematics of (a) non-integer-type and (b) integer-type occupations in each subsystem
29
Figure 2. Labeling of central regions of the BA and U polyene chains C40H42
30
Figure 3. Buffer-size (nb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD results of BA and U polyene chains C40H42
31
Figure 4. Conventional and FON-type DC-CCSD energies errors from standard CCSD energy of BA and U polyene chains C40H42 versus CPU times
32
Figure 5. Buffer-radius (rb) dependence of the conventional and FON-type DC-CCSD energy errors from the standard CCSD calculations of semiconducting (3,2) and metallic (3,3) SWCNTs
33
34
Highlights
The fractional-occupation-number based divide-and-conquer coupled-cluster (FT-DC-CC) is proposed.
The FON formalism solves the inconsistency of occupations for the DC-CC method.
The FON-DC-CC method reduces the energy errors from the DC-CC method.
FON-DC-CC reduces the computational costs by adopting smaller buffer sizes.
35