{'fraction_non_alphanumeric': 0.028053576295604097, 'fraction_numerical': 0.013355137639683837, 'mean_word_length': 4.744464325654216, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long-range Corrected Fragment Molecular Orbital Density-Functional Tight-binding Method for Excited States in Large Molecular Systems

Herein, we present a new method to efficiently calculate electronically excited states in large molecular assemblies, consisting of hundreds of molecules. For this purpose, we combine the long-range corrected tight-binding density-functional fragment molecular orbital method (FMO-LC-DFTB) with an excitonic Hamiltonian, which is constructed in the basis of locally excited and charge-transfer configuration state functions calculated for embedded monomers and dimers and accounts explicitly for the electronic coupling between all types of excitons. We first evaluate both the accuracy and efficiency of our fragmentation approach for molecular dimers and aggregates by comparing with the full LC-TD-DFTB method. The comparison of the calculated spectra of an anthracene cluster shows a very good agreement between our method and the LC-TD-DFTB reference. The effective computational scaling of our method has been explored for anthracene clusters and for perylene bisimide aggregates. We demonstrate the applicability of our method by the calculation of the excited state properties of pentacene crystal models consisting of up to 319 molecules. Furthermore, the participation ratio of the monomer fragments to the excited states is analyzed by the calculation of natural transition orbital (NTO) participation numbers, which are verified by the hole and particle density for a chosen pentacene cluster. The use of our FMO-LC-TDDFTB method will allow for future studies of excitonic dynamics and charge transport to be performed on complex molecular systems consisting of thousands of atoms.


I. INTRODUCTION
In materials science, solid-state physics, biochemistry and many other branches of molecular sciences, the study of excited states is essential for understanding photochemical and photophysical processes as well as the function of molecular materials. Simulations of the excited state properties necessitate an accurate description of the electronic structure of molecular systems, which was first accomplished with abinitio quantum chemical methods 1,2 . However, investigating excited states of extended systems in the framework of the ab-initio approach requires very long calculation times due to the high computational scaling with the system size. With the development of the time-dependent density functional theory (TD-DFT), a reasonably fast method of calculating the excited state properties of relatively large molecular systems (over 100 atoms) became available [3][4][5] . New developments of approximative methods for the prediction of the nature of the excited states, including linearscaling TD-DFT 6 , subsystem density functional theory 7 and multilayer QM/MM approaches (ONIOM) 8,9 , have further increased the applicability of TD-DFT to large polyatomic systems. However, despite the multitude of recent developments, the calculation of the excited states of molecular systems consisting of thousand atoms or more is still not possible even with TD-DFT.
The semiempirical quantum mechanical (SQM) methods, which approximate the wavefunction or density functional based methods, provide a further alternative for investigating the excited states of molecular systems. In the various SQM theories, a minimal basis set is used and the electronic integrals are approximated by the partial or complete neglect of the differential overlap between the atomic ba-sis functions, which leads to a sizable increase in the computational efficiency 10,11 . Whereas these approaches have been studied extensively in the 1970s-1990s [12][13][14][15][16][17] , the development of the self-consistent charge (SCC) density functional tight-binding (DFTB) method [18][19][20] has renewed the interest in semiempirical approaches in the last 20 years. Niehaus et al. 21 adapted the time-dependent density functional theory to the tight-binding formalism, enabling the calculation of excited state properties of large molecular assemblies of hundreds of atoms. The introduction of long-range corrections to the DFTB framework facilitated the investigation of excited state properties of molecular systems involving chargetransfer states [22][23][24] . The DFTB and TD-DFTB methods have been implemented in various software packages including DFTB+ 25 , DFTBaby 26 , ADF 27 , CP2K 28 , and hotbit 29 , which can be employed to predict the excited state properties of polyatomic systems of over 1000 atoms.
Recently, Grimme developed the simplified Tamm-Dancoff approximation 30 (sTDA) for the calculation of excited-state spectra of large molecules, which is based on a ground state DFT calculation, followed by the semiempirical treatment of the integrals involved in the linear-response TD-DFT calculation, similar to the TD-DFTB approach. Furthermore, this method was extended to the tight-binding formalism (sTDA-xTB) 31 by the parametrization of the ground-state Hamiltonian to reduce the computational requirements of the previous models. Further developments of this approach include the GFN-xTB methods by Grimme and coworkers, which can be employed to predict the ground state properties and structures of a multitude of molecular systems [32][33][34] . Additional notable semiempirical methods, which can be used to predict the excited state properties of large molecular systems are MNDO 12 , AM1 13 , PMx 14 and OMx 17 .
The computational efficiency and applicability of all above arXiv:2211.16248v1 [physics.chem-ph] 29 Nov 2022 methods can be further enhanced by using fragmentationbased approaches allowing the simulation of large molecular systems of over a million atoms. Prominent examples include the fragment molecular orbital [35][36][37][38][39][40][41][42][43] (FMO) method and the divide-and-conquer [44][45][46][47][48] (DC) scheme. In these fragmentation approaches, a partitioning scheme divides the system into multiple fragments. The properties of the complete system are then obtained by combining the results of the different fragments. In the case of the FMO method, the energy of the fragments -the monomers and fragment pairs -is calculated iteratively while considering the electrostatic embedding potential (ESP) resulting from the interaction of all fragments. The number of fragment pairs is determined by the closest atomic interfragment distance. The FMO and the DC approach have previously been employed to investigate complex molecular systems, such as proteins 49 , polymers 50 or optoelectronic materials 51 . Furthermore, Nishimoto, Fedorov and Irle combined the FMO method with the density functional tight-binding approach (FMO-DFTB) to enable the study of even larger systems 52,53 . The FMO-DFTB method has been extended to include analytical groundstate gradients 54 , the polarizable continuum model 55 , periodic boundary conditions 56 , the long-range correction 57 and many other theoretical concepts [58][59][60][61] . It has been succesfully employed to study the properties of large molecular systems, such as the charge transport in covalent organic frameworks 62,63 . While these fragmentation-based methods provide excellent tools to calculate the ground state properties of molecular aggregates, the study of excited states is limited to only local excitations within the fragments [64][65][66][67][68] . Recently, several approaches have been developed to investigate the excited states properties of molecular clusters using an excitonic Hamiltonian, consisting of local excitations (LEs) and chargetransfer states (CTs) [69][70][71] . Fujita and Mochizuki combined the multilayer FMO 38 and the transition density fragment interaction 72,73 method and introduced an excitonic Hamiltonian to calculate nonlocal excitations in molecular clusters. A basis consisting of fragment configuration state functions (CSFs), which define the locally excited and charge-transfer states, is used to construct the Hamiltonian 74 .
In this work, we combine the FMO-LC-DFTB method with the construction of an excitonic Hamiltonian to calculate the excited states of large molecular systems consisting of hundreds of molecules. As a basis for the description of the excited states of large molecular assemblies, we use locally excited and charge-transfer states, which are calculated for monomers and pair dimers employing the LC-TD-DFTB approach. While the diagonal elements of the Hamiltonian are represented by the energies of the basis states, the couplings between the basis states, which represent the off-diagonal matrix elements of the Hamiltonian, are calculated by utilizing the tight-binding formalism for the two-electron fourcenter integrals. The excited state energies of the full system are given by the diagonalization of the excitonic Hamiltonian. The accuracy and the efficiency of the proposed theory has been tested by the comparison with the full LC-TD-DFTB method for molecular dimers and large assemblies of molecules. The effective computational scaling was determined by the calculation of large clusters of anthracene and perylene bisimide aggregates. In addition, the applicability of our approach was confirmed by the calculation of the excited state properties of pentacene crystal models. The FMO-LC-TDDFTB method was implemented in our own software package DIALECT 75 , which is publicly available on Github and was written in the Rust programming language.
Although Rust is still a very young language -with the first stable version (1.0) released in 2015 -its potential for scientific computing applications becomes increasingly apparent 76 . The Rust programming language does not only have a strong focus on the code and memory safety, but also shows performance comparable to the C, C++ and Fortran languages. Additionally, Rust offers in our opinion a friendlier syntax and provides a strongly growing repository of development tools. Therefore, we expect Rust to gain high popularity also in the field of quantum chemistry in the coming years.
The present work is structured as follows: In Sec. II, the methodological framework of our approach is outlined. Subsequently, in Sec. III, the computational procedure and the different steps of a FMO-LC-TDDFTB calculation are illustrated. The application of the proposed theory to pentacene and the results of the benchmark calculations regarding the accuracy and the effective scaling of the FMO-LC-TDDFTB method are presented in Sec. IV. Finally, the conclusions and outlook are given in Sec. V.

II. METHODOLOGY
In the following, we will present a general formulation of the fragment molecular orbital (FMO) LC-TD-DFTB. As a starting point, we review the theory behind LC-TD-DFTB 23,24 . This is followed by the formulation of the FMOansatz for the description of the electronic ground state. Subsequently, we present how the molecular orbitals of the total system can be calculated and introduce a theoretical approach for calculating electronically excited states. This is concluded with a description of the necessary Hamiltonian matrix elements and the computational procedure that is used in our implementation. Throughout the paper we use atomic units and the following notation convention: atoms are denoted as uppercase letters from A-H and molecular fragments from I-Z without indices. Uppercase letters with indices denote matrix elements and bold uppercase letters matrices. For molecular orbital (atomic orbital) indices lowercase (greek) letters are used. It should be noted that in the following we use the terms fragment and monomer synonymously to refer to a molecular unit. Likewise, we use (molecular) cluster and aggregate synonymously to refer to a large number of molecular units that are not covalently linked.
A. Long-range corrected DFTB (LC-DFTB) As shown previously, the formalism of LC-DFTB can be derived by a Taylor expansion of a long-range corrected full density functional (e.g. LC-PBE) around a reference density that is given by a sum of atomic densities, as described in detail elsewhere 23,24 . In the following, we will only provide a brief summary of the working equations. After employing the usual tight-binding (TB) approximations 20 the total energy in SCC LC-DFTB is given by: Here P µν denotes the electron density matrix elements and H 0 µν the one-electron Hamiltonian matrix elements. The electron density difference matrix is defined as ∆P = P − P 0 , where P 0 is the diagonal reference electron density matrix. The deviation of the Mulliken charge q A on Atom A from the charge of the neutral Atom q 0 A is given by where S µν represents a matrix element of the overlap matrix.
Applying the tight-binding approximations to the twoelectron integrals in the Coulomb and exchange part of the energy and neglecting all 3-and 4-center integrals gives the following expressions: with the transition charges on atom A (in the atomic orbital basis) defined as The γ-matrices for charge fluctuation interactions can be calculated assuming that the charge fluctuations can be represented by spherical Gaussian functions, leading to where R is the distance between the atomic centers A and B and C AB = (2 σ 2 A + σ 2 B ) − 1 2 depends on the widths σ A and σ B of the charge clouds on the two atoms. The widths are determined by the atom-specific Hubbard parameters U A as In LC-DFTB the Coulomb potential is separated into a long-range and a short range part where the posi-tion of the smooth transition between the two regimes is controlled by the long-range radius R lr The short range term is already included by the local exchange-correlation functional employed in LC-DFTB. The electron integrals of the screened Coulomb potential for the long-range contribution can be calculated as The long-range γ-matrix is defined similarly to the γ-matrix in eq. 5, where depends on the range-separation parameter R lr . The transformation of the two-electron repulsion integrals into molecular orbital (MO) basis leads to the following expressions: and The atom-centered transition charges between MOs are defined as: Variational minimization of the LC-DFTB energy with respect to the molecular orbitals leads to the corresponding Hamiltonian, that is defined as: The ground state MOs and their energy are obtained by a self consistent procedure for solving the following general eigenvalue problem Being a semiempirical theory, (LC)-DFTB is heavily dependent on a parametrization for the electronic Hamiltonian, H 0 µν , the overlap matrix S µν and the repulsive potentials V rep AB . The parametrization usually starts with the computation of pseudoorbitals, the tabulation of Slater-Koster files and ends with the fitting of repulsive potentials 29,77,78 . A benchmark of the used LC-DFTB parametrization for organic and biological molecules is given in Ref. 79 . The fragment molecular orbital method in combination with DFTB was developed by Nishimoto, Fedorov and Irle 52 . Recently, it was adopted to LC-DFTB by Niehaus and Irle 57 and this theory forms the basis for our ground-state calculations. Therefore, we will only restate the most important equations in the following and refer to Ref. 52 and Ref. 57 for a more detailed description. Only the case of molecular clusters -molecules that are not covalently connected -is considered in this work, and therefore, we do not introduce the hybrid orbital projection operator 36 . The total energy of the entire system in the self-consistent FMO-LC-DFTB is given by Here E I (E J ) are the energies of fragment I (J) and E IJ is the energy of the pair IJ. The total number of fragments is denoted as N, while we will use N I for the number of atoms in fragment I. The energies E X where X is either I or IJ can be further separated into a part that depends only on the isolated molecule and another one that accounts for the environment of the molecule: The former is the so-called internal energy, E X , which is equivalent to the total energy calculated in conventional LC-DFTB of entity X. The embedding energy, E em X , accounts for the Coulomb interaction between the entity X and all other fragments: Inserting eq. 16 and eq. 19 in eq. 15 allows to rewrite the total FMO energy as: The last term on the r.h.s. is the difference in embedding energy of the pair and corresponding embedded but noninteracting monomers.
Here ∆∆q IJ A denotes the difference between charges of pair IJ and charges of fragments I and J for atom A: We employ the electrostatic dimer (ES-DIM) approximation 37 for pairs, in which the fragments are so far separated, that their orbital overlap will be zero. The energy of the far-separated pairs is given by It should be noted that although it is possible to separate the total energy into internal and embedding energies, the tightbinding Hamiltonian of an entity (cf. eq. 13) also depends on the electron density or, more specifically, the charges of all the other fragments where, V X µν , is the electrostatic potential (ESP) that acts on fragment X and is defined as: This makes it necessary to perform the SCC iterations for fragments self-consistently so that one iteration step per fragment is made to update all charges simultaneously. The obtained self-consistent charges are afterwards used for the calculation of the SCC iterations on pairs for which the ES-DIM approximation is not used.
Calculating the electronic ground state in the context of FMO(-DFTB) allows one to make accurate calculations of the energy, but one does not obtain MOs of the entire system. Instead, one obtains only the MOs for each of the fragments and all fragment pairs that are not affected by the ES-DIM approximation. Since they are obtained by independent calculations such MOs are not orthogonal. However, it is possible, as shown by Tsuneyuki et al. 40 , to construct orthogonal MOs for the whole system and to transform the Hamiltonian into such basis by using Löwdin's approach 80 Since the computation of the matrix inverse of the total overlap is computationally demanding and the overlap matrix in the basis of the fragments and pairs is almost diagonal, we approximate it in first order by where the overlap matrix matrix elements are given by and ϕ I p is the p-th MO of fragment I. The Hamiltonian of the whole system, H LCMO , can be constructed from the Hamiltonian matrices of the corresponding fragments and pairs as The ⊕ sign indicates that each block in the total Hamiltonian matrix is filled with the Hamiltonian of the corresponding fragment (H I ) or fragment pair (H IJ ).
In contrast to the diagonal Hamiltonian matrix of the fragments where ε a is the a-th MO energy of fragment X, the Hamiltonian of pair IJ accounts for the fact that the MOs between the fragments might not be strictly orthogonal where the sum runs over all MOs of the pair IJ. Of course, it would now be possible to determine the eigenvectors of the Hamiltonian and thereby obtain the molecular orbitals of the entire system. In principle, once the MOs of the whole system are available, the excited states could be calculated using standard TD-DFTB procedure. While this seems straightforward, the disadvantage is that the calculation of the excited states are as costly as in a non-FMO approach and one would lose the scalability of the method.

C. Excited states in the frame of FMO-LC-DFTB
An alternative way to compute excited states in the framework of the FMO-Ansatz is motivated by the idea that an excited state wavefunction of the whole system can be expressed as a linear combination of basis states of the respective fragments and fragment pairs. Since we restrict ourselves to the interaction between two fragments in FMO, we also use only locally excited (LE) states on one fragment and charge-transfer (CT) states between two fragments as basis states. The electronic excited state wavefunction, |Ψ , of a large molecule assembly can then be expressed as where the coefficients for fragment LE and CT configuration state functions (CSFs) are given by c m I and c m I→J , respectively. The coefficients are obtained by solving the eigenvalue problem Hc = Ec. Notice that we assume that the basis states are orthogonal such that S = 1. However, in the calculation of the Hamiltonian matrix elements the overlap matrix is explicitly included. This is in line with other semi-empirical methods that employ the same approximation. On the basis of FMO Hartree-Fock theory this approach was introduced by Fujita and Mochizuki 74 . However, we will show that our definition of the CT states differs and how this impacts not only the energies but also the scaling of the method. The ansatz in eq. 32 allows one to choose the energetically lowest basis states by setting the number of LE states per fragment (N LE ) and the number of CT states per pair (N CT ) to appropriate values. The number of eigenvalues for the corresponding eigenvalue problem then is given by N · N LE + (N 2 − N) · N CT . Furthermore, since these basis states are per definition locally excited or charge-transfer states the adiabatic states' eigenvectors (cf. 32) contain the state character information. In principle, it would also be possible to extend this approach to doubly or even higher excited states. The LE states, |LE m I , are calculated as singly excited states of a fragment. We restrict ourselves in this work to singlet states, but it should be noted that this approach can be extended to triplet excited states in a straightforward way. The m-th excited state (S m ) on fragment I defines the following LE basis state: Here, a † Ia,α and a Ia,β are the creation and annihilation operators of the α and β spin electron in the ath orbital of fragment I, and |G is the ground-state wavefunction. T m(I) is the oneparticle transition density matrix of the m-th excited state of fragment I in the MO basis and |Φ i→a I is the configuration state function of the excitation from the occupied orbital i to the virtual orbital a on fragment I. Similarly, the m-th intermolecular CT state from fragment I (hole) to fragment J (particle) can be defined as: The transition between the fragments I and J is restricted to the occupied orbitals of monomer I and the virtual orbitals of monomer J. The definition of the CT state differs to the approach of Fujita and Mochizuki 74 as they use the following single configuration state function as a CT state Both approaches have their merits, and we will briefly explain below why we define CT states as shown in equation 36. The latter necessitates solving the CIS matrix for the lowest N CT charge-transfer states of each pair. Therefore, the construction of the basis states will no longer scale linearly with the number of fragments but quadratically. Nevertheless, this approach has the advantage that all electronic couplings between CT states on a pair are omitted since these are eigenstates of the pair. Furthermore, in both approaches, a restriction of the considered CT states or occupied and virtual orbitals is necessary. Using equation 36, we can choose to use the N CT lowest energy CT states and ensure that all relevant orbitals for these states are considered, whereas with eq. 37, the selection of states is made by considering N hole occupied and N particle virtual MOs. A higher number of included orbitals quickly leads to very large eigenvalue problems as the number of considered CT states scales with (N 2 − N) · N hole · N particle . We will show in section IV A how the accuracy of the excitation energies and the computational scaling depend on using either eq. 36 or 37 for the definition of charge-transfer states and that the usage of eq. 36 is generally preferable.

D. Hamiltonian matrix elements
The definition of basis states in eq. 34 and 36 allows to express the matrix elements of the system Hamiltonian. We will first focus on the excitation energies (diagonal matrix elements) of LE and CT states. The energies are calculated by a LC-DFTB calculation of the fragment (pair) in the framework of linear-response LC-TD-DFT, which was adapted to tight-binding DFT by Niehaus et al. 21 and has been extended to include the LC correction by some of us 23 . In linearresponse TD-DFT excited states are computed by solving the non-Hermitian eigenvalue problem In the framework of LC-TD-DFTB the matrices A and B are defined as 23 where ε i is the energy of the i-th molecular orbital. The orbital energies ε i are represented by the diagonal elements of the orthogonalized total Hamiltonian H ii . Using the Tamm-Dancoff (TDA) approximation (B = 0), eq. (38) is reduced to a simple Hermitian eigenvalue problem which can be solved iteratively for the lowest eigenvalues by employing the algorithm developed by Davidson 81 . This is in principle equivalent to the configuration interaction singles (CIS) procedure. The energies of the CT states are calculated in analogy to the local excitations, where the excitations are restricted to the transition from the occupied orbitals of one monomer to the virtual orbitals of the other one.
To calculate the off-diagonal matrix elements of the Hamiltonian, the couplings between the basis states are required. In accordance with Ref. 74 , the Hamiltonian can be split into one electron and two-electron contributions which were derived using the usual Slater-Condon rules. The two-electron coupling matrix elements of the Hamiltonian require the calculation of the four-center electron repulsion integrals (cf. eq. 43). In contrast to the work of Fujita and Mochizuki 74 , the off-diagonal elements between the basis states that involve three or four fragments are not neglected.

LE-LE matrix elements
The coupling between locally excited states is given only by the two-electron part, as the one electron coupling between locally excited states vanish due to the Slater-Condon rule, where q m(I) is the transition charge of the m-th excited state of the fragment I on atom A. q i j A and q ab B are the atom-centered transition charges between the occupied and virtual orbitals of the two monomer fragements I and J (cf. eq. (12)). The first term of the LE-LE coupling (cf. eq. 44) characterizes the Coulomb interaction between the transition densities of the respective excitations. The second term represents the exchange interaction. In the case of LE states on far separated fragments, fragment pairs for which the ES-DIM approximation is used, the exchange contribution can be neglected as the overlap between both states will be zero: Therefore, only the long-range Coulomb interaction is taken into account in this case.

LE-CT matrix elements
The Hamiltonian matrix elements between an LE state on fragment I and a CT state which is formed by linear combination of holes on fragment J and electrons on fragment K contains a one-electron part as well as a contribution from the two-electron interaction The first and second term in eq. 50 correspond to the Coulomb and exchange interaction between the transition densities of the LE and CT state. The one-electron contribution is only non-zero if the CT state shares one of its fragment with the LE state and if the other fragment is in spatial proximity, so that the ES-DIM approximation is not used, because otherwise the matrix element H pq will be zero. If both fragments of the CT state are far apart (ES-DIM approx.) from the fragment of the LE state then the exchange term of eq. 50 will be zero and is neglected. In this case the matrix element becomes and only the Coulomb interaction is calculated.

CT-CT matrix elements
The one-electron CT-CT coupling vanishes for I = K or K = L (cf. eq. 42), and thus, is reduced to the diagonal contributions, which are included in the LC-TD-DFTB calculation of the CT state. The two-electron off-diagonal coupling between two charge-transfer states is given by As shown for the LE-LE and LE-CT matrix elements, we also neglect the exchange part of the coupling for the CT-CT couplings in case if either fragment I and K or J and L are different fragments that are far apart. Thus in this case the matrix element simplifies to

III. COMPUTATIONAL PROCEDURE
Herein, the sequence of the steps necessary to carry out a FMO-LC-TDDFTB calculation is briefly summarized. The In the first step, the input geometry is partitioned into the different monomer, pair and ES-DIM pair fragments. The ES-DIM approximation 37 is applied if the closest distance between two fragments exceeds the threshold value of twice the sum of the van der Waals radii of the closest atoms. Subsequently, the coupled lc-DFTB SCC iterations of the monomer fragments are performed until the convergence of all monomer calculations is achieved. As the electrostatic potential between all monomer fragments is required to calculate the monomer Hamiltonian (cf. eq. 24), the SCC iterations for all fragments are computed simultaneously in order to update the charges of all monomers in each step. The final self-consistent charges are stored, as they are necessary for the calculation of the pair, ES-DIM pair and embedding energies. The Anderson acceleration is used to speed up the convergence of the SCC routine 82,83 .
Then, using the electrostatic potential between all fragments (cf. eq. 25) obtained from the previous step, the pair energies are calculated. The charges of the pair fragments are stored so that they can be used in the calculation of the embedding energy. The ES-DIM pair energy (cf. eq. 23) and the embedding energy (cf. eq. 19) are determined, using the charges of the monomer and pair fragments obtained from the previous SCC routines. The results of the fragment calculations are combined according to eq. (18), which yields the ground-state energy of the FMO-LC-TDDFTB approach.
Thereafter, the LCMO-Hamiltonian matrix is built according to eq. 26-31, using the orbital energies, orbital coefficients and overlap matrices of the fragments. The H matrix is obtained by subsequent Löwdin orthogonalization. In order to construct the LE and CT basis states, LC-TD-DFTB excited state calculations are performed for both the monomers, the pairs and the ES-DIM pairs, where the orbital energies of the respective fragments are replaced by the diagonal matrix elements of the H matrix.
Subsequently, the off-diagonal matrix elements of the Hamiltonian are calculated according to the expressions for the couplings between the LE and CT states given in eq. 45 -54. As the calculation of the exchange terms in the eq. (45), (50) and (53) includes all possible transitions between the orbitals of the monomer pairs, a screening of the one-particle transition density matrices T was introduced, where only matrix elements over a certain threshold are considered. Thus, the number of transitions is limited to the actually contributing occupied and virtual orbitals of the fragments.
At the end of the FMO-LC-TDDFTB calculation, the Davidson diagonalization is utilized to obtain the excited state energies and coefficients, which are then used to calculate the oscillator strengths.
The LC-TD-DFTB and FMO-LC-TDDFTB methods were implemented in our own software package DIALECT, which is available on Github 75 .
All calculations regarding the benchmarks and scaling tests were performed on a single core of a computing node with dual E5-2680 Xeon CPUs (2.4 GHz, 14 cores each), 188 Gb of DDR4 memory and a SATA hard drive.
The DFTB parameter set ob2-split 79 was employed in all calculations. In the case of the benchmark calculations, a long-range radius of 3.03 a 0 was used. In order to compen-sate the overestimation of the charge-transfer energies in the pentacene clusters, IP-EA tuning 84 was employed and an optimized long-range radius was used for these systems. This procedure will be described in greater detail in section IV C.
The Mercury 85 program was employed to generate the various clusters of the anthracene, pentacene and perylen bisimide systems from their respective crystal structures [86][87][88] .
In case of the scans of potential energy curves of the pyrene dimer, the DFT-D3 dispersion correction 89,90 was calculated using the simple-dftd3 91 program and the ob2-split dispersion parameters 79 .

A. Accuracy of FMO-LC-TDDFTB excited states
The accuracy of the FMO-LC-TDDFTB methodology is evaluated calculating the excitation energies of π-stacked dimers of pyrene. The potential energy curves of the first 6 excited states were calculated as a function of the π-π stacking distance R z and the parallel shift coordinate R x . As a reference we performed calculations with the conventional LC-TD-DFTB approach. The main purpose of this comparison is to test the FMO approach and the validity of the employed approximations and not to benchmark the general performance. The latter has been thoroughly benchmarked in previous works. 79,92,93 In the case of the FMO-LC-TDDFTB calculation, 20 locally excited states for each monomer and 15 charge-transfer states for each pair were included. In Fig.  2a, the electronic state energies of the FMO-LC-TDDFTB and LC-TD-DFTB methods are shown for a scan of the π-π stacking distance R z . While there is a slight deviation from the LC-TD-DFTB energies at distances from 2.75 Å to 3.0 Å, the deviation between the two methods decreases with increasing interfragment distances. From a distance of 3.5 Å onwards, the results of the FMO-LC-TDDFTB calculation show an excellent agreement with the LC-TD-DFTB energies. Below distances of 2.75 Å, which already belong to the repulsive part of the potential energy curves, the deviation from the reference LC-TD-DFTB calculation results from the non-orthogonality of the MOs. Fig. 2b shows the potential energy curves of the pyrene dimer for a scan of the parallel shift coordinate R x at a π-π stacking distance of 3.1 Å. Here, the results of the FMO-LC-TDDFTB calculation are in good agreement with the LC-TD-DFTB reference.
As the excited state Hamiltonian is constructed using LEs and CTs as basis states, the analysis of the excited state composition becomes easily accessible. The contributions of the LE and CT states to the first 4 excited states for different parallel shift distances are shown in Fig. 2c. As the shift along R x decreases, the contribution of the charge-transfer states increases. While the S 1 , S 2 and S 3 state show CT character of up to 30%, 20% and 12%, the S 4 onlyconsists of local excitations for all R x values.
In addition, the influence of the number of basis states on the accuracy of the excited states energies was investigated. To this end, the mean absolute errors (MAEs) of the first 6 excited states were calculated for different numbers of LE and CT states. Table I   As expected, increasing the number of basis states improves the accuracy of the excited states energies. However, Table I shows that a further increase from 20 LE and 15 CT to 30 LE and 20 CT states yields only a marginal improvement of the MAEs as the additional states have a negligible electronic coupling to the first 6 excited states. At all physically relevant intermolecular separations, the FMO-LC-TDDFTB method al-most perfectly reproduces the PES for the pyrene dimer. The deviations start to be visible at distances around 3 Å. Below 2.75 Å, which is already in the strongly repulsive part of the potential, the non-orthogonality of the MOs causes the deviation from the results of the LC-TD-DFTB reference as the error in the couplings between the LE and CT states increases. Nevertheless, at distances at around 2.75 Å the MAEs are still relatively small; thus, the model should be sufficient for most physically realistic situations.
To evaluate the accuracy of our implementation of the charge-transfer states (cf. eq. 36) in comparison to the approach of Fujita et al. 74 (cf. eq. 37), the mean absolute errors of the first 40 excited states from the LC-TD-DFTB reference of a π-stacked system of four pyrene monomers were calculated within both implementations. Fig. 3 shows the results for stacking distances of 3.5 and 4.0 Å. For both methods, the number of LE states was set to 10 and the number of CT states was varied. As explained in section II C, in our implementation, the dimension of the Hamiltonian grows linearly as we consider the lowest N CT states, and in the case of of Fujita's and Mochizuki's approach, the dimension of the Hamiltonian grows quadratically as all transitions between N hole occupied and N particle virtual orbitals are considered as separate CT states. As the amount of the charge-transfer states is increased for both implementations, the MAEs decrease in value. While both approaches achieve approximately the same accuracy, our method necessitates a much smaller size of the Hamiltonian.
In order to investigate the accuracy of the FMO-LC-TDDFTB method regarding the energies and oscillator strengths of larger systems, we compared the absorption spec- tra of an anthracene cluster composed of 48 fragments (1152 atoms) including the first 40 excited states for both theoretical methods. The basis was constructed using two LE states for each monomer and one CT state for each pair. The individual vertical transitions were convolved using a Lorentzian line-shape function with a width of 2 meV. The comparison between the calculated spectra is depicted in Fig. 4 and shows a slight blue shift of the FMO-LC-TDDFTB energies by about 20-30 meV. However, since the general error in TD-DFT and, in particular, in LC-TD-DFTB excitation energies is on the order of 0.5 eV 79,93 , the additional error due to the FMO approach is hardly significant. Aside from slightly underestimating the intensity of the first peak, our approach could sufficiently reproduce the anthracene cluster's LC-TD-DFTB absorption spectrum.

B. Computational scaling of FMO-LC-TDDFTB
To investigate the computational efficiency of the developed method for the calculation of excited states, we compare the wall times of the method against conventional LC-TD-DFTB for anthracene clusters containing up to 360 fragments (8640 atoms). In the case of the LC-TD-DFTB approach, the number of fragments was limited to 48 (1152 atoms) due to the rapid increase in computational demand. For each cluster, the first 40 excited singlet states were calculated, and as in the previous anthracene calculation (cf. Fig. 4) two LE states per monomer and one CT state per pair were used for FMO-LC-TDDFTB computations. The wall times of the FMO-LC-TDDFTB and LC-TD-DFTB calculations are compared in Fig. 5. In the case of the anthracene cluster containing 48 fragments (1152 atoms), a speedup factor of approx. 1.5×10 4 was achieved. Whereas the LC-TD-DFTB calculation, using a full active space, took almost 11 days, the fragment approach finished in 60 seconds. Furthermore, we analyzed the timings of the LC-TD-DFTB method using only 20 % of the active orbital space. As expected, this reduction in size of the A matrix in the TDA-DFTB approach significantly accelerates the computation of excited states for the anthracene cluster consisting of 48 fragments which took 6293 seconds. However, the FMO-LC-TDDFTB approach is still around 100 times faster for this system and the speedup grows to factor of more than 200 for the next bigger anthracene cluster containing 70 fragments (1680 atoms). It should be noted that a restriction of the active spaces only lowers the prefactor of the method and not scaling. In addition, the restriction of the active orbitals leads to a loss of accuracy in (LC)-TD-DFTB as shown in Fig.  6. While the energies and oscillator strength of FMO-LC-TDDFTB with only two LE states and one CT state has shown a good agreement with LC-TD-DFTB, the reduction of the active orbital space to 20 % leads to significantly larger errors of around 150 meV. In order to gain an insight into the time requirements of the different steps involved in a FMO-LC-TDDFTB calculation, the total wall times have been split into separate parts. Fig. 7a shows the timings of five steps: In addition, the asymptotic scaling of the respective parts of the calculation was evaluated. With the increase in the size of the anthracene clusters, the Hamiltonian size grows up to a dimension of 129960 2 for the system of 360 fragments (8640 atoms). Thus, the most time demanding steps are the calculation of the electronic couplings between the basis states and the following diagonalization for the lowest eigenvalues, which scale with factors of O(N 3.6 ) and O(N 5.1 ), respectively. The SCC routine implemented in our code shows the same nearly linear scaling of O(N 1.1 ), which has been reported previously 52 . As the construction of the Fock matrix is limited by the scaling of the matrix multiplications during the Löwdin orthogonalization, it scales with a factor of O(N 2.9 ). Compared to the other parts of the calculation, the scaling of the construction of the basis states is negligible with a factor of O (N 1.9 ). The combination of all steps of the procedure results in an effective scaling of O(N 3.4 ) for the total run time.
A molecular system involving larger fragments was investigated to verify the effective scaling of the FMO-LC-TDDFTB method. To this end, different clusters of substituted perylene bisimides (PBIs), where each monomer contains 142 atoms, were created, and the first 40 singlet excited states were calculated. Herefore, a basis of 2 locally excited states for each monomer and one charge-transfer state for each pair was employed. As shown in Fig. 7b, the asymptotic scaling of the FMO approach remained approximately the same, if the Davidson diagonalization and the electronic couplings are excluded. Since the PBI clusters consist of distinctly less monomers than the anthracene clusters, the Hamiltonian is also significantly smaller. Subsequently, the fit of the scaling, which is strongly determined by the timings of the bigger clusters, results in smaller factors for the PBI systems, and thus, the total run time only scales with a factor of O(N 2.2 ). As the system size grows above 100 fragments (approx. 15000 atoms), the calculation of the Hamiltonian becomes the most time consuming factor of the whole calculation. Subsequently, a further increase in the number of fragments would also result in a higher effective scaling factor as the Hamiltonian grows and the calculation of the couplings becomes the dominant aspect.

C. Application to Pentacene
After demonstrating that FMO-LC-TDDFTB can simulate excited states and absorption spectra with significantly reduced effort and, in return, at nearly the same accuracy as LC-TD-DFTB, we wish to show the scope of possible applications by performing calculations on large pentacene clusters which serve as a model for bulk pentacene. This system has recently gained a lot of experimental and theoretical attention since it serves as a prototype for organic materials whose properties are determined by an interplay between local (Frenkel) and CT excitations. [94][95][96] The FMO-LC-TDDFTB method was used to simulate the absorption spectrum of pentacene clusters of different sizes, which were generated from crystal structure  data. As a first step, an adjustment of the long-range radius was required since the energies of the pentacene aggregates' charge-transfer states are overestimated with the default radius of 3.03 a 0 . Subsequently, ionization potential and electron affinity (IP-EA) tuning 84 were employed to estimate a reasonable value for the long-range radius by minimizing the following expression A scan of the radius yielded a minimum value of 6.875 a 0 as the optimal long-range radius which was used for all subsequent pentacene simulations. The ionization potential of 6.59 eV and the electron affinity of 1.35 eV were taken from experimental data of pentacene 97,98 . A value of 6.67 eV and 1.77 eV was obtained using the optimally tuned radius for ε HOMO and ε LUMO , respectively. For each of the pentacene clusters, which contain up to 319 fragments (11484 atoms), a calculation of the singlet excited state energies was performed. The basis states were constructed using three LE states for each monomer and one CT state for each pair. Depending on the system size, up to 3000 excited states have been calculated. The comparison between the experimental spectrum of a pentacene film 99 and our results is shown in Fig, 8. The calculated oscillator strengths of the different pentacene clusters have been convolved with a Lorentzian line shape function of 10 meV width. In addition, the spectra were red-shifted by 0.7 eV to compensate for the error of the absorption energy. In the literature similar shifts are used even in the case of TD-DFT. 69 Overall, the calculated oscillator strengths for the first two bands of the absorption spectrum are in good agreement with the experimental spectrum. As the size of the pentacene clusters increases, the intensity of the first band grows in strength, while the oscillator strength of the second band shows a slight decrease in magnitude. Comparing the third and fourth absorption bands of the spectra, the deviation in the energy and intensity are apparent. While the intensity of the third peak is overestimated by a factor of two, the magnitude of the forth band is severely underestimated by FMO-LC-TDDFTB. However, these discrepancies may partially stem from the fact that the influence of the vibrational modes is fully neglected in the present theory, which play a significant role in the case of the pentacene clusters. [102][103][104] Studying the excited states of large molecular systems leads to the question of how many monomers in a cluster contribute to an excited state. To this end, the participation numbers of the natural transition orbitals (NTOs) 105,106 of the excited states of the pentacene clusters were calculated according to the expression where λ i are singular values which were obtained by a singular value decomposition (SVD) of the transition density matrix between the electronic ground state and a singlet excited state 107 . Due to the fact that the pentacene monomers in the crystal structure all share the same internal coordinates, the locally excited states of the various monomers are nearly degenerate, which results in unrealistically high participation numbers for a perfectly arranged structure. In order to model the structural disorder that is present in a real system, we optimized the monomer geometry in the ground state using LC-DFTB and sampled coordinates of the monomers from the ground-state canonical harmonic Wigner distribution. 108 Subsequently, the NTO participation numbers were calculated for the singlet excited states. In order to interpret the data, which showed strongly fluctuating participation numbers of the various states, the simple moving average was calculated by convolving the data with a rectangular window function. To guarantee a consistent approach, the width of the window was chosen to be 20% of the complete number of data points of the respective system. The resulting participation numbers of the different pentacene clusters are shown in Fig. 9. As expected, as the size of the pentacene clusters increases, the number of monomers contributing to the excited states grows. However, a convergence of the participation numbers would also be expected after a specific system size is reached. In the case of the pentacene aggregates, the increase in the average participation number seems to be slowing down gradually.
In order to visualize electron and hole delocalization the particle and hole density of a chosen pentacene system, containing 99 fragments (3564 atoms), was computed. The densities, which are shown in Fig. 10, were selected due to their relative large participation number of 22.7. Both local excitations and charge-transfer states are discernible.
Thus, we have proven the applicability of our method to large molecular systems. The FMO-LC-TDDFTB method is sufficient to estimate the excited state properties of molecular aggregates and can be used to analyze the composition of the excited states.

V. CONCLUSIONS AND OUTLOOK
In this work, we have developed a new method to calculate the excited states in large molecular assemblies, consisting of hundreds of molecules. To this end, we based our approach on the fragment orbital density-functional tightbinding method 52,57 (FMO-DFTB) and employed an excitonic Hamiltonian constructed from locally excited and charge-transfer configuration state functions. Our method has been implemented in the software package DIALECT, publicly available on Github 75 .
We have proven that the accuracy of the proposed theory is sufficient to reproduce the excited state LC-TD-DFTB potential energy curves of π-stacked molecular pyrene dimers. The mean absolute errors are below 80 meV at all physically relevant intermolecular separations for the pyrene systems, showing that for distances of 2.75 Å and above, our method produces satisfactory results. The applicability of the present theory was confirmed by the calculation of the excited state energies and oscillator strengths of a large anthracene cluster. The spectrum of the FMO-LC-TDDFTB approach was in good agreement with the LC-TD-DFTB reference.
The comparison of the wall times of the LC-TD-DFTB and FMO-LC-TDDFTB calculation showed excellent speedup factors for different anthracene cluster sizes of up to 1.5×10 4 . In addition, the computational effective scaling of the proposed theory was evaluated for molecular aggregates involving differently sized monomers. The results showed that in the case of small monomers, the scaling is highly dependent on the calculation of the matrix elements of the Hamiltonian. However, if larger monomers like the substituted perylene bisimide are to be considered, the calculation of the matrix elements of the Hamiltonian only becomes the dominant factor of the effective scaling after the systems grows to a size of over 15000 atoms. The effective scaling of the anthracene and PBI systems was determined as O(N 3.4 ) and O(N 2.2 ), respectively.
At last, the application of our method to pentacene crystals confirmed its validity to large molecular systems. The calculated spectra of the pentacene clusters showed reasonable agreement to the experimental reference spectrum. We were able to analyze the number of monomers that contribute to the excited states of the system by calculating the NTO participation numbers, which were verified by the particle and hole density of a chosen pentacene cluster.
In the future, we plan to extend our approach to quantumclassical dynamics simulations of large molecular systems, like organic semiconductors or optoelectronic materials. To this end, we will implement the analytical gradients of the excited states to enable the investigation of exciton dynamics and charge transport simulations. Concerning the efficiency, graph theory methods will be used to accelerate the diagonalization of the Hamiltonian by partitioning the matrices involved in the Davidson algorithm into separate blocks, and thus reduce the dimensionality of the problem. In addition, we intend to apply the parametrization of the DFTB framework to a wider range of elements to facilitate the investigation of different molecular systems.

CONFLICTS OF INTEREST
There are no conflicts to declare.