Understanding excited state properties of host materials in OLEDs: simulation of absorption spectrum of amorphous 4,4-bis(carbazol-9-yl)-2,2-biphenyl (CBP) †

4,4-Bis(carbazol-9-yl)-2,2-biphenyl (CBP) is widely used as a host material in phosphorescent organic light-emitting diodes (PhOLEDs). In the present study, we simulate the absorption spectra of CBP in gas and condensed phases, respectively, using the eﬃcient time-dependent long-range corrected tight-binding density functional theory (TD-LC-DFTB). The accuracy of the condensed-phase absorption spectra computed using the structures obtained from classical molecular dynamics (MD) and quantum mechanical/molecular mechanical (QM/MM) simulations is examined by comparison with the experimental absorption spectrum. It is found that the TD-LC-DFTB gas-phase spectrum is in good agreement with the GW-BSE spectrum, indicating TD-LC-DFTB is an accurate and robust method in calculating the excitation energies of CBP. For the condensed-phase spectrum, we find that the electrostatic embedding has a minor influence on the excitation energy. Computing accurate absorption spectra is a particular challenge since static and dynamic disorders have to be taken into account. The static disorder results from the molecular packing in the material, which leads to molecule deformations. Since these structural changes sensitively impact the excitation energies of the individual molecules, a proper representation of this static disorder indicates that a reasonable structural model of the material has been generated. The good agreement between computed and experimental absorption spectra is therefore an indicator for the structural model developed. Concerning dynamic disorder, we find that molecular changes occur on long timescales in the ns-regime, which requires the use of fast computation approaches to reach convergence. The structural models derived in this work are the basis for future studies of charge and exciton transfer in CBP and related materials, also concerning the degradation mechanisms of CBP-based PhOLEDs.


Introduction
Nowadays, organic light-emitting diodes (OLEDs) have been considered as a new generation of technologies due to their high brightness, low-cost, low energy consumption, light weight and mechanical flexibility. 1-4 They can be used in many devices such as televisions, mobile phones, and other lighting resources. 5,6 Among different types of OLEDs, the phosphorescent OLEDs (PhOLEDs), serving as the second-generation OLEDs, have recently attracted much attention because they can generate light from both singlet and triplet states and thus allow the internal quantum efficiency to reach nearly 100%. 7,8 In PhO-LEDs, the heavy metal-based phosphorescent emitter is doped in an appropriate host material, which is used to transfer excitation energy to the emitter and confine the exciton in the emissive layer. Amorphous materials have been widely used as host materials in PhOLEDs, because the emitter is uniformly distributed in the amorphous matrix to minimize the concentration quenching effect. 9 In addition, the absence of grain boundaries in amorphous materials further minimizes the exciton recombination. 10 The excited-state properties of the host material play an important role in photoluminescence characteristics of OLEDs. For instance, it has been found that the excited-state molecular aggregation in the vicinity of host material/emitter layer interfaces limits the operation lifetimes of PhOLEDs. 11 Therefore, the fundamental understanding of the excited-state properties of host materials is the key for the design of highly efficient PhOLEDs. In the present work, we explore the excited-state properties of a representative host material using reliable and efficient quantum chemistry methods. Hence, we obtain an accurate description of excited states energies and exciton couplings, where the phase space has been sufficiently sampled. To achieve this, we first benchmark semi-empirical quantum chemical methods, the so-called density functional based tight-binding method (DFTB), as well as its long-range corrected version called LC-DFTB, which allows us to describe excited state properties with reasonable accuracy.
In the past decades, DFTB has been used as an alternative to the widely used density functional theory (DFT) to compute molecular properties. DFTB is 2-3 orders of magnitude faster than DFT and can reach the accuracy of DFT methods by using careful parameterization approaches. 12,13 However, DFTB uses semilocal PBE functionals 14,15 for calculations of the Hamiltonian matrix and thus inherits the errors of typical DFT functionals using generalized gradient approximations (GGA). 16,17 Recently, the LC-DFTB approach was developed by the incorporation of arbitrary LC functionals in the DFTB to partly cure the DFT-GGA errors. 18 The time-dependent scheme of LC-DFTB (TD-LC-DFTB) has been successfully parameterized and benchmarked for a test set of small organic molecules involving charge-transfer excitations, 19 rhodopsins and light-harvesting complexes 20 and biological molecules. 18 In the present study, we will examine the quantitative accuracy of TD-LC-DFTB in calculating the absorption spectra of a most widely used host material in PhOLEDs, 21,22 4,4 0bis(9-carbazolyl)-biphenyl or CBP.
By combining classical MD simulations with quantum mechanical/molecular mechanical methods (QM/MM), where the DFTB methods are used in the QM part, we will be able to address the effects of static and dynamic disorders on the absorption spectra of CBP. The paper will proceed as follows: first, the benchmark calculations of TD-LC-DFTB approach are performed by comparing the excitation energies of CBP molecules using various quantum chemistry methods. Second, the absorption spectra of CBP in the gas and condensed phase are computed using TD-LC-DFTB and different structure sampling methods are discussed. Finally, the aggregation effect on the optical properties of amorphous CBP has been addressed.

Benchmark of TD-LC-DFTB
The bond length alternation (BLA) is defined as the difference in average single and double bond lengths along extended p systems. In the present study, the p-conjugated chain marked with red color in Fig. 1(b) is used to evaluate BLA of CBP. To compare the accuracy of different quantum chemistry methods, geometry optimizations of CBP molecule at the ground state are performed using HF, DFT functionals (B3LYP, PBE, oB97XD) and MP2 with 6-31+g* basis set, respectively. The geometry optimization using DFTB is performed using the third-order DFTB (DFTB3) with 3OB parameter set. 23 In the constrained optimizations, the central, side and improper dihedral angles are fixed at 451, 601 and 01, respectively, which are obtained by fully optimization by oB97XD/6-31+g*, while the other degrees of freedom are relaxed. The excitation energies are computed using CIS, TD-DFT, ADC(2), CC2 and the Bethe-Salpeter equation (BSE) from the GW approximation (GW-BSE) method with 6-31+g* basis set and TD-LC-DFTB with mio-0.8 parameter set, 24 respectively. All DFTB and TD-LC-DFTB calculations in the present study are performed using the DFTB+ program suite, 25 if not specified otherwise. The GW-BSE and ADC (2) calculations are performed using the TURBOMOLE v7.4 computational chemistry package. 26 All other calculations are carried out with the GAUSSIAN 09 package. 27 Recently, the gas phase single-point absorption spectrum of CBP has been studied using GW-BSE method. 28 The GW-BSE method is an accurate and reliable approach for calculating the excited-state properties of extended systems 29,30 and thus it is used to assess the performance of TD-LC-DFTB in this work. Therefore, the geometry is adopted from ref. 28 which is optimized at the ground state using PBE0 functional 31 with the def2-TZVP basis set 32,33 . Afterward, vertical excitation energies and oscillator strengths of the lowest 25 excited singlet states are computed using TD-LC-DFTB with mio-0.8 parameter set. The energy lines are broadened by a Lorentzian function with a width of 0.05 eV and summed up to provide the absorption spectra normalized to the highest peak. Note that the inclusion of the dielectric screening into LC corrected functional can lead to significant shifts in the condensed-phase adsorption spectra. However, we would expect that the inclusion of the dielectric screening makes negligible differences in the adsorption spectra for the systems with strongly localized excitations. The small exitonic couplings between the CBP monomers (shown below) indicates that the localized electronic structure of CBP, which justifies the use of LC corrected functional with bare (i.e., unscreened 1/r) electron-hole interaction.
The self-consistent eigenvalue-only GW method with def2-TZVP basis set is performed to benchmark the accuracy of TD-LC-DFTB for the calculation of the adsorption spectrum of CBP. We use the spectral representations for the response and Green's function, 30,55 which enables an analytic evaluation of energy integrals and derivatives. The advantage of this approach is that it avoids the uncertainties that may arise from the widely used plasmon-pole models and prevents additional computational parameters needed in numerical frequency integration methods.

Molecular dynamics simulations
The force field parameters are taken from the general AMBER force field (GAFF), 34,35 where the atomic charges are generated by the restrained electrostatic potential 36,37 fitting procedure (RESP) calculated at the HF/6-31g* 38,39 level using the Gaussian 09 package. 27 The force field parameterizations are imposed on two types of dihedral angles (see Fig. 1). The first are the proper dihedral angles between two central phenyl groups and the dihedral angle between carbazole moiety and phenyl group (henceforth, referred to as ''central angle'' and ''side angle'' or a and b) and the second type is improper dihedral angles (referred to as g). The dihedral constants for a, b and g are parametrized by fitting to the relaxed potential energy scan along each angle using oB97XD/6-31+g* method. The forcefield parameterization procedures as well as the fitting parameters can be found in the ESI † (S-1). The van der Waals (vdW) interactions are represented by a 6-12 Lennard-Jones potential, where the parameters e and R min,ij are derived from the GAFF force field. Since we take the structures from the MD-trajectories, vdW effects are included. This also holds for the QM/MM simulations, where the vdW interactions are included in the QM/MM Hamiltonian. In gas-phase classical MD simulation, a CBP molecule is first minimized using the steepest descent algorithm and then equilibrated at 300 K for 5 ns with a time step of 1 fs. The Nose-Hoover thermostat 40 is employed to control the temperature. 5000 structures are saved during the equilibrium simulation. For each CBP structure, 25 excitation energies are computed using TD-LC-DFTB and the absorption spectrum is given by the sum of all energy lines broadened by a Lorentzian function with a broadening width of 0.05 eV. In the present study, all classical MD simulations are carried out using the GROMACS (version 2016.3) program package 41 unless stated otherwise.
A disordered thin-film morphology of 5000 CBP-molecules is generated by simulating physical vapor deposition, based on a Monte-Carlo (MC) protocol. During the deposition process, molecules are added to the simulation box (L x = L y = L z = 90 Å) one after another where the flux of molecules is parallel to the z-axis and the xy-plane represents the surface of the substrate. The deposition simulation starts with a molecule at a random initial position, followed by 32 simulated annealing (SA) cycles. Within each SA-cycle, the molecules' conformation space is sampled by 130 000 MC-steps, where only molecule translations or rotations of a dihedral (alpha or beta) are considered. Based on a Metropolis criterion, MC moves are accepted and stored or rejected, where the temperature at MC-step n is given by T n = T s c n . The c = (T e /T s ) 1/N determines the cooling rate as per MC step, T s = 4000 K and T e = 300 K denote the start-and end-temperature, respectively. After the SA-cycles are performed, the conformation with the lowest energy is accepted as the final conformation. Subsequently, the electrostatic and Lennard-Jones potential are updated. This procedure is repeated 5000 times until all molecules are deposited. 42 . After deposition, an equilibrium simulation of 2 ns at 300 K and 1 bar is performed to equilibrate the temperature and pressure. The Nose-Hoover thermostat 40 and Parrinello-Rahman barostat 43 are employed to obtain a reasonable canonical ensemble.
The absorption spectra of CBP in the condensed phase are simulated using two initial sampling methods, ensemble sampling and time-series sampling. In the former approach, the vertical excitation energies of 25 excited states are computed using TD-LC-DFTB for each monomer (5000 in total) in the equilibrated CBP supercell and summed up to generate the absorption spectra. In the time-series sampling, the absorption spectra are simulated using the structures of a specific CBP molecule in 5000 snapshots uniformly selected along 5 ns MD simulation. In both cases, the electrostatic effects from the environment on the excitation energies are considered by including the MM point charges in the DFTB calculation.
In the quantum mechanics/molecular mechanical molecular dynamics (QM/MM MD) simulations, the whole system is divided into a quantum mechanical (QM) region, containing only one CBP molecule treated by DFTB3 with 3OB parameter set, and molecular mechanical (MM) region, where the remaining CBP molecules are treated on a molecular mechanics (MM) level using the modified GAFF. The equilibrated supercell obtained from the classical MD simulation is used as a starting point for a 1 ns QM/ MM simulation with a time step of 1 fs, where 5000 snapshots are taken as reference structures for the excited-state calculations of the CBP molecule in the QM region. The QM/MM MD simulations are performed with the GROMACS program package. 44 The excitation energy calculations including the effects of the electrostatic interactions from the environment are carried out with the DFTB+ program suite. 25 The experimental data 45 are digitized using the open-source tool WebPlotDigitizer. 46 The static and dynamic properties of molecular aggregations are described by the excitonic coupling between the excited states of monomers. Excitonic couplings are conventionally composed of the long-range Coulomb contributions and the short-range (exchange and overlap) contributions. In the present work, only the Coulomb contributions are considered and the complete evaluation of excitonic coupling is planned for the future. The Coulomb coupling can be determined as the interactions between the atomic transition charges of monomers. 47 The Coulomb coupling with the nearest neighbour is considered in the present study. The atomic transition charges are computed using the TD-LC-DFTB method implemented in the DFTB+ program suite. 25 Recently, TD-LC-DFTB excitonic couplings have been successfully benchmarked for rhodopsins and light-harvesting complexes. 20

Results and discussion
Optical properties of CBP have been studied in detail in a recent work, 28 where the performance of DFT approaches has been compared with GW-BSE estimates. As described for many other molecules, standard DFT-GGA and even several hybrid DFT functionals such as B3LYP or PBE0 suffer from the delocalization error, which renders their prediction for many molecules of interest unreliable. Especially for charge-transfer excited states, excitation energies are severely underestimated, which is also the case for CBP. Range-separated DFT functionals, however, can systematically improve the description of these states. In the following, TD-LC-DFTB will be tested and applied for the description of CBP excited states.

Bond length alternation and excitation energies
Excitation energies of conjugated organic molecules are very sensitive to the molecular structure. Dihedral angles and the BLA are two parameters which have a sizable impact on excitation energies. Recently, the importance of these structural details for excitation energies has been discussed in detail. 20 For CBP, different quantum chemical methods predict slight differences in the central geometrical parameters, as can be seen from Table 1. Further information regarding the bond lengths is provided in the Table. S1 (ESI †).
DFTB gives the lowest value of BLA while HF leads to the highest BLA. The values for BLA obtained by hybrid methods lie in between due to the varying moderate amount of exact exchange. 20,48 This finding is in line with an earlier work on other conjugated molecules. 20 The method for optimization of the ground state structure has a notable impact on excitation energies, as can be seen in Fig. 2. The vertical excitation energies computed with various quantum chemistry methods for ground state geometries are provided in the Table S2 (ESI †).
All methods show a slight variation with optimized structures, but differ significantly in absolute excitation energy values. As expected, HF/CIS overestimates the excitation energies for all optimized geometries, while the local and semi-local methods underestimate them due to the (semi-) local nature of the functionals. CC2, ADC(2) and GW-BSE are in excellent agreement with each other. The excitation energies computed by the rangeseparated oB97XD method are slightly higher, as reported before, probably due to the amount of exact exchange used in this functional. LC-DFTB, similar to the LC-DFT methods, tends to overestimate excitation energies systematically, as has been shown for both ground state 19 as well as excited-state 48 optimized geometries, representing absorption and fluorescence energies. Different range-separated functionals have been shown to differ slightly in their predictions, and the excellent absolute agreement of LC-DFTB in this matter can be regarded as fortuitous. For the calculation of absorption spectra, this will simplify the comparison with experimental data. However, more important is the fact that the effect of structural changes on excitation energies is represented accurately, which will be investigated in the following.
As can be seen from Table 1, the different methods predict variations in both, BLA and dihedral angles. As it will be discussed later, the values of the dihedral angles in the bulk are determined by the intermolecular interactions. For the isolated molecule, the rotational barriers (6.5 kcal mol À1 and 3 kcal mol À1 for central and side angles, respectively, see Fig. S1, ESI †) are low and can be overcome at room temperature; however, rotations may be constrained by the surrounding molecules in the bulk which stabilizes the rotamers. On the other hand, the BLA is found to impact excitation energies in conjugated molecules. 20 To estimate the effect of BLA on excitation energies, the CBP molecule in the gas phase is optimized with the side-, central-(a, b) and improper dihedral (g) angles constrained to 451, 601 and 01, respectively. The optimizations are performed by oB97XD/6-31+g*. As shown in Table 2, a clear pattern of increasing BLA with increasing amount of exact exchange is observed in the DFT methods, while the BLA tends to be more precisely predicted with MP2, since HF tends to overestimate this value. 28 The excitation energies follow the BLA trends, i.e. a blue shift with increasing BLA is observed, because the dihedral angles are the same for all methods and the vertical excitation energies only depend on the BLA. The vertical excitation energies can be found in the Fig. S2 and Table S3 (ESI †). A closer inspection of the bond lengths shows that the force field gives the least reliable estimates. In particular, the C-N bond-lengths deviate from the values computed by the other methods, which leads to a smaller excitation energy, as would be expected from the value of the BLA. Therefore, for the final estimation of absorption spectra, the QM/MM geometries using DFTB for the QM region are used. DFTB uses the PBE functional, and its geometrical prediction is therefore quite close to that of the other DFT methods, as discussed in detail in recent works. 20,49 Table 1 Bond length alternation and dihedral angles for the CBP molecule optimized by various quantum chemistry methods. The bond length alternation is defined in Fig. 1

Absorption spectra in gas phase
Recently, a detailed study on CPB electronic structure has been published based on the GW-BSE method, which is taken as a reference here. 28 The lower panel of Fig. 3 shows a comparison of the gas-phase absorption spectra of CBP using TD-DFTB, TD-LC-DFTB and the GW-BSE reference (using the global minimum geometry taken from ref. 28), respectively. As expected, TD-DFTB largely underestimates the excitation energies by about 1 eV, compared to GW-BSE, while TD-LC-DFTB reproduces the GW-BSE absorption spectrum well. The strong red-shift of TD-DFTB spectrum can be related to the well-known notorious self-interaction error of GGA functionals. It is not surprising that the introduction of long-range (LC) functionals in DFTB can reduce the errors of GGA, which leads to a good agreement between TD-LC-DFTB and GW-BSE spectra. The upper panel of Fig. 3 shows the main orbital transitions in the gas-phase absorption spectrum of CBP. Some other electron transitions are shown in Table S4 in the ESI. † The electron density of HOMO is delocalized over the whole molecule, while the lower occupied orbitals (HOMOÀn, n = 2, 3, 5) exhibit a strong electron localization on the carbazole groups. The electron density of LUMO is localized on biphenyl rings, while the electron density of higher unoccupied orbitals (LUMO-n, n = 1, 2, 6) is localized on two carbazole groups. The peak at 308 nm corresponds to the charge transfer (CT) from the carbazole groups to biphenyl rings, which may be sensitive to the rotation of the central dihedral angle. Two higher energy excitations (o250 nm) are related to the locally excited states in the carbazole groups. Therefore, the rotation of central and side dihedral angles may have little impact on high energy excitations, while it affects the lowest excitation energy significantly. 28 To estimate the effect of dihedral angles on excitation energies, we compute the 25 lowest vertical excitation energies and consequently the absorption spectra for different central and side angles, respectively. The dihedral angles are incremented in 101 from 0 to 901 using the ground-state optimized structure at the PBE0/def2-TZVP as the starting structure, taken from ref. 28. For these structures, single-point energy calculations are carried out with TD-LC-DFTB and GW-BSE, respectively. Fig. 4(a) shows that for the central dihedral angles, the lowest-energy peak shifts from 3.91 eV to the higher energy at 4.15 eV and the intensity becomes lower with the increasing central angles, while the location and intensity of the high-energy peak at about 5.2 eV change slightly. Further information regarding the comparison between TD-LC-DFTB and GW-BSE excitation energies of the lowest-energy peaks, as well as the corresponding Kohn-Sham orbitals, is provided in the Tables S5 and S6 in ESI. † It is noted that for the central angle from 0 to 901, the lowest-peak exhibits the excitation transition from charge-transfer excitation to local excitation.
Furthermore, Fig. 4(b) shows a weak, blue-shifted absorption band at the lowest-energy peak in comparison to the central dihedral angle. To quantify the variation found for the different dihedral angles, the standard deviation of the lowest-energy peak has been calculated for central dihedral angle (4.04 AE 0.09 eV), which is about 2 times larger than that for side dihedral angle (4.02 AE 0.05 eV). This finding is consistent with the recent study on the excited-state properties of CBP that the rotational conformation has significant impacts on the low-energy absorption band, while it has minor effects on high-energy bands. 28 Therefore, the CT-type excitation is more sensitive to the central dihedral angle than the side dihedral angle.

Absorption spectra in condensed phase using force field geometries
The absorption in the bulk phase critically depends on the distribution of dihedral angles of the individual molecules, which results from the intermolecular interactions. Different values of the dihedral angles can be stabilized by steric interactions, leading to static disorder in the system. Furthermore, it is expected that dynamic disorder will lead to an additional broadening of the absorption spectrum. In principle, both effects could be included by computing excitation energies Table 2 Bond lengths and bond-length alternation (BLA) for the constrained optimizations with the different methods. Bond length and BLA are given in Å

Bond length [Å]
BLA C 1 QC 2 C 2 -C 3 C 3 QC 4 C 4 -C 5 C 5 QC 6 C 6 -C 7 C 7 QC 8 FF  oB97XD 1.392 1.386 1.396 1.483 1.396 1.386 1.392 0.024  B3LYP 1.398 1.390 1.402 1.484 1.402 1.390 1.397 0.  along the trajectories of all molecules in the simulation box. Since this is computationally demanding, a stepwise procedure is followed, as described below. In order to get a sufficiently large sample of molecules, force field trajectories are used, and absorption energies along each trajectory are computed with TD-LC-DFTB. Then, the trajectories using DFTB/MM are computed in order to obtain spectra based on more accurate geometries.
For one MD-snapshot containing all 5000 molecules in the supercell, the absorption energy of the bulk is then computed by successively evaluating the excitation energy of every molecule in the sample using TD-LC-DFTB, surrounded by the others using electrostatic embedding. The total spectrum is then the overlay of all 5000 single molecule spectra, as shown in Fig. 5(a). It is noted that the absorption spectrum is identical with the spectrum computed without considering the electrostatic interactions from surrounding molecules, which indicates that the electrostatic effects have a negligible influence on the absorption spectrum (see Fig. S4 in the ESI †). As shown in Fig. 5(b), the static disorder is characterized by a significant spread in the distribution of the dihedral angles. The maxima of the aand b-dihedral angles are around 60 and 40 degrees, respectively, consistent with the values of a and b for the structure optimized by oB97XD method in the gas phase. As expected from Fig. 4, planarization of the molecule will lead to a red shift in the lowest absorption band, which is reflected in the absorption spectrum of the ensemble snapshot (red dashed line in Fig. 5(a)). It is red-shifted by 0.09 eV compared to the gas-phase spectrum of the optimized molecule in Fig. 3. 3.3.1 Absorption spectra along trajectories of individual molecules. To investigate the impact of dynamic disorder, three molecules (referred to as M1, M2 and M3) are selected at random positions of the supercell. Afterwards, absorption spectra are computed on a 5 ns classical MD trajectories and excitation energies are evaluated for snapshots every 1 ps, resulting in 5000 conformations in total. As shown in Fig. 5(a), the absorption spectra of M1, M2 and M3 are quite close to the absorption spectrum computed based on the ensemble sampling. Fig. 5   Note that the molecule is not symmetric and therefore a presentation 01-901 is not appropriate. and the fluctuations are constrained to a limited range. It is interesting to see that individual monomers sample only a restricted part of phase space. The results on the ns time scale indicate that the computation of absorption spectra needs to be performed using longer-simulation-time trajectories until convergence. Therefore, using the same starting configurations, the classical MD trajectory is extended to 10 ns and excitation energies are evaluated for snapshots every 2 ps, leading to 5000 conformations in total. The convergence behavior is observed for M2 and M3 using the 10 ns-trajectory, whereas M1 is not converged and its respective angle distribution (mean value = 301) is different from those obtained for M1 in the 5 ns trajectory (mean value = 401) (see Fig. 6) which will lead to a stronger red shift and increased intensities of spectrum. In order to get further insight into the sampling methods, the absorption spectra as well as central angle distribution are examined for snapshots of molecule M1 at each nanosecond (Fig. 7). The importance of the central dihedral angle compared to the side dihedral angle on the CT-type excitation is demonstrated in Section 3.2. As shown in Fig. 7(b), the mean value of the central dihedral angle between 1-5 ns of the trajectory is around 351 whereas between 5-10 ns is around 201.
As shown in Fig. 3, the lowest absorbance peak is located at 4 eV for the optimized structure in the gas phase, while in the condensed phase it is red-shifted by about 0.3 eV (Fig. 5a). In principle, there are several factors determining the absorption maximum, which are (i) the dihedral angles, (ii) the electrostatic environment, (iii) and the exciton couplings. It is shown in ESI † (see Fig. S4) that the electrostatic environment has a minor effect on the absorption spectrum. In the following, the effect of the other two factors will be assessed individually.
3.3.2 Comparison to gas-phase MD spectrum. To estimate the effect of the molecular interactions within the bulk, a MD simulation of CBP in the gas phase is performed, and the excitation energies along this trajectory are computed. Fig. 8(a) shows a comparison of absorption spectra of CBP for 5000

View Article Online
conformations from gas-phase MD and 5000 structures from one condensed-phase MD-snapshot. The lowest-absorbance peak in the condensed-phase spectrum exhibits a red shift of 0.05 eV in comparison with the gas-phase MD spectrum. This is consistent with the central/side dihedral angle distribution shown in Fig. 4, i.e., more conformations with smaller angles occur in the condensed phase than in the gas phase. This indicates that intermolecular interactions in the bulk have a moderate influence on the absorption spectrum of CBP. The difference is not due to the electrostatic effect, but due to the accessible phase space. In the gas-phase MD, the molecule can sample the full angular degrees of freedom, since rotational barriers can be overcome at room temperature. However, in condensed phase molecules may be more restrained, leading to a narrower angular distribution. Since in the bulk, the molecules sample more planar structures, a red-shift of the spectrum is observed (see Fig. 4).

3.3.3.
Effects of rotational conformation. The distribution of dihedral angles for the amorphous ensemble containing 5000 molecules (Fig. 5(b)) is split into 6 categories, where each category represents molecules with dihedral angles in increment of 301. To obtain good statistics, the absorption spectra of each category from 01 to 901 and their complementary categories from 90 to 180 degrees are merged, e.g., structures with dihedral angles of 01-301 and 1501-1801 are pooled together. Fig. 9 shows the partial absorption spectra of dihedral angles for all 3 pooled categories: Category 1: 01-301 and 1501-1801, Category 2: 301-601 and 1201-1501 and Category 3: 601-901 and 901-1201. For the central dihedral angle at the peak with the lowest energy, as the angle gets closer to 01 (i.e., the planarity of biphenyl increases), the spectrum becomes more red-shifted and the intensity increases. Furthermore, it seems there are additional peaks around 4.4 and 5.4 eV for Category 3. For the side dihedral angle, it seems there is no dependency between intensity and dihedral angles. It is also noted that the curve for the side dihedral angle in Category 1 is not smooth due to the low sample size.

Absorption spectra for QM/MM sampling
As discussed above, the structures of a single molecule over a long simulation time are unable to represent all structures in the phase space. The snapshot of the ensemble, on the other hand, includes static disorder but misses the explicit treatment of dynamic disorder, which can be obtained by analysing MD trajectories. To investigate how many molecules are sufficient to represent the ensemble, 3 sets of molecules containing 10, 20 and 50 molecules are randomly selected from different locations of the supercell. A 5 ns long NPT MD simulation is performed. From this trajectory 500, 250 and 100 snapshots are taken for these 3 sets of molecules, respectively, to maintain the same number of structures in every test set. The convergence test with respect to the number of molecules (Fig. S5, ESI †) shows that 50  molecules is sufficient to obtain converged angle distributions as well as a converged absorption spectrum.
As discussed in Section 3.1, DFTB gives a reasonably accurate estimation for the bond lengths in CBP molecules. Here, the QM/MM simulations are performed for amorphous CBP, where 50 molecules randomly distributed in the CBP supercell are treated by DFTB/3OB while the remaining molecules of the system are treated by the modified force field. Fig. 10(a) shows the comparison of the experimental absorption spectrum of CBP 45 with the absorption spectrum using the structures from QM/MM simulations and 5000 conformations from one MD-snapshot, respectively. The QM/MM simulation reproduces the peaks of the experimental absorption spectrum of CBP quite well. The experimental spectrum is much broader than the computed spectrum., which may be due to the defect absorption in the CBP film. A key ingredient of this good agreement for the absorption peaks is the (fortuitous) accuracy of LC-DFTB for the CBP excited states, as discussed above, which allows us to  compute absolute excitation energies in this case. In particular, the BLA is sufficiently well represented by the ground state DFTB/MM calculations. In comparison to the spectrum based on the FF geometries, the lowest-energy peak of the QM/MM absorption spectrum exhibits a red shift of 0.2 eV. The red-shift of excitation energy can be attributed to the smaller BLA using GGA-type DFT method, 50 in comparison to the BLA obtained by the parameterized force field. (see Table 2). In addition, there is a slight deviation in the description of the angular distribution. Fig. 10(b) shows the angle distributions of CBP monomers obtained by MD and QM/ MM simulations, respectively. MD simulation with modified force field generates fewer molecules with angles smaller than 301 (larger than 1501) than DFTB/MM simulation, which indicates a more planarity of CBP conformations obtained by DFTB than force field.
3.3.5 Effects of aggregation. It is well known that the aggregation of chromophores has a significant influence on the absorption spectrum. [51][52][53][54] Delocalization of excited states over two or more molecules could lead to a further red-shift of adsorption spectrum. To examine the aggregation effect on the absorption spectrum of CBP, we compute the nearest-neighbor excitonic coupling, which is an important quantity to describe the extent of molecular aggregation. For all 5000 nearest neighbours, those dimers which are doubly counted (28 dimers), are removed. As shown in Fig. 11(a) The absorption spectra are shown in Fig. 11(b). The low energy part of the monomer absorption spectra is very close to that of the dimer absorption spectra, for all coupling strengths, i.e., there is no indication for a substantial delocalization of excited states. The small couplings of J o 100 meV is in the same range with the dynamical disorder, i.e., the dynamical disorder enforces localization. This is indicated by the fact that the excitation energies computed for monomers and dimers in every class are very similar. If a coupling of the low energy states occurred, an additional red-shift would be expected. Nevertheless, it is very interesting to see that a red-shift occurs which correlates with the coupling strength, although this cannot be induced by electronic delocalization. The results can be interpreted as follows: High excitonic couplings may indicate that dimers are stacked very well. This stacking induces conformations more shifted towards planar conformations, as shown in Fig. 11(c), and more planar conformations lead to a red-shift in the absorption spectrum.

Conclusions
The absorption spectrum of CBP molecules is significantly altered in the amorphous phase due to interactions among the monomers. While electrostatic interactions as well as excitonic couplings turn out to be of minor importance, structural effects, in particular torsion around the main dihedral angles play a major role. Steric interactions induce a planarization of the molecules in the ensemble, which leads to an absorption red-shift of about 0.3 eV w.r.t. a gas phase absorption spectrum.
Conformational sampling is a key to capture the static and dynamic disorder in the system. Due to the packing in the amorphous phase, the individual monomers show only restricted motions, however, relevant conformational changes appear also on the nanosecond timescale, which requires efficient sampling techniques.
Standard force field parameters have to be taken with some care; we found that a major reparametrization of dihedral potentials is required for a realistic model of the amorphous structure. If these important degrees of freedom are modelled too soft or too stiff, the molecular ensembles may be largely misrepresented, since certain parts of the conformational space may be under-or oversampled. This, however, would show up in the absorption spectrum, since excitation energies are very sensitive to these degrees of freedom.
Our simulations could successfully reproduce the experimental absorption spectrum which may indicate that our structural model is quite accurate. We can deduce this from the computed spectrum, since the applied LC-DFTB method very accurately reproduces the excitation energies of high-level quantum chemical methods as GW-BSE, CC2 or ADC (2). Further, the ground state structures have been sampled with DFTB, which shows a very good agreement for the molecular structures, particularly, for the BLA.
The structural model derived here will be used in future work to investigate charge and exciton transfer processes in CBP, which will allow study of the processes behind the degradation of this material.

Conflicts of interest
There are no conflicts to declare.