^{1}Division of Medicinal Safety Science, National Institute of Health Sciences, Kawasaki, Kanagawa 210-9501, Japan^{2}Department of Chemistry and Research Center for Smart Molecules, Faculty of Science, Rikkyo University, Toshima, Tokyo 171-8501, Japan^{3}Department of Physics, College of Science and Technology, Nihon University, Chiyoda, Tokyo 101-8308, Japan^{4}Graduate School of System Informatics, Kobe University, Kobe 657-8501, Japan
Received November 17, 2020; Accepted March 31, 2021; Published May 12, 2021
A dependable description of the frontier orbitals around the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) is essential for the theoretical elucidation of exciton and charge-transfer dynamics in molecular systems. In this study, an efficient scheme for constructing molecular orbitals (MOs) of large molecular systems is proposed within the framework of the fragment molecular orbital (FMO) method. In contrast to the previously proposed energy-based methodology, called the linear combination of molecular orbitals of FMO (FMO-LCMO), the present scheme employs the density matrices of monomer and dimer MOs of fragments to select the monomer MOs as the basis set for constructing the Hamiltonian matrix of the total system. The accuracy and computational cost of the proposed scheme are assessed in comparison with those of the conventional energy-based methods for three model systems of DNA stacked base pairs, pseudo-glycine pentamer, and water clusters, thus demonstrating its efficiency in terms of data compression.
This article is published by the Physical Society of Japan under the terms of the Creative Commons Attribution 4.0 License. Any further distribution of this work must maintain attribution to the author(s) and the title of the article, journal citation, and DOI.
A number of fragment-based methods have been proposed to conduct quantum chemical calculations of large molecular systems, including proteins, nucleic acids, and nanomaterials.^{1}^{)} One of the most successful examples of such attempts is the fragment molecular orbital (FMO) method,^{2}^{,}^{3}^{)} which applies fragmented molecular orbital (MO) calculations in the embedding electrostatic field of the entire system. The FMO method has been used primarily for calculations of the energies and gradients in geometry optimizations and molecular dynamics,^{4}^{,}^{5}^{)} and because the method is an energy-based fragmentation approach, obtaining the MOs and their energies for whole molecular systems is not easy. The simplest way to obtain the MOs and their energies in an FMO scheme is to take the union of fragment (monomer) MOs, which corresponds to the one-body expansion, FMO1. Pair corrections can then be added to the FMO-MO method,^{6}^{)} in which FMO provides an initial density for a couple of iterations in regular MO schemes, or in the linear combination of MOs of fragments (FMO-LCMO) method,^{7}^{)} where one can obtain a subset of orbitals and energies around the Fermi level or the frontier orbitals.^{8}^{–}^{14}^{)} In the present study, we focus on this FMO-LCMO method. Relevant descriptions of frontier orbitals are essential for understanding charge transfer and exciton dynamics in large molecular systems.^{15}^{,}^{16}^{)}
The FMO-LCMO method was originally proposed to obtain canonical MOs of large systems within the framework of the FMO scheme^{7}^{)} (see Sect. 2 for details). The FMO method first decomposes the total system into fragments and calculates the MOs of each fragment self-consistently under the electrostatic field of the other fragments. Dimer, trimer,^{17}^{,}^{18}^{)} or tetramer^{19}^{)} calculations are then carried out to account for the charge transfer and exchange interactions over multiple fragments. With the FMO-LCMO method, to remove the excess energy, the diagonal elements of the Hamiltonian matrix of the total system are computed from the results of these fragment calculations, whereas the off-diagonal elements are constructed from the dimer or trimer matrices projected to the monomer MO space.^{7}^{,}^{8}^{)} The total Hamiltonian matrix is then diagonalized to determine the canonical MOs of the entire system. It is thus supposed that the Hamiltonian matrix before diagonalization provides appropriate diagonal site energies and off-diagonal transfer integrals, including the exchange and charge-transfer effects.^{9}^{–}^{13}^{)} Although the FMO-LCMO method has been reported to attain a high accuracy at the FMO3 level,^{8}^{)} inaccuracies in the virtual orbitals have been considerable for some systems at the FMO2 level.^{20}^{)} The reason for this difficulty is likely to be the lack of higher-order exchange effects in the embedding potential. The accuracy of the FMO-LCMO method can then be assessed in terms of the orbital energies and density of states for a representative set of materials compared with full MO results without fragmentation.
The FMO-LCMO method is regarded as a computational scheme that can efficiently describe the MOs of large-scale molecular systems by using fragment monomer orbitals as the bases. In general, there is arbitrariness in the selection of the monomer MO bases; because interest is primarily in the canonical orbitals near the highest occupied molecular orbital (HOMO)/lowest unoccupied molecular orbital (LUMO), the HOMO/LUMO vicinities of the monomer MOs are uniformly selected as the bases. Therefore, our question is the following: Is there a basis selection scheme that most efficiently reproduces the canonical orbitals over this arbitrariness?
Thus, the current problem that should be resolved is how to optimize the selection of relevant fragment MOs for constructing the Hamiltonian matrix of the total system such that the physicochemical behaviors of the frontier orbitals are appropriately described with a relatively lower computational cost. Herein, we consider that one can construct a density matrix obtained by projecting a monomer MO onto a dimer MO. Each matrix element consists of an overlap integral between them, representing the contribution of a monomer MO to make up a dimer MO. To realize this concept, we diagonalize the density matrix and select the monomer MOs according to their eigenvalues and eigenvector components.
One possible solution is thus the diagonalization of the density matrix (referred to as TFDM below) between the monomer orbitals and the higher-order polymer orbitals. From the eigenvalues (occupation numbers) and eigenvectors (where their constituents are monomer orbitals), the monomer orbitals important for constituting the polymer orbitals are determined. That is, the idea here is that the monomer orbitals are superposed onto the higher-order polymer orbitals, and the monomer orbitals having a large overlap are essential for reproducing the higher-order polymer orbitals. In this study, we propose a new calculation scheme that extracts the effective basis space by efficiently compressing the unoptimized basis space, such as HOMO/LUMO neighborhoods that have been employed thus far.
In the following Sect. 2, we illustrate the basic theory and associated computational scheme of the present proposal. In Sect. 3, the test calculations for three molecular systems are shown to assess the usefulness of the present method. Some concluding remarks are given in Sect. 4.
2. Theory
In this section, we first review the fundamental formulation of the FMO method and the conventional monomer MO-based FMO-LCMO method. We next introduce a novel density-matrix-based scheme and related considerations. Furthermore, we show the computational levels and settings used in this study.
FMO-LCMO scheme
The FMO method^{2}^{–}^{5}^{)} defines the total energy of each subsystem X (\(X=I,IJ,IJK,\ldots\)) as \(E_{X}\), including the environmental electrostatic potential from other fragments, where \(I,J,K\) represent fragment indices and their combinations refer to the dimer and higher-order fragments. Once each monomer energy within the environmental potential is self-consistently obtained, the total energy of the entire system can be expressed in a many-body expansion manner with a correction energy of \(\Delta E_{X}\): \begin{equation} E^{\text{FMO}} = \sum_{I}E_{I}+\sum_{I>J}\Delta E_{IJ}+\sum_{I>J>K}\Delta E_{IJK}+\cdots. \end{equation} (1) In the pairwise approximation of FMO (FMO2), dimer energies are required in the correction formula \begin{equation} \Delta E_{IJ} = E_{IJ}-E_{I}-E_{J}, \end{equation} (2) which indicates the difference in energy from the monomers.
In the FMO-LCMO method, the total Hamiltonian matrix of the entire system is analogically constructed in the monomer MO representation^{7}^{)} \begin{equation} H^{\text{FMO}} = H^{(1)}+\Delta H^{(2)}+\Delta H^{(3)}+\cdots \end{equation} (3) by adding the ζ-th order corrections \(\Delta H^{(\zeta)}\) to the first-order Hamiltonian matrix \(H^{(1)}\). Each Hamiltonian matrix consists of its submatrices of each subsystem X having fragment MOs \(|\phi^{X}\rangle\) and associated energies \(\varepsilon^{X}\) with their indices of \(p, q, r, s\) as follows: \begin{align} H_{pq}^{Y[X]} &= \langle\phi_{p}^{I}|\hat{H}^{X}|\phi_{q}^{J}\rangle\\ &= \sum_{r}\sum_{s}\langle\phi_{p}^{I}|\phi_{r}^{X}\rangle H_{rs}^{X}\langle\phi_{s}^{X}|\phi_{q}^{J}\rangle\\ &= \sum_{r}\varepsilon_{r}^{X}\langle\phi_{p}^{I}|\phi_{r}^{X}\rangle\langle\phi_{r}^{X}|\phi_{q}^{J}\rangle, \end{align} (4) where Y represents a block of the Hamiltonian matrix based on the monomer MOs, taking only a diagonal block (\(Y = I\), \(I = J\)) or off-diagonal block (\(Y = IJ\), \(I\neq J\)). The superscript \(Y[X]\) denotes the contribution of the Hamiltonian components from a subsystem X to a block Y of the entire system. When Y is equal to X, \(H_{pq}^{Y[X]}\) is simply expressed as \begin{equation} H_{pq}^{Y[X]} = \varepsilon_{p}^{X}\delta_{pq}, \end{equation} (5) where \(\delta_{pq}\) is Kronecker's delta. The matrix element for the first-order Hamiltonian is then represented as \begin{equation} H_{pq}^{(1)} = \sum_{I}H_{pq}^{I[I]}, \end{equation} (6) and that for the second-order correction is also given as \begin{equation} \Delta H_{pq}^{(2)} = \begin{cases} \displaystyle\sum_{J \neq I}(H_{pq}^{I[IJ]}-H_{pq}^{I[I]}) &\text{(diagonal block)}\\ H_{pq}^{IJ[IJ]} &\text{(off-diagonal block)} \end{cases}, \end{equation} (7) as required for the pairwise approximation of FMO-LCMO (FMO2-LCMO). Figure 1 depicts the contribution of matrix elements above to the total Hamiltonian matrix. After the construction of the total Hamiltonian matrix, solving the generalized eigenvalue problem for the matrix yields delocalized LCMOs with one-electron orbital energies and the corresponding coefficients in the monomer MO representation. Because there exists arbitrariness in spanning the basis space of the monomer MOs, their HOMO/LUMO vicinities are uniformly selected as the bases in a conventional scheme. In the present research, as shown below, we confine ourselves to the consideration of FMO2-LCMO as a first step in the formulation.
Figure 1. (Color online) Construction of the total Hamiltonian matrix in the FMO2-LCMO method for a system containing N fragments.
Basis selection based on density matrix
Herein, we express the MOs of the fragment monomer and dimer as \(|\phi_{p}\rangle\) and \(|\psi_{t}\rangle\), respectively. A matrix element of the pseudo density matrix \(\mathbf{D}\) for each monomer, \(\rho_{pq}\), is defined as follows: \begin{equation} \rho_{pq} = \sum_{t}\langle\phi_{p}|\psi_{t}\rangle\langle\psi_{t}|\phi_{q}\rangle, \end{equation} (8) where the target range of t for dimers can be taken arbitrarily and may be near the frontier orbitals of the dimers. To emphasize the inter-relationship between the fragment monomer and dimer MOs, we hereafter call this \(\mathbf{D}\) the trans-fragment density matrix (TFDM).^{21}^{)} Upon diagonalizing \(\mathbf{D}\), we obtain the matrix \(\mathbf{L}\), which satisfies the following equation with the identity matrix \(\mathbf{I}\): \begin{equation} \mathbf{L}^{-1}\mathbf{DL}= \vec{\Lambda}\mathbf{I}, \end{equation} (9) where \(\mathbf{L}\) constitutes eigenvectors \(\vec{a}_{1},\vec{a}_{2},\ldots,\vec{a}_{j}, \ldots \) associated with eigenvalues \(\Lambda_{1},\Lambda_{2},\ldots,\Lambda_{j}, \ldots \), respectively. The eigenvalues correspond to the occupation numbers, and the associated eigenvectors consist of the coefficients for the MOs of each monomer. The most efficient way to compress the LCMO space based on the monomer MOs is to select the monomer MOs that contribute significantly to the eigenvalues and eigenvector components. For this selection, we introduce the following contribution index for the i-th monomer MO: \begin{equation} c_{i} := \sum_{j}^{N}\Lambda_{j}^{\ell} |a_{ij}|^{m} > c^{\text{min}}, \end{equation} (10) where N is the number of eigenvalues, \(\ell\) is the power of the eigenvalues \(\Lambda_{j}\), and m is the power of the absolute value of the eigenvector components \(a_{ij}\). The monomer MO is thus adopted as the basis when \(c_{i}\) is larger than the threshold \(c^{\text{min}}\). Once the essential monomer MOs have been extracted for each dimer, their aggregation is considered as the basis set for the LCMO. We call such an adaptive selection based on TFDM a DM scheme, and distinguish it from the conventional scheme based on a predetermined range of MO energy levels, which herein is called the EL scheme. More specifically, the DM scheme first targets the dimer MOs and then selects the monomer MO bases according to the fitness for the dimer MOs, whereas the EL scheme directly targets the monomer MOs used as the bases as is. In the present study of the DM scheme, we fix the parameter values as \(\ell = 1\) and \(m = 2\). For \(c^{\text{min}}\), we investigated its effect within the range of 0.01, 0.05, 0.1, 0.2, and 0.5.
We make some supplementary remarks on the present scheme. Let us consider a density operator in terms of the following wave functions \(|\Psi_{\nu}\rangle\): \begin{equation} \hat{W} = \sum_{\nu}w_{\nu}|\Psi_{\nu}\rangle\langle\Psi_{\nu}|, \end{equation} (11) with the coefficients \(w_{\nu}\). As a special case for this, we consider a projection operator as \begin{equation} \hat{P} = \sum_{\nu\in \textit{target}}|\Psi_{\nu}\rangle\langle\Psi_{\nu}|, \end{equation} (12) where \(w_{\nu} = 1\) when ν belongs to the target states and is \(w_{\nu} = 0\) otherwise. If we consider a density matrix with the basis set of fragment monomer MOs as \begin{equation} \rho_{pq} = \langle\phi_{p}|\hat{P}|\phi_{q}\rangle = \sum_{\nu\in \textit{target}}\langle\phi_{p}|\Psi_{\nu}\rangle\langle\Psi_{\nu}|\phi_{q}\rangle, \end{equation} (13) with “target” being the target (e.g., frontier) orbitals of the dimer MOs, it recovers the density matrix of Eq. (8). Furthermore, if ν in Eq. (12) forms a complete system, we find \begin{equation} 1 = \sum_{\nu}|\Psi_{\nu}\rangle\langle\Psi_{\nu}|, \end{equation} (14) and Eq. (13) becomes \begin{equation} \rho_{pq} = \delta_{pq}, \end{equation} (15) thus implying that the selection of the target states is essential. Whereas the earlier FMO-LCMO scheme corresponds to an approximation that the monomer MOs around the HOMO and LUMO are selected as the target states, the present scheme focuses on the inclusion of inter-fragment orbitals associated with the charge transfer with the exchange effects. This idea thus complies with the essence of the density matrix renormalization group (DMRG)^{22}^{–}^{24}^{)} concerning the inclusion of boundary effects upon the expansion of the system size.
In addition, the present renormalization procedure is intimately associated with the singular value decomposition (SVD),^{25}^{–}^{28}^{)} which is represented as follows: \begin{equation} \langle\phi_{p}|\psi_{r}\rangle = \sum_{\alpha}U_{p\alpha}\lambda_{\alpha}V_{r\alpha}, \end{equation} (16) where \(U_{p\alpha}\), \(V_{r\alpha}\), and \(\lambda_{\alpha}\) refer to the left-hand matrix, right-hand matrix, and eigenvalue of the SVD, respectively. If Eq. (16) is substituted into Eq. (8), we find \begin{equation} \rho_{pq} = \sum_{r}\sum_{\alpha}\sum_{\beta}U_{p\alpha}\lambda_{\alpha}V_{r\alpha}V_{r\beta}\lambda_{\beta}U_{q\beta}, \end{equation} (17) and hence \begin{equation} \mathop{\text{Tr}}\nolimits\rho_{pq} = \sum_{p,q,r}\sum_{\alpha,\beta}\delta_{pq}U_{p\alpha}U_{q\beta}V_{r\alpha}V_{r\beta}\lambda_{\alpha}\lambda_{\beta}. \end{equation} (18) If we assume the monomer MOs \(\{\phi_{p}\}\) form a complete system and use \begin{equation} \sum_{p}U_{p\alpha}U_{p\beta} = \delta_{\alpha\beta}, \end{equation} (19) we are led to \begin{equation} \mathop{\text{Tr}}\nolimits\rho_{pq} = \sum_{\alpha}\sum_{r}\lambda_{\alpha}^{2}V_{r\alpha}^{2}. \end{equation} (20) When the states r form a complete system, we find that \(\sum_{r}V_{r\alpha}^{2}=1\) and \(\lambda_{\alpha}^{2}\) refer to the eigenvalues of the density matrix \(\rho_{pq}\). The idempotence of the density matrix is retained when the basis set forms a complete system, but not when some truncations are introduced.
Computational details
In this study, we employ the Hartree–Fock (HF) approximation with the 6-31G^{*} basis set for the FMO2-LCMO analysis. Although more accurate descriptions based on the density functional theory (DFT) are feasible, we suppose that the prototypical HF-level analysis will provide essential insights into the development of the FMO-LCMO methodology. We used the developer version of ABINIT-MP^{5}^{,}^{29}^{)} for the following analysis.
In the following, we compare the results of the FMO-LCMO scheme obtained by the earlier EL-based method and the present DM-based method with those of the canonical orbitals obtained using a conventional MO method. When employing all monomer MOs as bases of FMO-LCMO, the scheme is called full FMO-LCMO, which indicates the ultimate performance of the FMO-LCMO method. The resultant orbital of the full FMO-LCMO should be distinguished from the canonical orbital. Regarding the choice of target orbitals (“window”), we consider three schemes: 1) full, 2) HOMO−n:LUMO+n, and 3) LOMO:LUMO+n. For the HOMO−n:LUMO+n scheme, we account for HOMO−n, …, HOMO−1, HOMO, LUMO, LUMO+1, …, LUMO+n. For the LOMO:LUMO+n scheme, we consider all occupied orbitals, where LOMO is the lowest occupied molecular orbital. To distinguish between localized monomer MOs and delocalized MOs such as LCMOs that linearly combine the localized MOs, or canonical MOs, we use the superscript “(d)” for the latter MOs, as MO\(^{\text{(d)}}\). The redundancy of the bases derived from the bond detached atom (BDA)^{5}^{)} is removed by the diagonalization of overlap integrals when we do not employ the minimal basis set.^{8}^{)}
3. Results and Discussion
In this section, we demonstrate the results of the DM-based LCMO scheme using stacked base pairs of deoxyribonucleic acid (DNA), pseudo-glycine pentamer, and water clusters as test systems in comparison with those of the conventional EL-based scheme.
Stacked base pairs of DNA
First, we consider the stacked Watson–Crick (WC) base pairs of DNA with four types of AA, AT, CG, and GA.^{30}^{)} Figure 2 shows the structure of CG stacked base pairs, where each base is a single fragment without BDA. We carried out the DM (HOMO−n:LUMO+n) calculations in comparison with the EL scheme. Figure 3 illustrates the comparison results concerning the distribution of selected monomer MOs in the case of the CG stacking, where the corresponding orbital energies of each monomer are shown in Fig. S1 in the Supplemental Material section.^{31}^{)} For the EL scheme, only the frontier MOs (HOMO:LUMO) of all monomers are uniformly selected as the bases in the case of \(n = 0\), and each monomer MO located just at the lower and upper energy positions of the above bases are added in the case of \(n = 1\), and further monomer MOs are added in the case of \(n = 2\); the total numbers of selected bases are 8, 16, and 24, respectively. For the DM scheme, the selected bases vary depending on the overlapped dimer orbitals and the threshold value of the contribution index. For the basis selected with the threshold of 0.5, other monomer MOs are added stepwise to the bases according to the reduction of the threshold of 0.2, 0.1, 0.05, and 0.01, and the total numbers of selected bases with the threshold of 0.1 are 3, 6, and 7 for the targeted dimer MO windows (\(n = 0\), 1, and 2), respectively. As seen in this comparison, the present DM scheme adaptively selects the monomer MO bases by fitting to the designated dimer orbitals with the energy width n and given threshold parameters, whereas the EL scheme homogeneously selects the monomer MO bases around HOMO/LUMO according to the designated value of n for the monomer MO energy levels. Thus, in the present DM scheme, the essential monomer MOs are efficiently selected according to the threshold parameters such that the number of monomer MO bases is adaptively changed. We then compare the resultant FMO-LCMO orbitals with the canonical orbitals in Table I. In contrast to the full 610 bases, both the EL and DM schemes well reproduce the orbital energies around HOMO\(^{\text{(d)}}\)/LUMO\(^{\text{(d)}}\) with only tens of monomer MO bases. However, unnatural orbitals appear at LUMO\(^{\text{(d)}}\) and LUMO+1\(^{\text{(d)}}\) in the case of \(n = 2\) of the EL scheme, leading to a reduction in the HOMO\(^{\text{(d)}}\)/LUMO\(^{\text{(d)}}\) gap. In principle, such unnatural orbitals must be shifted to higher energy levels owing to their orthogonality with lower-energy-level occupied orbitals that should canonically exist. In the FMO-LCMO scheme, which drastically limits the basis space of the monomers, the occupied orbitals in depth are often omitted, resulting in such energy level discrepancies. This difficulty can be attributed to the fact that the required orthogonality to the occupied orbitals is not attained in the extended basis sets beyond the minimal set.^{8}^{)}
Figure 3. (Color online) Monomer MO distribution as selected bases in CG stacked base pairs, determined by EL/DM-based FMO-LCMO schemes targeting the MO energy window of HOMO−n:LUMO+n (\(n=0,1,2\)). The basis labels of MO energy levels below the HOMO and above the LUMO are abbreviated as H−k and L+k (\(k=1,2,\ldots\)), respectively. For EL, only the frontier MOs (HOMO:LUMO) of all monomers are uniformly selected as the bases in the case of \(n = 0\) (the bases are depicted in a color labeled as \(n = 0\)), one monomer MOs located just at the lower and upper energy positions of the above bases are added in the case of \(n = 1\) (the increment bases for \(n = 0\) are depicted in a color labeled as \(n = 1\)), and further ones are similarly added in the case of \(n = 2\) (the increment bases for \(n = 1\) are depicted in a color labeled as \(n = 2\)). For DM, the selected bases vary depending on the overlapped dimer orbitals and the threshold value of the contribution index. When some bases selected with the threshold of 0.5 are depicted in color and labeled as 0.5, other monomer MOs are added stepwise to the bases according to the reduction of the threshold of 0.2, 0.1, 0.05, and 0.01 (the increment bases for 0.5, 0.2, 0.1, and 0.05 are depicted in colors labeled as 0.2, 0.1, 0.05, and 0.01, respectively).
Table I. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of HOMO−n:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for CG stacked base pairs.
To cope with this orthogonality problem, we further conducted a DM (LOMO:LUMO+n) analysis in comparison with the corresponding EL calculation, where all occupied orbitals down to the LOMO and the unoccupied orbitals near the LUMO are considered. The distribution of the monomer MOs selected as bases is shown in Fig. 4, and the corresponding orbital energies of the monomer MOs are shown in Fig. S1 in the Supplemental Material section,^{31}^{)} where all occupied orbitals are employed as bases in both the EL and DM schemes. This comparison shows that the unoccupied orbitals of the monomers selected in the DM scheme are distributed discontinuously and broadly according to the criteria of the target width and contribution threshold. Figure 5 then illustrates the dependence of LCMO energy levels near the LUMO\(^{\text{(d)}}\) on the threshold of the contribution index in the DM scheme. Thus, we observe that the energy levels of LUMO\(^{\text{(d)}}\) and LUMO+1\(^{\text{(d)}}\) converge well within the range of 0.01–0.5 of the threshold parameter. In this case, for example, the LUMO of fragment 4 does not contribute to the LUMO\(^{\text{(d)}}\) of LCMO in the case of an \(n = 0\) target fitting. However, a target value of \(n = 1\) or higher and a threshold value of 0.1 or less are required to reproduce the LUMO+2\(^{\text{(d)}}\) orbital energy of LCMO. Table II shows a comparison of the LCMO energy levels by EL and DM methods with those of the canonical orbitals in the case of targeting LOMO:LUMO+n for CG stacked base pairs. Thus, we find the resolution of the orthogonality difficulty, leading to an accurate reproduction of the canonical orbitals in both the EL and DM schemes, whereas the required numbers of monomer bases are substantially increased up to 140–150.
Figure 4. (Color online) Monomer MO distribution as selected bases in CG stacked base pairs, determined through EL/DM-based FMO-LCMO schemes targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)). The color coding follows the legend in Fig. 3. The basis labels of monomer MO energy levels below the HOMO and above the LUMO are abbreviated as H−k and L+k (\(k = 1, 2,\ldots\)), respectively.
Figure 5. (Color online) Dependence of delocalized MO energy levels near LUMO on the threshold of contribution index in DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) for CG stacked base pairs.
Table II. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of LOMO:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for CG stacked base pairs.
Furthermore, we carried out analogous analyses for AA-, AT-, and GA-stacked WC base pairs, the structures of which are depicted in Figs. S2, S3, and S4, respectively, in the Supplemental Material section.^{31}^{)} The dependence of the FMO-LCMO energy levels near the LUMO\(^{\text{(d)}}\) on the threshold parameter \(c^{\text{min}}\) in the DM-based LOMO:LUMO+n target scheme are illustrated in Figs. S5–S7 in the Supplemental Material section^{31}^{)} in comparison with the full FMO-LCMO and canonical MO results. Tables S1–S6 in the Supplemental Material section^{31}^{)} show comparisons of the orbital energy levels between the EM/DM based FMO-LCMO results and the canonical MO results in the cases of HOMO−n:LUMO+n and LOMO:LUMO+n target schemes. Taken together, the results for AA, AT, and GA stacked pairs are similar to those for the CG case above but indicate that the orthogonality problem may occur because of the increase in the number of MO bases even in the DM and EL schemes. That is, whereas the HOMO−n:LUMO+n target scheme seems to work well in the DM case for the stacked base pairs of DNA, we should focus on the appearance of unnatural orbitals because of the orthogonality difficulty when the number of monomer MO bases increases.
Pseudo-glycine pentamer
Next, we consider a pseudo-glycine pentamer molecule with the BDA, where the N-terminal glycine residue is modified from NH_{2}CH_{2}CO– to HCONHCH_{2}CO–, and the C-terminal glycine residue is modified from –NHCH_{2}COOH to –NHCH_{3}.^{7}^{)} Figure 6 shows its structure and fragmentation. In this system, it is remarkable that the orthogonality difficulty owing to a reduction of occupied orbitals in the target vividly appears, as shown in Table III. In other words, when we employ the HOMO−n:LUMO+n target scheme, large deviations in the LCMO orbitals from the canonical orbitals are observed on the LUMO\(^{\text{(d)}}\) side in both the EL and DM schemes. It should also be noted that there are small errors in the orbital energies at the HOMO\(^{\text{(d)}}\) (−0.012 a.u.) and LUMO\(^{\text{(d)}}\) (−0.007 a.u.) levels, even in the full LCMO scheme.
Figure 6. (Color online) Structure and fragmentation of pseudo-glycine pentamer. The glycine residues at the N- and C-termini are modified as HCONHCH_{2}CO– and –NHCH_{3}, respectively, and thus called “pseudo” glycine residues. This peptide is divided into five fragments at the sites of each Cα atom. Residue codes are written with the fragment number in parentheses.
Table III. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of HOMO−n:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for pseudo-glycine pentamer.
Hence, we need the LOMO:LUMO+n target scheme to reproduce the canonical orbitals. Figure 7 shows the distribution of the monomer MOs selected as bases, and the corresponding orbital energies are shown in Fig. S8 in the Supplemental Material section.^{31}^{)} The distribution is nontrivial and discontinuous, which is analogous to the DNA base pair case. For example, LUMO+36 in monomer fragment 5 selected for \(c^{\text{min}} = 0.05\) provides a prototypical case. Figure 8 then illustrates the distribution of the LCMO energy levels near the HOMO\(^{\text{(d)}}\)/LUMO\(^{\text{(d)}}\) in the DM scheme. Although the occupied orbitals are almost identical in all cases, the distribution on the unoccupied side varies considerably according to the conditions of the target and threshold values, but with the converged LUMO\(^{\text{(d)}}\) level below \(c^{\text{min}} = 0.1\). Figure 9 illustrates the dependence of the LCMO energy levels near LUMO\(^{\text{(d)}}\) on the threshold of the contribution index in the DM-based LOMO:LUMO+n target schemes, which shows the relevance of \(c^{\text{min}}\leq 0.1\) in comparison with the energies of the canonical MOs. By considering all occupied MOs of the monomer bases, the FMO-LCMO orbitals near HOMO\(^{\text{(d)}}\)/LUMO\(^{\text{(d)}}\) are drastically improved in the DM scheme. Figure 10 shows the dependence of the delocalized MO shapes near LUMO\(^{\text{(d)}}\) on the threshold of the contribution index in the DM based LOMO:LUMO target schemes, and Fig. 11 shows those in the DM-based LOMO:LUMO+n target schemes with \(c^{\text{min}} = 0.1\). The delocalized MO shapes with \(c^{\text{min}} = 0.1\) reproduce those of the canonical shapes well, which also indicates the relevance of the threshold parameter. The delocalized MO energy levels also appear to be recovered through the EL scheme by including all occupied orbitals into the bases (see Table IV). However, the actual shapes of the delocalized orbitals, particularly LUMO\(^{\text{(d)}}\) and its uppers, are different from those in the DM scheme, indicating that the canonical orbitals are not well reproduced, as shown in Fig. 12. This shows, as suggested by the distribution of the monomer bases in the DM scheme in Fig. 7, that the monomer bases necessary for the description of the orbitals near LUMO\(^{\text{(d)}}\) are insufficiently considered. Thus, the EL scheme cannot provide a relevant solution to how many and in what manner high-energy orbitals of monomer bases should be considered, which can be replicated in the DM scheme.
Figure 7. (Color online) Monomer MO distribution as selected bases in pseudo-glycine pentamer, determined through DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)). The color coding follows the legend in Fig. 3. The basis labels of monomer MO energy levels below the HOMO and above the LUMO are abbreviated as H−k and L\( + k\) (\(k = 1, 2,\ldots\)), respectively.
Figure 8. (Color online) Delocalized MO energy levels determined through DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) for pseudo-glycine pentamer.
Figure 9. (Color online) Dependence of delocalized MO energy levels near LUMO on the threshold of contribution index in DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) for pseudo-glycine pentamer.
Figure 10. (Color online) Delocalized MOs near frontier MOs generated by DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO while changing threshold parameter \(c^{\text{min}}\) for pseudo-glycine pentamer. The overlap integrals with canonical MOs are shown below each molecular orbital.
Figure 11. (Color online) Delocalized MOs near frontier MOs generated by DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) with \(c^{\text{min}}=0.1\) for pseudo-glycine pentamer. The overlap integrals with canonical MOs are shown below each molecular orbital.
Figure 12. (Color online) Delocalized MOs near frontier MOs generated by EL-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) for pseudo-glycine pentamer. The overlap integrals with canonical MOs are shown below each molecular orbital.
Table IV. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of LOMO:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for pseudo-glycine pentamer.
To advance the discussion, we now perform a model calculation in which the target orbitals are confined to those near the LUMO, which is referred to as the None:LUMO+n target scheme. The distribution of the monomer MOs selected as bases is shown in Fig. 13, and the energy levels of corresponding monomer MOs are found in Fig. S8 in the Supplemental Material section.^{31}^{)} The monomer MO bases selected in the DM scheme are then only unoccupied orbitals. However, as shown in Fig. 14, some of the LCMO orbitals obtained with these monomer MO bases appear in the occupied orbital region, thus reflecting the orthogonality problem. Through a comparison with Fig. 8, we can see that the orbitals well above LUMO\(^{\text{(d)}}\) are brought into the occupied orbital region owing to a lack of orthogonality in the occupied orbitals.
Figure 13. (Color online) Monomer MO distribution as selected bases in pseudo-glycine pentamer, determined through DM-based FMO-LCMO scheme targeting the MO energy window of None:LUMO+n (\(n = 0, 1, 2\)). The color coding follows the legend in Fig. 3. The basis labels of MO energy levels below the HOMO and above the LUMO are abbreviated as H−k and L+k (\(k = 1, 2,\ldots\)), respectively.
Figure 14. (Color online) Delocalized MO energy levels determined by DM based FMO-LCMO scheme targeting the MO energy window of None:LUMO+n (\(n = 0, 1, 2\)) for pseudo-glycine pentamer.
Water clusters
Finally, we consider water clusters of dimer, trimer, and tetramer,^{32}^{)} as shown in Fig. 15. Note that the water molecules as fragments are situated in a heterogeneous environment in the dimer and trimer cases.
Figure 15. (Color online) Structures of water dimer, trimer, and tetramer. They are in the \(C_{\text{s}}\) symmetry with a hydrogen bond, in the \(C_{1}\) symmetry with three hydrogen bonds, and in the \(S_{4}\) symmetry with four hydrogen bonds, respectively. Water molecules are written with the fragment number in parentheses.
Figures 16–18 illustrate the calculated results for LCMO energy levels near LUMO\(^{\text{(d)}}\) as functions of the threshold of the contribution index (\(c^{\text{min}}\)) in the DM scheme for LOMO:LUMO+n targets. It is apparent that different behaviors were observed between the water dimer and trimer/tetramer cases, and the \(c^{\text{min}}\) dependence in the latter is not smooth but somewhat discrete. This suggests that there are MO bases that are not regarded as important in the fragment dimer level but play a significant role in the LCMO orbitals in the cases of trimers and tetramers. This fact may also indicate the necessity of a description beyond the FMO2 approximation for water trimer and tetramer systems.^{18}^{)}
Figure 16. (Color online) Dependence of delocalized MO energy levels near LUMO on the threshold of contribution index in DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) for water dimer.
Figure 17. (Color online) Dependence of delocalized MO energy levels near LUMO on the threshold of contribution index in DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) for water trimer.
Figure 18. (Color online) Dependence of delocalized MO energy levels near LUMO on the threshold of contribution index in DM-based FMO-LCMO scheme targeting the MO energy window of LOMO:LUMO+n (\(n = 0, 1, 2\)) for water tetramer.
We carried out the HOMO−n:LUMO+n and LOMO:LUMO+n target calculations in both the EL- and DM-based schemes, the results of which are compiled in Tables V –X. Thus, it can be observed that the reproduction of the canonical results is unsatisfactory near LUMO\(^{\text{(d)}}\) when we employ the HOMO−n:LUMO+n target. However, we can achieve a satisfactory reproduction if we employ the LOMO:LUMO+n target to accommodate LCMO orbitals that should be appropriately described. Because there still remain certain ambiguities in the present study owing to the neglect of higher-order effects beyond the fragment dimer approximation, we may need further investigations concerning this issue.
Table V. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of HOMO−n:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for water dimer.
Table VI. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of LOMO:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for water dimer.
Table VII. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of HOMO−n:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for water trimer.
Table VIII. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of LOMO:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for water trimer.
Table IX. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of HOMO−n:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for water tetramer.
Table X. Comparison of delocalized MO energy levels in a.u. by EL/DM-based FMO-LCMO schemes targeting the MO energy window of LOMO:LUMO+n (\(n=0,1,2\)) to those by full FMO-LCMO scheme using all monomer MO bases, and canonical MO scheme, for water tetramer.
4. Conclusion
Whereas the present basic study has focused on relatively small molecular systems, it is becoming increasingly difficult to apply diagonalization with monomer-based orbitals when one constructs an FMO-LCMO for large-scale molecular systems, in which the full inclusion of the monomer MOs will be intractable. One of our strategies is to confine the basis space to that near the HOMO/LUMO, although we need to fully account for the occupied basis orbitals when the orthogonality with respect to the delocalized occupied orbitals is important. Under these circumstances, it is essential to efficiently select the basis space necessary for the relevant description of an LCMO. As observed above for the distribution of monomer MO bases selected in the DM scheme, the basis selection in the EL scheme lacks a rational justification for reproducing the orbitals near HOMO\(^{\text{(d)}}\)/LUMO\(^{\text{(d)}}\); That is, the orbital energy range necessary for a relevant description varies according to the specific molecular systems. By contrast, the DM scheme can provide a more robust methodology for selecting the monomer MO bases because it relies on a universal threshold of contribution index, irrespective of molecular systems. Thus, we found an advantage of the DM scheme such that it enables a rational selection of monomer MO bases to accurately construct the canonical orbitals for the entire system. Through the present analysis, we found a threshold value of approximately \(c^{\text{min}} = 0.1\) to be a reasonable choice in light of the compromise between the accuracy and cost of diagonalization.
Thus, the present density-matrix-based FMO-LCMO method is superior to the energy-based method in terms of the efficiency of data compression of monomer MO bases adopted to construct the delocalized MOs for the entire system. Our criterion for selecting the localized MOs of fragment monomers can provide a systematic way to control the cost performance of an approximation.
Acknowledgements
We thank Professor Takatoshi Fujita for the useful discussions. S.T. would like to acknowledge the Grants-in-Aid for Scientific Research (Nos. 17H06353 and 18K03825) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT).
References
1 Fragmentation: Toward Accurate Calculations on Complex Molecular Systems, ed. M. S. Gordon (Wiley, Chichester, U.K., 2017). Crossref, Google Scholar
2 K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayasi, Chem. Phys. Lett. 313, 701 (1999). 10.1016/S0009-2614(99)00874-X Crossref, Google Scholar
3 T. Nakano, T. Kaminuma, T. Sato, Y. Akiyama, M. Uebayasi, and K. Kitaura, Chem. Phys. Lett. 318, 614 (2000). 10.1016/S0009-2614(00)00070-1 Crossref, Google Scholar
4 D. G. Fedorov, T. Nagata, and K. Kitaura, Phys. Chem. Chem. Phys. 14, 7562 (2012). 10.1039/c2cp23784a Crossref, Google Scholar
5 S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, and K. Fukuzawa, Phys. Chem. Chem. Phys. 16, 10310 (2014). 10.1039/C4CP00316K Crossref, Google Scholar
6 Y. Inadomi, T. Nakano, K. Kitaura, and U. Nagashima, Chem. Phys. Lett. 364, 139 (2002). 10.1016/S0009-2614(02)01291-5 Crossref, Google Scholar
7 S. Tsuneyuki, T. Kobori, K. Akagi, K. Sodeyama, K. Terakura, and H. Fukuyama, Chem. Phys. Lett. 476, 104 (2009). 10.1016/j.cplett.2009.05.069 Crossref, Google Scholar
8 T. Kobori, K. Sodeyama, T. Otsuka, Y. Tateyama, and S. Tsuneyuki, J. Chem. Phys. 139, 094113 (2013). 10.1063/1.4818599 Crossref, Google Scholar
9 H. Nishioka and K. Ando, J. Chem. Phys. 134, 204109 (2011). 10.1063/1.3594100 Crossref, Google Scholar
10 H. Kitoh-Nishioka and K. Ando, Chem. Phys. Lett. 621, 96 (2015). 10.1016/j.cplett.2014.12.057 Crossref, Google Scholar
11 H. Kitoh-Nishioka, K. Welke, Y. Nishimoto, D. G. Fedorov, and S. Irle, J. Phys. Chem. C 121, 17712 (2017). 10.1021/acs.jpcc.7b05779 Crossref, Google Scholar
12 R. Sato, H. Kitoh-Nishioka, T. Yanai, and Y. Shigeta, Chem. Lett. 46, 873 (2017). 10.1246/cl.170161 Crossref, Google Scholar
13 H. Kitoh-Nishioka, Y. Shigeta, and K. Ando, J. Chem. Phys. 153, 104104 (2020). 10.1063/5.0018423 Crossref, Google Scholar
14 T. Fujita and Y. Noguchi, Phys. Rev. B 98, 205140 (2018). 10.1103/PhysRevB.98.205140 Crossref, Google Scholar
15 Y. Suzuki, K. Ebina, and S. Tanaka, Chem. Phys. 474, 18 (2016). 10.1016/j.chemphys.2016.05.001 Crossref, Google Scholar
16 Y. Suzuki, H. Watanabe, Y. Okiyama, K. Ebina, and S. Tanaka, Chem. Phys. 539, 110903 (2020). 10.1016/j.chemphys.2020.110903 Crossref, Google Scholar
17 D. G. Fedorov and K. Kitaura, Chem. Phys. Lett. 433, 182 (2006). 10.1016/j.cplett.2006.10.052 Crossref, Google Scholar
18 T. Fujita, K. Fukuzawa, Y. Mochizuki, T. Nakano, and S. Tanaka, Chem. Phys. Lett. 478, 295 (2009). 10.1016/j.cplett.2009.07.060 Crossref, Google Scholar
19 T. Nakano, Y. Mochizuki, K. Yamashita, C. Watanabe, K. Fukuzawa, K. Segawa, Y. Okiyama, T. Tsukamoto, and S. Tanaka, Chem. Phys. Lett. 523, 128 (2012). 10.1016/j.cplett.2011.12.004 Crossref, Google Scholar
20 D. G. Fedorov and K. Kitaura, J. Chem. Phys. 147, 104106 (2017). 10.1063/1.5001018 Crossref, Google Scholar
(21)
The idempotency is not necessarily required for this TFDM. Google Scholar
22 S. R. White, Phys. Rev. Lett. 69, 2863 (1992). 10.1103/PhysRevLett.69.2863 Crossref, Google Scholar
23 S. R. White, Phys. Rev. B 48, 10345 (1993). 10.1103/PhysRevB.48.10345 Crossref, Google Scholar
24 T. Xiang, Phys. Rev. B 53, R10445 (1996). 10.1103/PhysRevB.53.R10445 Crossref, Google Scholar
25 V. Klema and A. Laub, IEEE Trans. Autom. Control 25, 164 (1980). 10.1109/TAC.1980.1102314 Crossref, Google Scholar
26 L. De Lathauwer, B. De Moor, and J. Vandewalle, SIAM J. Matrix Anal. Appl. 21, 1253 (2000). 10.1137/S0895479896305696 Crossref, Google Scholar
27 F. Kleibergen and R. Paap, J. Econ. 133, 97 (2006). 10.1016/j.jeconom.2005.02.011 Crossref, Google Scholar
28 M. Brand, Linear Algebra Appl. 415, 20 (2006). 10.1016/j.laa.2005.07.021 Crossref, Google Scholar
29 T. Nakano, Y. Mochizuki, K. Fukuzawa, S. Amari, and S. Tanaka, Developments and Applications of ABINIT-MP Software Based on the Fragment Molecular Orbital Method, in Modern Methods for Theoretical Physical Chemistry of Biopolymers, ed. E. B. Starikov, J. P. Lewis, and S. Tanaka (Elsevier, Amsterdam, 2006) Chap. 2, p. 39. 10.1016/B978-044452220-7/50066-6 Crossref, Google Scholar
30 K. Fukuzawa, C. Watanabe, I. Kurisaki, N. Taguchi, Y. Mochizuki, T. Nakano, S. Tanaka, and Y. Komeiji, Comput. Theor. Chem. 1034, 7 (2014). 10.1016/j.comptc.2014.02.002 Crossref, Google Scholar
31 (Supplemental Material) The monomer MO energy levels of CG stacked base pair of DNA; the molecular structures and FMO-LCMO results of AA, AT, GA stacked base pairs of DNA; the monomer MO energy levels of pseudo-glycine pentamer are provided online. Google Scholar
32 S. Maheshwary, N. Patel, N. Sathyamurthy, A. D. Kulkarni, and S. R. Gadre, J. Phys. Chem. A 105, 10525 (2001). 10.1021/jp013141b Crossref, Google Scholar