Getting excited: Challenges in quantum-classical studies of excitons in polymeric systems

A combination of classical molecular dynamics (MM/MD) and quantum chemical calculations based on the density functional theory (DFT) was performed to describe conformational properties of diphenylethyne (DPE), methylated-DPE and poly para phenylene ethynylene (PPE). DFT calculations were employed to improve and develop force field parameters for MM/MD simulations. Many-body Green's functions theory within the GW approximation and the Bethe-Salpeter equation were utilized to describe excited states of the systems. Reliability of the excitation energies based on the MM/MD conformations was examined and compared to the excitation energies from DFT conformations. The results show an overall agreement between the optical excitations based on MM/MD conformations and DFT conformations. This allows for calculation of excitation energies based on MM/MD conformations.


Introduction
Multiscale modelling has become one of the leading themes in modelling materials and their different properties. The most famous use of this term is in the 2013 Nobel Prize in Chemistry ''for the development of multiscale models for complex chemical systems''. Its increasing popularity can also be seen in the titles of published peer reviewed papers: according to the Web of Science, in 2015, 1157 articles had the term ''multiscale'' in their title, a decade earlier the number was 520 and in 1990 it was only 41.
''Multiscale'' is often used synonymously with ''coarse graining''. Coarse graining typically refers to obtaining interaction potentials and parameters for a higher level system from structural equilibrium data. Examples of such are force matching, 1,2 and inverse Boltzmann 3 and inverse Monte Carlo methods. 4,5 The latter two are based on the Henderson theorem 6 which is essentially the Hohenberg-Kohn theorem 7 for classical systems; for a brief derivation and discussion of the relationship between them, see the study of Murtola et al. 8 More heuristic approaches such as the Martini model 9 and the PLUM model 10 are other common approaches; for a comparison between PLUM, Martini and atomistic models, see the study of Bereau et al. 11 One of the leading ideas is that systematic coarse graining allows, at least in principle, also fine graining, that is, remapping the higher level model back to the original more microscopic one. 12 Backmapping procedures exist for the Martini model 13,14 and such approaches have been proven to be useful in modelling lipids and proteins, see e.g. the study of Pannuzzo et al. 15 One somewhat less considered but an important issue is sampling at different levels of coarse-graining. 16 For more details on coarse graining, see recent reviews in ref. 8, 17-19. Multiscale modelling is a much broader concept. For example, instead of linking scales via deriving new interaction potentials, in hybrid simulations some part of the system is described with a different resolution from the rest and information is transmitted between the two different regions. Examples of such are QM/MM (quantum-molecular mechanics), MM/CG (molecular mechanics-coarse grained) and even QM/QM (quantum-quantum). The idea is that more detailed information in some well defined region is sought after and the crucial issue is how to couple the main system and the subsystem. This has been discussed extensively, see, e.g., ref. 20, but the essence is that both the dynamic and static properties must be communicated between the systems. This may include polarization, changes in the charge-state of the system, and so on. Yet another multiscale approach is the so-called adaptive resolution method, or AdResS. In this case, run-time information is transmitted between layers of description ranging from atomistic even up to continuum. 21,22 Electronic excitations pose a significant challenge for multiscaling 23-31 since typical DFT methods describe the ground state. An assessment of the interplay between the molecular electronic structure, the morphological order, and thermodynamic properties requires the knowledge of the material morphology at atomic resolution, as well as strategies to couple quantum mechanical techniques to classical environments for accurate evaluation of electronic excitations. [32][33][34] Morphology can have several characteristic length scales and be kinetically arrested. Employing all-atom molecular dynamics is limited to a few microseconds, which might be too short to fully relax molecular positions and orientations during aggregation 35 in polymer melts, [36][37][38] or polymers in miscible solutions. 39 In such cases more coarse representations might be helpful to overcome the limitations of atomistic models. [40][41][42][43] Using empirical atomistic potentials in multiscale simulations of excitations based on quantum calculations requires that the structural description at different levels of resolution are compatible with each other. For example, bond length deviations or fluctuations in angles and torsions can lead to substantial artifacts if the backmapped/ fine-grained geometries do not match the potential energy surfaces (PESs) of the underlying quantum mechanical system. Such a situation regularly arises for conjugated polymers since conjugation can depend sensitively on conformation. In (semi)flexible polymers, conjugation along a single chain can be broken due to large out-of-plane torsions between two repeat units. Broken conjugation and wave function localization [44][45][46] are often intuitively interpreted based on a simple empirical criterion, the dihedral angle between two adjacent repeat units. 28 In general, details are specific to the backbone chemistry, functionalization by side chains, and solute-solvent interactions. Characteristics of conjugation also directly influence the localization behavior of electronic excitations and hence the electronic and optical properties of the polymer.
In this paper, some of the underlying challenges pertaining to the transfer of structural atomistic details between quantum and all-atom resolutions are demonstrated. As an example, we consider the calculation of the optical properties of poly-PPE (see the chemical structure in Fig. 1), a relatively rigid conjugated polymer consisting of aromatic phenyl rings bridged by alternating single and triple carbon bonds. PPEs can be prepared in a variety of morphologies, ranging from extended single chains to polydots and their optical properties make them particularly attractive for use in fluorescence imaging and sensing. [47][48][49] Due to the importance of backbone torsions for conjugation and hence excitations, we compare PESs for phenyl rotations in diphenylethyne (DPE, see Fig. 2) obtained using density functional theory (DFT) to the ones from all-atom simulations using a standard force field and experimental data. Significant discrepancies were found and as a result, the atomistic force field was re-parameterized. With this modified force field, ground state geometries are optimized for n-PPE oligomers with n = 1,. . ., 10 and then used in GW-BSE calculations. The associated excitation energies are benchmarked with results from quantum-mechanical treatment, revealing qualitatively similar characteristics as a function of n but deviations at the quantitative level. Finally, conformations from MD simulations of 2,5-dinonyl-10-PPE solvated in toluene are used in a QM/MM setup to evaluate absorption and emission spectra.

Methodology
MM/MD calculations were performed by a force field of OPLS (optimized potential for liquid simulation) 50-52 form using GRO-MACS simulation software version 4. 53 The force field parameters are taken from the polymer consistent force field 54,55 (PCFF) as converted to the OPLS form in ref. 56 and 57. We refer to it as PCFF* from now on. The OPLS potential energy function consists of harmonic bond stretching (V bond ), angle bending potential (V angle ), and non-bonded terms (V non-bonded ) including Lennard-Jones (LJ) and electrostatics, and proper and improper dihedral potential terms (V torsion ): 50-52 The parameters k b,i and k y,i are the bond force constant for bond i and the angle force constant for angle i, respectively. r 0 and y 0 are the initial (reference, equilibrium) bond distance and angle bending, respectively. k 1,i , k 2,i ,. . . are the torsional force constants for each dihedral i. q i e is the partial atomic charge of atom i in which e is the charge of one electron, s ij are the LJ radii, e ij are the LJ energies (well-depth) and r ij are the distances between atoms i and j. The geometric combination rules were used following the convention adopted in the OPLS force field s ij ¼ s ii s jj À Á 1 2 and e ij ¼ e ii e jj À Á 1 2 ! . The intramolecular non-bonded interactions were evaluated for atom  pairs separated by three or more bonds. The 1,4-intramolecular interactions were reduced 50-52 by a factor of 1/2. To obtain relaxed scans of potential energy surfaces (PESs) from MM/MD, energy minimization of the DPE molecule in a vacuum was performed using the conjugate gradient method followed by a short MD run (100 ps) with constant particle number (N) and temperature (T). A Langevin thermostat 58 with a time step of 1 fs and open boundary conditions was applied. Temperature was maintained at 10 K with 10 fs damping constant. All LJ interactions were cut-off at 1.2 nm. A plain cut-off scheme was used for electrostatic interactions with 2.0 nm real space cut-off: with open boundary conditions plain cut-off can be used. For systems with periodic boundary conditions, the particle-mesh Ewald (PME) 59,60 should be used instead. For more discussion about the importance of electrostatic interactions, please see ref. 61. The cut-off distance for the short-range neighbor list was 1.2 nm and the neighbor lists were updated at every step. The intention was to evaluate the ground state energies of DPE molecules with different torsional angles between aromatic rings. To do that, after the first energy minimization step, a short MD run at very low temperature was used to bring the system out of possible local minima. Then a second conjugate gradient energy minimization was performed to obtain the ground state MM/MD PES.
DFT optimizations and relaxed PES scans were performed using the B3LYP exchange correlation functional 62-65 and the def2-TZVP basis set 66 as implemented in the Orca package. 67 Due to the lack of van der Waals (dispersion) interactions in standard DFT, Grimme's DFT-D3 method 68 was employed.
In order to calculate electronically excited states, many-body Green's function theory in the GW approximation with the Bethe-Salpeter equation (GW-BSE) 69 was employed, since static DFT 70 cannot describe coupled electron-hole excitations. For details of the application to molecular systems, the reader is referred to ref. [71][72][73][74][75][76]. The GW-BSE method is based on a set of Green's function equations of motion which contain electronhole interaction (BSE) leading to the formation of excitons. It utilizes the DFT molecular orbitals and energies to calculate the one-particle Green's function (G) and screened Coulomb interaction (W) to obtain single-particle excitations within the GW approximation as introduced by Hedin and Lundqvist. 69 An electron-hole excitation cannot be described in an effective single-particle picture but instead requires explicit treatment of a coupled two-particle system. The electron-hole amplitudes and associated transition energies can be obtained by solving the Bethe-Salpeter equation [72][73][74] . For the calculation of excitation energies according to the GW-BSE method, first DFT calculations were performed using the Orca package, 67 the B3LYP functional, 62-65 the effective core potentials of the Stuttgart/ Dresden type, 77 and the associated basis sets that are augmented by additional polarization functions 78 of d symmetry. The specialized GW-BSE implementation for isolated systems [71][72][73]79 available in the XTP module of the VOTCA software package 80,81 is used in all further steps related to the excitations. For molecular visualizations, Visual Molecular Dynamics (VMD) 82 and Jmol 83 were used.

Force field parametrization
Due to the influence of conformational details on the optical properties 84,85 in PPEs, one needs to determine if the force field yields reliable minimum energy configurations. Hence, relaxed scans of the potential energy surface (PES) were obtained using both MM/MD and DFT.
The resulting PESs are shown in Fig. 3. The PCFF* result (red triangles) shows a minimum at 901, corresponding to twisted phenylene rings. In contrast, the result of the DFT-based scan (black squares) indicates a minimum energy configuration in which the two phenyl rings are co-planar, which is also extracted from experiments. 48,86 The force field predicts a practically free rotation of phenylenes for T Z 0, while a barrier of around B4 kJ mol À1 (B1.5 k B T) is obtained using DFT. The latter is comparable to the one reported in ref. 87. The experimental potential barrier is around B2 kJ mol À1 . 48,86,88 Overall, the scans imply that the PCFF* force field does not correctly model the ground state conformations of DPE, which can have severe implications for the derived optical properties.
To remedy this situation, the existing force field is refined by the addition of a torsional potential between the two adjacent phenylenes (see Fig. 4 for definition of involved atoms). By fitting eqn (4) to the differences of DFT and PCFF* potential energy surfaces, corresponding Ryckaert-Bellemans 89 force parameters, provided in Table 1, were obtained. The PES is re-calculated with MM/MD using the modified force field, yielding the scan as shown in Fig. 3

(blue circles). It is in good agreement with the DFT result.
To assess the transferability of the modified force field, we repeat the above scans of the torsional potential for para methylated-DPE (see the chemical structure shown in Fig. 2). The PESs resulting from both MD and DFT calculations are shown in Fig. 5. With the modified force field (blue circles) one can observe a good agreement with the DFT data (black squares). Both approaches predict a minimum energy configuration at 1801 twist. The energetic preference of this cis conformation of Me-DPE over the trans conformation (01) is driven by attractive dispersion interaction among the two CH 3 . While this preference is also obtained with the original PCFF* force field (red triangles), no barrier between cis and trans configurations is found. In terms of obtaining minimum energy configurations and energy barriers in the PES, the modified PCFF is clearly more reliable.

Optical excitations in single molecules
For systems such as solvated polymer chains, the system size makes use of classical simulations inevitable to obtain structural information. Even with the modified force field at hand, it is not automatically guaranteed that the use of the MM/MD geometries in QM/MM schemes does not lead to spurious errors in the computed excitations. To further assess the level of reliability of such calculations, the evolution of the optical absorption properties of Me-DPE is examined as a function of phenyl torsions based on the respective optimized geometries.
The optical absorption spectra resulting from the GW-BSE are shown in Fig. 6. DFT optimized geometries were used in (a) and MD energy minimized geometries using the modified force field were used in (b). The height of the curves indicates the strength of the excitation. Comparing both results, it is evident that the same dependency on the torsional angle is obtained by both approaches. With increasing twist from 01 to 901 the main absorption peak gradually shifts to higher energies while its strength decreases at the same time until it vanishes at 901. Inspection of the electron and hole densities of the excitations for co-planar and perpendicular (see bottom of Fig. 6) also reveals no localization of the excitation during the rotation confirming that the conjugation via the CRC bond is indeed strong. The identical behavior of the lowest energy excitations (which are typically those of interest) for both MM/MD and DFT conformations indicates that the modified force field is suitable for use in QM/MM calculations.
So far the analysis has been limited to small model systems. As a next step towards more realistic system sizes, para phenylene ethynylene (PPE) oligomers (see Fig. 1) are investigated. Increasing the number of repeat units n from one to 10, lowest optically active excitation energies were computed with the GW-BSE for geometries optimized using quantum (GW-BSE@QM) and force field (GW-BSE@MM) approaches, respectively. The results shown in Fig. 7 exhibit a monotonous decrease with n for both approaches. Such a strong size-dependence can be traced back to an increase in the size of the conjugated system. From the particle-in-a-box model, one can estimate, e.g., the optical excitation energy of an infinitely long chain via O(n) = O N À a/n. By fitting the data for n 4 3 to this model, a value of O QM N = 3.08 eV is obtained for QM geometries. For n 4 7 the respective excitation energies vary only slightly and approach the region in which experimental absorption is measured in   N is a cumulative result of slight discrepancies in the bond length within the phenylenes and the C-C bridge bonds. In conclusion, the use of geometries determined using MM/MD in GW-BSE calculations can be expected to lead to slight quantitative overestimates of excitation energies. Qualitatively, however, a satisfying agreement is found.

Absorption and emission of solvated 2,5-dinonyl-10-PPE
As an illustrative example of a typical application of a combined quantum-classical simulation of optical excitations, a single chain of 10-PPE was functionalized by nonyl side chains in 2,5 positions of the phenyl rings and solvated by toluene. The modified PCFF force field for the PPE backbone is used in combination with OPLS for nonyl and toluene. For technical details of the simulations, see ref. 90. A sample configuration of the backbone and side chains is taken from the 56 ns long trajectory and is shown in Fig. 8(a). Since toluene is a poor solvent for the nonyl, one can observe extended and partially strongly interacting side chains. As a consequence, the backbone is under considerable non-uniform stress leading to the overall curvature of the usually rigid polymer.
Starting from the same initial configuration, sets of 11 snapshots with time steps of 10 fs, 100 fs, 1 ps, and 10 ps, respectively, have been taken from the classical MD trajectory and optical excitation energies were calculated for them. For snapshots with a time step of 10 fs and 100 fs, no variation in the excitation energies was observed. We noticed that the excitation energies start to change once there is a 1 ps time step between the snapshots. A plot for excitation energies of these 11 snapshots for different time steps is given in the ESI. † These 11 snapshots with a time step of 1 ps are taken from the classical MD trajectory and each of the snapshots is partitioned into a quantum (the backbone) and a classical region comprising the side chains and the solvent molecules. QM and MM regions interact via static partial charge distributions. The aim of this setup is to evaluate the excitations of the polymer backbone taking its curved conformation into account while reducing discrepancies between the force-field and QM geometries, as much as possible. At the same time, the bridging carbon-carbon bond between the phenyl 2 and 5 positions and the nonyl side chain, defining the boundary between QM and MM regions, needs to be broken and the dangling bond is saturated by the hydrogen atom. This can be achieved with the help of a re-mapping scheme based on the definition of molecular fragments. Using centers of mass and gyration tensors, fragments of optimized QM configurations were mapped onto the orientation and alignment of the corresponding fragments in the MD configurations. Fig. 8(b) illustrates the re-mapping scheme for PPE. Each phenyl ring (PHE), ethyne pair (ETH), and terminal methyl group (CH3) are defined as a unique fragment. A 10-PPE backbone is hence subdivided into a total of 23 fragments (10 PHE, 11 ETH, and 2 CH3) for mapping purposes.
With the re-mapped conformations at hand, the coupled GW-BSE/MM system is solved and the absorption spectrum is determined as an average over the eleven snapshots. Individual spectra are broadened by Gaussian functions with a FWHM of 0.3 eV and the resulting spectrum is shown as a blue line in Fig. 9. It is characterized by a single peak at an energy of 3.64 eV, which is larger than the value of 3.11 eV obtained for an isolated single oligomer. This spectral blue shift is a direct result of the polymer's curvature constrained by the side chain interactions. With the re-mapping scheme it is also possible to approximate emission spectra by using excited state QM geometries as a reference. 91 Upon excitation, electrons are promoted to higher, often anti-bonding, molecular orbitals causing an extension of bonds. Constrained by side chains, a more general modification of the overall conformation (at least on short time scales) cannot be expected. Solving the GW-BSE/MM system based on excited  state re-mapping yields the emission spectrum shown as a red line in Fig. 9. While no changes in the spectral shape can be observed, the peak position of the emission is red-shifted by 0.29 eV compared to the absorption peak. This Stokes shift is in good agreement with experimental data in the range of 0.3-0.4 eV. 48

Conclusions
A combination of atomistic (MM/MD) and DFT calculations was performed to describe the conformational properties of diphenylethyne (DPE), methylated-DPE and poly para phenylene ethynylene (PPE). MM/MD simulations based on the PCCF* force field were not able to provide a good description of the ground state conformation of the DPE molecule. Due to this, DFT calculations were employed to develop force field parameters to improve the MM/MD simulations. The modified force field was able to describe the conformation of methylated-DPE in agreement with DFT results. The GW-BSE method was utilized to describe excited states of the methylated-DPE and n-PPE polymer with n = 1, 2,. . ., 10. Optical excitations were obtained for the methylated-DPE and nPPE based on MM/MD energy minimized structures using the modified force field and DFT optimized geometries. The results for methylated-DPE show that the lowest energy excitations based on the MM/MD conformations and DFT optimized geometries follow the same pattern. This nearly identical behavior for the lowest energy excitations indicates that one can describe optical excitations using the GW-BSE method based on MM/MD conformations. The results for the excitation energies for nPPE indicate that there is an overall agreement between the results of the GW-BSE based on MM/MD energy minimized structures and DFT optimized geometries. There is a discrepancy of around 0.25 eV between the two. This discrepancy is a cumulative result of geometric differences between MM/MD and DFT structures. Overall agreement between MM/MD and QM based excitations is enough to validate the use of MM/MD conformations as the basis for the calculation of optical excitations using the GW-BSE method. Fig. 9 Optical absorption (blue) and emission (blue) spectra of 2,5-dinonyl-10-PPE in toluene as obtained from GW-BSE/MM calculations using fragment based re-mapping. Both spectra are averages over 11 snapshots taken every 1 ps, respectively. A Stokes shift (red shift) of 0.29 eV is found, which compares well to experimental observations of 0.3-0.4 eV. 48