Evaluating First-Order Molecular Properties of Delocalized Ionic or Excited States in Molecular Aggregates by Renormalized Excitonic Method

Ⅰ.   INTRODUCTION

Theoretical understanding of various excited state processes in the condensed phases, ranging from optoelectronic conversion in material science to light harvesting in biological photosynthetic system, usually depends on an accurate description of their electronic structure features. However it is far beyond the capability of conventional ab initio excited state quantum chemistry methods such as configuration interaction singles (CIS) and complete active space self-consistent field (CASSCF) as well as time-dependent density functional theory (TDDFT), due to their too huge computational costs for the multi-chromophore systems [1-4]. Meanwhile, adopting a simple concept of Frenkel exciton (FE, a form of excitation wave) [5, 6] to approximately describe the essential features of tightly bound states of electron-hole pairs as a quasiparticle, model Hamiltonians utilizing FE model and/or its variants to incorporate charge-transfer (CT) excitons or interactions with environments can provide an economic and powerful tool to characterize the collective excitations in crystals and aggregates at a qualitative or semi-quantitative level, beneficial for the microscopic interpretation of electronic spectroscopy and ultrafast dynamics experiments in multi-chromophore systems [7-12]. For example, Sarovar, et al. presented a first rigorous quantification of quantum entanglement in excitation energy transport for a biological Fenna-Matthews-Olson protein system by virtue of modelling excitation dynamics within the FE model [13]. Recently, we also reported an interesting super-radiance behavior in a class of monolayer $J$-aggregates for two-dimensional (2D) molecular crystals through simulating the aggregate emission spectrum by a FE-vibration interaction model with experimental verification [14].

Despite the great success of its wide applications in excited state problems, empirical FE model still suffers from a few drawbacks. First of all, the definition of model basis and the construction of the model are not automatic. The choice of its parameter values sometimes relies heavily on experimental fitting or experience setting, greatly restricting its further extension to new systems. Secondly, the accurate derivation of parameter values from quantum chemical calculations is also challenging [15-18]. Not only long-range electrostatic couplings but also short-range exchange ones are necessary to be adequately considered for evaluating inter-chromophore coupling terms. However, it is nontrivial but especially crucial for small molecule aggregates, in which the inter-chromophore distance is short [19]. On the other hand, more than one excited state at each chromophore with degenerate or slightly higher excitation energies may also get involved in the photo-physics of the whole system, but FE usually considers only one local excitation per site as the model basis which might be not optimal. For example, it is revealed that the energetic position of the second lowest excited state in high-symmetry molecules such as fullerene and moderately to largely sized $\pi$-conjugated macromolecules can remarkably modulate the exciton dissociation rate in organic semiconductors by up to more than two orders of magnitude [20, 21]. How to account for the influence by these energetically higher excited states at each site with an economic price is still an open question.

Describing the inter-chromophore interactions by employing the Bloch's effective Hamiltonian theory [22, 23] may provide an efficient way to improve the accuracy of evaluating couplings and account for the effect by other excited states. By projecting the exact solution of the chromophore dimer at a specified level (e.g. CIS or TDDFT) onto the model space built by product states of a small number of local monomer states (obtained from the monomer solution at the same level), the effective Hamiltonian based on this small model space is able to exactly reproduce the dimer's low-lying excitation energies with a price of losing the accuracy of wavefunctions and accordingly cover the inter-chromophore exchange interactions and partly the effect by other monomer excitations. Applying this idea to larger systems, contractor renormalization group (CORE) [24] and renormalized excitonic method (REM) [25] were proposed by assuming the low-lying excited states of the whole system will be a linear combination of various multiple or only single monomer excitations. Since 2012, we [26-29] implemented REM calculations at various ab initio levels of full configuration interaction (FCI), CIS, symmetry adapted cluster configuration interaction (SAC-CI) and TDDFT, and the tests of hydrogen chains, polyene, polysilenes and water chains, benzene aggregates as well as general aqueous systems with polar and nonpolar solutes and tetracene/PTCDA interface indicated that the ab initio REM method can give good descriptions of excitation energies for lowest electronic excitations and similarly ionic potentials for delocalized ionic states in large systems with economic computational costs. Nevertheless, evaluating wavefunction properties is important for investigating the spectrum intensity, population analysis as well as geometry optimization, but is completely absent in those works, greatly hindering the popularity of REM's applications.

In this work, we further extend the REM calculations for the evaluations of first-order molecular properties (e.g. charge population and transition dipole moment) of delocalized ionic or excited states in molecular aggregates, through generalizing the effective Hamiltonian theory to effective operator representation. Results from the test calculations for four different kinds of one dimensional (1D) molecular aggregates suggest that our new scheme provides an efficient computational tool for describing both the energetic and wavefunction properties of low-lying excited states in large-sized complicated systems.

Ⅱ.   METHODOLOGY AND COMPUTATIONAL DETAILS

A   $\mathit{\boldsymbol{Ab}}$ $\mathit{\boldsymbol{initio}}$ REM scheme

In our previous works [26-29], we gave detailed descriptions for REM at different ab initio levels of FCI, CIS, SAC-CI, and TDDFT. Here, we just give a brief introduction to the basic idea of REM. Ab initio REM is a fragment-based excited state quantum chemical method, by virtue of approximating the excited states of the whole system as only the linear combination of various single chromophore excitations. In REM, the whole system can be divided into $N$ blocks (i.e. chromophores), as illustrated in FIG. 1. Here $I$, $J$, $K$, and $L$ denote for block monomers, and additionally, the adjacent monomer pairs form the dimers.

Figure  1.  The partition of the whole system into block monomers and dimers.

DownLoad: Full-Size Img PowerPoint

Assuming one has an orthonormal set of molecular orbitals (MOs), which are localized in each block monomer, standard quantum chemical calculations of each monomer with this MO set at the levels of FCI, CIS, SAC-CI or TDDFT can obtain a few low-lying eigenstates for each monomer. For building the excitonic model, usually only a small number of these local eigenstates will be kept. For the simplest case, one can preserve only two states per block (the neutral ground state $\left| {\psi _I^0} \right\rangle$ and one excited state, or one ionic state for ionic systems $\left| {\psi _I^1} \right\rangle$ with ${\hat H_I}\left| {\psi _I^n} \right\rangle$ = $E_I^n\left| {\psi _I^n} \right\rangle$). Please note that the states on different blocks are orthogonal to each other because we use an orthonormal MO set. Then the direct product states of these local states (allowing only single block excitation/ionization) constitute a model space

and wavefunction will be the linear combination of bases in this model space

To obtain the coefficients $\{c_I\}$, one needs to build the effective Hamiltonian ($\hat H^{\text{eff}}$) matrix with model space bases. If only one-body and two-body inter-block interactions are considered, operator could be written as

where $\hat H_{I, J}^{\text{eff}}$ is effective interaction between block $I$ and $J$,

$\hat H_{IJ}^{{\text{eff}}}$ matrix can be evaluated by using Bloch effective Hamiltonian theory [22, 23].

Let us consider we calculate the first few excited states for the block dimer $IJ$ at the same quantum chemical level as done for monomers, i.e.

Considering the simplest case, $\left| {\psi _{IJ}^1} \right\rangle$ and $\left| {\psi _{IJ}^2} \right\rangle$ constitute a target space for dimer $IJ$. At the same time, a tiny model space for the block dimer system can also be constructed by using two product bases $\left| {\psi _I^1\psi _J^0} \right\rangle$ and $\left| {\psi _I^0\psi _J^1} \right\rangle$. The project operator for this dimer model space could be expressed as

The effective Hamiltonian of the dimer $\hat H_{IJ}^{{\text{eff}}}$ is expected to satisfy the following equations

Then, defining the basis transformation matrix between the iso-dimensional target space and model space for the dimer

The matrix representation of Eq.(6) can be rewritten as

where

Therefore, the matrix representation of $\hat H_{IJ}^{\text{eff}}$ within the model space shall be obtained by

However, matrix $\tilde{\mathit{\boldsymbol{U}}}$ is usually not orthonormal because both the target space and model space are not complete. To avoid $\mathit{\boldsymbol{H}}_{IJ}^{{\text{eff}}}$ becoming non-Hermitian, an extra orthonormalization process is necessary to be performed on matrix $\tilde{\mathit{\boldsymbol{U}}}$ before calculating $\mathit{\boldsymbol{H}}_{IJ}^{{\text{eff}}}$ using Eq.(10).

After obtaining all one-body and two-body effective Hamiltonians, we could calculate the effective Hamiltonian ${\hat H^{{\text{eff}}}}$ of the whole system using Eq.(2). Then the energy and the wavefunction of low-lying excited states ($\left| {\psi _{{\text{REM}}}^m} \right\rangle$) could be easily obtained by a standard diagonalization of a small $\mathit{\boldsymbol{H}}^{{\text{eff}}}$ matrix ($N$$\times$$N$). The excitation energy can also be computed by taking the ground state energy $E_0$ as the reference, where $E_0$ is

And the ground state wavefunction could be approximated by

Now, let us explain what kind of orthonormal MO set we use in our above REM calculations. In our previous works [26-29], we employed orthogonal localized molecular orbitals (OLMOs) and block canonical molecular orbitals (BCMOs), which function well in conjunction with REM framework. However, both of these two MO sets suffer from a few inconveniences, i.e. the expensive costs of a Hartree-Fock (HF) calculation of the whole system and the difficulty in localizing a large number of unoccupied MOs for OLMOs, and the numerical complexity in treating the nonorthogonality in BCMOs. In this work, we use symmetric orthogonalized block canonical molecular orbitals (SYM-BCMOs), which are generated by a simple symmetric orthogonalization treatment of BCMOs and maintain the features of locality and canonicality as well as orthogonality. The BCMOs of the model system is an assembly of the canonical MOs (CMOs) of HF calculations for different block monomers, i.e.

where $\left| {\phi _m^I} \right\rangle$ is the $m$th canonical MO of block monomer $I$. Then, the overlap between different BCMOs can be evaluated by

where $c_\mu^m (I)$ and $c_\nu^n (J)$ are elements of atomic orbital (AO) coefficient matrix $\mathit{\boldsymbol{C}}_0$ of BCMOs by assembling CMOs of each block monomer, denoting the coefficient of the $\mu$th ($\nu$th) AO in the $m$th ($n$th) CMO of block monomer $I (J)$. $\chi_{\mu\nu}$ is the overlap between the two AOs. Having the overlap matrix $\mathit{\boldsymbol{S}}$ for SYM-BCMOs, the AO coefficient matrix $\mathit{\boldsymbol{C}}$ of SYM-BCMOs can be easily gained by a standard symmetric orthogonalization process, i.e.

B   Evaluation of first-order properties within REM framework

In this work, we further use REM framework to evaluate molecular properties based on first-order reduced density matrix (1-RDM), whose operator can be described as

in which $p$ and $q$ denote for SYM-BCMO indices. However, the REM wavefunction is not accurate because the monomer state bases are incomplete and not optimal for inter-block interactions. Therefore, instead of applying 1-RDM operator to REM wavefunction ${\left| {\psi _{{\text{REM}}}^n} \right\rangle }$ directly, we generalize Bloch's effective Hamiltonian theory to other operators including $\hat D_{pq}^{{\text{eff}}}$. The total $\hat D_{pq}^{{\text{eff}}}$ is similarly approximated as a summation of one-body and two-body terms as

Likewise, suppose one gets 1-RDMs for the two eigenstates of dimer $IJ$ ($\left| {\psi _{IJ}^1} \right\rangle$ and $\left| {\psi _{IJ}^2} \right\rangle$) after a standard quantum chemical calculation, i.e.

If we define

one can easily calculate the $D_{pq}^{{\text{eff}}}\left( {IJ} \right)$ matrix in model space basis by

where matrix $\tilde{\mathit{\boldsymbol{U}}}$ is defined by Eq.(7).

Taking ($p$, $q$) pair as an example, after obtaining all one-body and two-body effective 1-RDMs ($\mathit{\boldsymbol{D}}_{pq}^{{\text{eff}}}\left( I \right)$, $\mathit{\boldsymbol{D}}_{pq}^{{\text{eff}}}\left( {IJ} \right)$), we can build ($N$$\times$$N$)-dimensional matrix $\mathit{\boldsymbol{D}}_{pq}^{{\text{eff}}}$ using Eq.(16). Then the element for the final 1-RDM of the $n$th excited state of the whole system can be calculated by

Finally, by assembling all ($p$, $q$) pairs we can construct the (${N_{{\text{basis}}}}$$\times$${N_{{\text{basis}}}}$)-dimensional total 1-RDM of the whole system in SYM-BCMO basis (${\mathit{\boldsymbol{D}}^{{\text{eff}}}}$(SYM-BCMO)), where ${N_{{\text{basis}}}}$ is the number of total AO basis in the whole system. For the purpose of performing Mulliken population analysis, we shall transform ${\mathit{\boldsymbol{D}}^{{\text{eff}}}}$(SYM-BCMO) back into that in atomic basis (${\mathit{\boldsymbol{D}}^{{\text{eff}}}}$(AO)) by

where $\mathit{\boldsymbol{C}}$ is the coefficient matrix of SYM-BCMOs obtained by Eq.(14). And population matrix $\mathit{\boldsymbol{P}}$ could then be easily got by calculating each element according to

where ${\chi _{\mu \nu }}$ represents the overlap of atomic orbitals. Then the net charges at atom $i$ can be computed by

where $Z_i$ is nuclear charges of atom $i$.

As for calculating dipole moments properties of the aggregate, we can utilize the dipole moment operator

Similar to effective operators for Hamiltonian and 1-RDM, we can also define ${\hat \mu ^{{\text{eff}}}}$ as the following if we consider only up to two-body interactions

However, here the construction of the dimer model space will be slightly modified because we will also consider the transition dipole moments (TDMs) between the ground state and the excited states. Taking dimer $IJ$ as an example, three states ($\left| {\psi _{IJ}^0} \right\rangle$, $\left| {\psi _{IJ}^1} \right\rangle$ and $\left| {\psi _{IJ}^2} \right\rangle$) are preserved in the target space after a standard quantum chemical calculation and the dipole moments of various states and TDMs in the target space can be easily evaluated by

which will constitute a (3$\times$3)-dimensional matrix ${\boldsymbol\mu _{IJ}}$.

At the same time, the model space becomes span$\left\{ {\left| {\psi _I^0\psi _J^0} \right\rangle , \left| {\psi _I^1\psi _J^0} \right\rangle , \left| {\psi _I^0\psi _J^1} \right\rangle } \right\}$ and thus its project operator ${\hat P_{IJ}}$ should be redefined as

Similarly, one can easily calculate the matrix $\boldsymbol\mu _{IJ}^{{\text{eff}}}$ in model space basis by

where basis transformation matrix is defined as

After obtaining all one-body and two-body effective dipole moment matrices, we can finally build (($N$+1)$\times$($N$+1))-dimensional matrix ${\boldsymbol\mu ^{{\text{eff}}}}$ using Eq.(25). Then the dipole moments of various ground and excited states as well as TDMs between different states for the whole system can be easily calculated according to

C   Computational details

In this work, we take 1D chain systems of ammonia, formaldehyde, ethylene and pyrrole to show the REM strategy of calculating atomic Mulliken charges and TDMs in molecular aggregates at the level of CIS/STO-3G as illustrative simple examples. FIG. 2 shows the aggregate structure derived by translation of a single molecule along a specified direction, and the translational distance is marked. In ammonia aggregates, the ammonia molecules are aligned along one of the N$-$H bonds; in formaldehyde aggregates, the formaldehyde molecules are aligned along the direction perpendicular to one of the C$-$H bonds; in ethylene aggregates, the ethylene molecules are aligned along the C = C bond; in pyrrole aggregates, the pyrrole molecules are aligned along the direction perpendicular to the N$-$H bond. Please note that the intermolecular spacing of H and N or O atom shown as FIG. 2(a) for ammonia aggregates and FIG. 2(b) for formaldehyde aggregates respectively are both 1.85 Å, the typical hydrogen-bond length [30]. The translational distances in FIG. 2 (c) and (d) are empirically chosen to avoid steric hindrance effect between neighboring ethylene and pyrrole molecules. The monomer structures are optimized based on density functional theory (DFT) at the M062X/6-31G(d) level using Gaussian 09 package [31]. CIS and REM-CIS calculations are implemented using our in-house code, and the AO integrals and BCMO generation by monomer HF calculation are obtained using Gaussian 09 package.

Figure  2.  1D molecular aggregates calculated in this work: (a) ammonia, (b) formaldehyde, (c) ethylene, and (d) pyrrole.

DownLoad: Full-Size Img PowerPoint

Ⅲ.   RESULTS AND DISCUSSION

A   Mulliken charge population in delocalized ionic states

As illustrated in our previous works [26-29], REM can be successfully generalized to delocalized ionic states or charge transfer states. Here, we firstly present the Mulliken charge population analysis in delocalized cationic states of 1D ammonia and formaldehyde aggregates (shown in FIG. 2). The calculated ionization potentials (IPs) and the differences of cationic atomic charge are summarized in Table Ⅰ, and the Mulliken charges ($\rho$) of atoms in ammonia and formaldehyde aggregates are displayed in FIG. 3 and FIG. 4 respectively, taking the pentamer as an example.

Table  Ⅰ.  Calculated IPs, the differences of IPs ($\Delta$) and the maximum differences $\Delta\rho_{\max}$ and average differences $\Delta\rho_{\text{average}}$ of cationic atomic charges between CIS with REM-CIS for 1D ammonia and formaldehyde aggregates.

 | Show Table

DownLoad: CSV

Figure  3.  (a) Structure of 1D pentamer of ammonia. (b-d) The Mulliken charges of N and H atoms in (NH$_3$)$_5$$^+$.

DownLoad: Full-Size Img PowerPoint

Figure  4.  (a) Structure of 1D pentamer of formaldehyde. (b-e) The Mulliken charges of O, C and H atoms (CH$_2$O)$_5$$^+$.

DownLoad: Full-Size Img PowerPoint

For ammonia aggregates, it is shown that the differences between IPs calculated by REM-CIS and those from standard CIS are all less than 0.01 hartree. At the same time, one may also notice that these IPs converge rapidly with the increasing oligomer size, implying a localized ionization in 1D ammonia aggregates. Such a charge locality is also verified by the Mulliken charge population analysis, taking an ammonia pentamer, (NH$_3)_5$$^+$ (see FIG. 3) as an example. The net charge dominates mostly at the left-ending monomer, and the largest net positive charge on it also induces a negative charge on the nitrogen atom in the left-ending monomer with a decreased magnitude. This is consistent with our REM wave function analysis, which express the ground state of (NH$_3)_5$$^+$ as a linear combination of

In contrast to ammonia, from Table Ⅰ one may notice that IPs of formaldehyde calculated by REM-CIS show a much slower decrease from (CH$_2$O)$_3$ to (CH$_2$O)$_8$ and have larger differences from those by standard CIS. To explain the phenomenon, we represent the HOMOs (highest occupied molecular orbitals) of (NH$_3$)$_5$ and (CH$_2$O)$_5$ in FIG. 5, which are illustrated using Multiwfn [32]; the monomers Mulliken charges of (NH$_3$)$_5$$^+$ and (CH$_2$O)$_5$$^+$ are represented in FIG. 6. From the HOMO of pentamer, the distribution of ionized electron of ammonia aggregates is mainly from the left-most monomer, while for formaldehyde aggregates, the distribution of ionized electron, besides the left-most monomer, is also from the rest. From the values of positive charge calculated by CIS and REM-CIS, the charge distribution of ammonia aggregates is more localized than formaldehyde aggregates. Therefore, the results by REM-CIS agree with the CIS at least qualitatively. However, because of the two-body approximation, more accurate results calculated by REM-CIS stem from localized systems than delocalized systems.

Figure  5.  HOMO of pentamer by HF. (a) Ammonia pentamer and (b) formaldehyde pentamer.

DownLoad: Full-Size Img PowerPoint

Figure  6.  The Mulliken charges of each monomer. (a) Ammonia pentamer and (b) formaldehyde pentamer.

DownLoad: Full-Size Img PowerPoint

The charge distribution may be affected by two effects, polarization and delocalization. According to chemical intuitions, strong hydrogen bond existed in ammonia aggregates (the typical hydrogen-bond length of 1.85 Å is chosen as molecular spacing), leading to the inhomogeneous charge distribution; despite of typical hydrogen-bond length chosen in formaldehyde aggregates, the atom connected to H atom is the less electronegative C atom. Therefore, the polarization in formaldehyde aggregates is weaker than that in ammonia aggregates. Meanwhile, the parallel $\pi$-orbitals exist in formaldehyde aggregates and cause more delocalized charge distribution. Consequently, these two effects affect the charge distribution together. These effects may be also expressed by two-body inter-block interactions in REM framework. We take interaction of block $I$ and block $J$ as an example. The item of $v_{IJ}^{I^+}$ represents the level of positive charge captured by block $I$ and block $J$ staying neutral, which reflects the level of polarization, and it's similar to $v_{IJ}^{J^+}$. And the item of $h_{IJ}$ reflects the level of positive charge delocalized in block $I$ and block $J$. From the calculated results of effective Hamiltonian, the values of $v_{IJ}^{I^+}$ and $v_{IJ}^{J^+}$ range from $-$0.06 hartree to $-$0.19 hartree for all indexes of $I$ and $J$ in ammonia pentamer, and the maximum presents in the left-most dimer, while the values of $v_{IJ}^{I^+}$ and $v_{IJ}^{J^+}$ range from $-$0.04 hartree to $-$0.13 hartree for all indexes of $I$ and $J$ in formaldehyde pentamer, and the maximum presents in the left-most dimer similarly. The magnitude of $h_{IJ}$ is 10$^{-3}$ hartree for all indexes of $I$ and $J$ in ammonia pentamer, while the magnitude of $h_{IJ}$ is 10$^{-2}$ hartree for all indexes of $I$ and $J$ in formaldehyde pentamer. It is in agreement with our above qualitative analysis.

For the wave function derived by REM-CIS, we utilize two schemes to calculate the reduced density matrix of molecular orbits, in which scheme A applies the operator $a_p^\dagger a_q$ to the model space basis function directly, while scheme B applies the effective operator ($a_p^\dagger a_q$)$^{\text{eff}}$ to the model space basis function at first, then indicate the whole density matrix as the linear combination of density matrix of monomers and dimers, which are introduced in Section II.

The maximum and average differences of cationic atomic charges of ammonia aggregates and formaldehyde aggregates are also represented as Table Ⅰ, and the charge comparison of pentamer between two schemes is represented in FIG. 3 and FIG. 4 respectively. From Table Ⅰ, we may notice that in both ammonia and formaldehyde aggregates, the calculated results by scheme B are closer to the standard CIS than those by scheme A. However, the IPs and charges of formaldehyde aggregates calculated by REM-CIS have larger deviations from standard CIS method than ammonia aggregates. In the fragment-based approaches, there are mainly two approaches for acquiring higher accuracy [25, 33-35]: (ⅰ) consider larger blocks and only block dimers which have two-body effective interactions only, or (ⅱ) use smaller blocks but with considerations of block trimers or even tetramers which have three- or four-body effective interactions. In most cases, the former way is more convenient and efficient. Here, we resort to the former scheme and use (CH$_2$O)$_2$ as a monomer, then (CH$_2$O)$_4$ as a dimer. This new scheme is then called as scheme C. The results of atomic charges of formaldehyde's octamer are represented in FIG. 7. The difference of ionization potentials between REM-CIS$^{\text{B}}$ and CIS is $-$0.02191 hartree, and the difference between REM-CIS$^{\text{C}}$ and CIS is $-$0.01571 hartree. The maximum difference and average difference of atomic charges between REM-CIS$^{\text{B}}$ and CIS are $-$0.1033 and 0.0240 unit charge respectively, while the maximum difference and average difference between REM-CIS$^{\text{C}}$ and CIS are 0.0528 and 0.0170 unit charge respectively. Meanwhile, we can also intuitively observe the improvement by increasing the size of block from FIG. 7.

Figure  7.  (a) Structure of 1D octamer of formaldehyde. (b-e) The Mulliken charges of O, C and H atoms (CH$_2$O)$_8^+$.

DownLoad: Full-Size Img PowerPoint

B   Transition dipole moments

Now we move to test our REM method to the calculations of the lowest singlet excited state S$_1$ of the 1D neutral molecule aggregates at CIS level. Here we kept two states in each monomer, the ground state S$_0$ and the singlet excited state S$_1$ as model states, and took the S$_1$ and S$_2$ excited states as target states for dimers. 1D ethylene aggregates and 1D pyrrole aggregates are selected as the tested systems, whose geometrical configurations are represented in FIG. 2 (c) and (d). These systems show large value of transition dipole moments in a certain direction, thus the calculated systems are ideal systems to test the REM wave function for calculating transition dipole moments. Similarly, we also utilize two schemes to calculate the transition dipole moments, in which scheme A applies the coordinate operator $\hat X$ to the model space basis function directly, while scheme B applies the effective operator ($\hat X$)$^{\text{eff}}$ to the model space basis function at first and then compute the transition dipole moments of whole system as the linear combination of transition dipole moments of monomers and dimers, which introduced in Section II. The values of transition dipole moments between ground state and S$_1$ state in the other directions except the direction in which molecules are arranged in one dimension are zero, so we only list the values in this direction.

Calculated results of pyrrole and ethylene aggregates are listed in Table Ⅱ. For ethylene aggregates and pyrrole aggregates, the maximum differences of excited energies between REM-CIS and CIS is 8.353$\times$10$^{-3}$ and 5.691$\times$10$^{-3}$ hartree respectively, which are in agreement with the standard CIS quite well. And apparently, when we use the REM wave function to calculate the transition dipole moments, the results by scheme B is closer to that of the standard CIS method than scheme A, which shows that scheme B can give more reasonable results based on REM wave function than scheme A in practice again. And the values of transition dipole moments of 1D pyrrole aggregates and 1D ethylene aggregates calculated by REM-CISB are in agreement with the standard CIS quite well.

Table  Ⅱ.  Calculated excited energies $E$, transition dipole moments (TDM) between S$_0$ to S$_1$ and differences $\Delta$ of them between CIS with REM-CIS for linear pyrrole and ethylene aggregates.

 | Show Table

DownLoad: CSV

The correlation diagrams of calculated results of atomic charges and transition dipole moments between CIS and REM-CIS are summarized and shown in FIG. 8. It is evident that, the REM-CIS wavefunction gives a reasonable description of wave function properties of large systems, for both atomic charges and transition dipole moments, with a linear correlation coefficient higher than 0.995.

Figure  8.  Correlation diagrams of calculated charge density (a) and transition dipole moment (b) between CIS and REM-CIS.

DownLoad: Full-Size Img PowerPoint

Before the end, we represent the computational time scaling of ionic 1D formaldehyde aggregates and ethylene aggregates in FIG. 9. FIG. 9 (a) and (b) only represent the time for energy calculations, while FIG. 9 (c) and (d) represent the total time of both energy and property calculations. The computational costs of REM-CIS show a tendency of linear growth basically with the enlargement of monomers, and can be fitted by a linear function shown as FIG. 9, while a polynomial growth is shown for CIS. Therefore, with the enlargement of systems, the REM-CIS show an economic computational costs compared with standard CIS. For CIS, the most computation-intensive part is utilizing the wavefunction derived by CIS to calculate the properties, such as atomic charges and transition dipole moments, and the computational costs are greatly decreased by REM.

Figure  9.  Computational time of 1D aggregates between CIS with REM-CIS, (a) IPs of formaldehyde aggregates, (b) excited energies of ethylene aggregatesm, (c) IPs and atomic charges of formaldehyde aggregates, (d) excited energies and TDMs of ethylene aggregates.

DownLoad: Full-Size Img PowerPoint

Ⅳ.   CONCLUSION

In this work, we proposed a new fragment-based approach to calculate first-order molecular properties (e.g. charge population and transition dipole moment) of delocalized ionic or excited states in molecular aggregates, by integrating the renormalization group theory and the effective operator representation. Testing results for four kinds of 1D molecular aggregates (ammonia, formaldehyde, ethylene and pyrrole) indicate that our new scheme can efficiently describe not only the energies but also wavefunction properties of the low-lying delocalized electronic states in large systems. This provide a useful low-scaling electronic structure approach for excited state studies in various large chemical systems.

In this work, we clarified that the numerical tests are only proof-of-principle studies. The future applications of the approach for realistic molecular systems requires an incorporation of more advanced excited state electronic structure methods (such as TDDFT, SAC-CI and CASSCF) and larger basis set as well as higher-body effect. These works are under way in our group.

Ⅴ.   ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (No.22073045), and the Fundamental Research Funds for the Central Universities.