Solvent effects on optical excitations of poly para phenylene ethynylene studied by QM/MM simulations based on Many-Body Green's Functions Theory

Electronic excitations in dilute solutions of poly para phenylene ethynylene (poly-PPE) are studied using a QM/MM approach combining many-body Green's functions theory within the $GW$ approximation and the Bethe-Salpeter equation with polarizable force field models. Oligomers up to a length of 7.5\,nm (10 repeat units) functionalized with nonyl side chains are solvated in toluene and water, respectively. After equilibration using atomistic molecular dynamics (MD), the system is partitioned into a quantum region (backbone) embedded into a classical (side chains and solvent) environment. Optical absorption properties are calculated solving the coupled QM/MM system self-consistently and special attention is paid to the effects of solvents. The model allows to differentiate the influence of oligomer conformation induced by the solvation from electronic effects related to local electric fields and polarization. It is found that the electronic environment contributions are negligible compared to the conformational dynamics of the conjugated PPE. An analysis of the electron-hole wave function reveals a sensitivity of energy and localization characteristics of the excited states to bends in the global conformation of the oligomer rather than to the relative of phenyl rings along the backbone.


Introduction
Conjugated polymers have attracted a lot of interest due to their variable functional properties and application potential, e.g., in biochemical sensors [1,2,3,4,5,6], lasers [7], light emitting diodes (LEDs) [8,9], organic transistors [10], and photovoltaic devices [11,12]. The ease of processability and band gap tunability of polymeric semiconductors facilitates the realization of this potential, since it provides the opportunity for a targeted manipulation of electronic and morphological properties of single polymer chains and their aggregates. This, in turn, can be achieved by synthetic strategies, exploitation of properties of functional side chains, and/or solvent-induced transitions. One specific example in which both morphological and electro-optical properties of polymers are purposefully modified is the formation of highly fluorescent conjugated polymer dots for fluorescence imaging in live cells [13]. Less toxicity together with flexibility and biocompatibility make these polydots attractive substitutes for their inorganic counterparts [14,15].
Among fluorescent polymers, poly para phenylene ethynylenes (poly-PPE) (see the chemical structure in Figure 1) are a class of strongly conjugated polymers with a rigid backbone and absorption and emission of light tunable from the ultraviolet (absorption) to the visible (emission) range [14]. In particular, PPEs can be used as fluorescence sensors because their fluorescence intensity is sensitive to the presence of other co-solutes. Due to the significance of polymer conformations for the photophysical properties of functionalized PPEs, it is desirable to steer the tendency of the polymer to stay in single strands or to aggregate in a particular solvent by rational manipulations [16]. It has been demonstrated that modification of side chain functionalization by, e.g., substituting surfactants with varying degree of hydrophobicity and solute-solvent or solvent-air interactions, can lead to controlled morphologies and, concomitantly, optical properties. In flexible and semi-flexible polymers conjugation along a single chain may be broken [17,18,19] due to a substantial out-of-plane torsion angle between two adjacent repeat units. However, such a simple criterion may be insufficient to describe and interpret the complex interplay between the actual chemistry of the backbone, functionalization by side chains, and solute-solvent interactions and their cumulative effect on the localization characteristics of excitations and hence the electronic and optical properties of the polymer.
In this study, a combined quantum-classical (QM/MM) approach is used to investigate optical properties of 2,5-dinonyl poly para phenylene ethynylene oligomers in dilute solutions with toluene and water, respectively, as model systems. Conformations of the PPE chains are explored by using MD simulations. Electronic excitations are calculated based on many Body Green's Functions theory within the GW approximation and the Bethe-Salpeter equation (GW -BSE). The use of the latter technique in traditional quantum-chemical applications has recently increased [20,21,22,23], not least due to its accurate prediction of both localized Frenkel and bimolecular charge transfer excitations [24,25,26]. Linking GW -BSE to a classical environment, represented at atomistic resolution by a polarizable force field, allows for the determination of optical properties in realistic environments from the self-consistent solution of the coupled QM/MM system. With this approach, it is possible to disentangle the conformational (as a result of side chain-solvent interactions) and electronic (due to local electric fields and polarization effects) contributions to the absorption spectra.
The rest of this paper is organized as follows: The methodologies and computational details, are described in sec. 2. The resulting optical properties for isolated oligomers with up to ten repeat units and their sensitivity to structural details are presented in sec. 3.1 and structural properties of dinonyl-10-PPE in solution are discussed in sec. 3.2. In sec. 3.3, the respective optical absorption spectra resulting from GW -BSE/MM calculations are discussed with a focus on the electronic contributions of the environment and the conformational dynamics. A brief summary in sec. 4 concludes the paper.

Methodology
Classical molecular dynamics (MM/MD) simulations were performed using an OPLS type (optimized potentials for liquid simulations) [27,28,29] force field. The parameters were taken from Refs. [30,31] in which PCFF (polymer consistent force field) force field parameters were converted to OPLS form. Modification for torsional potential parameters of the phenylene rings as in Ref. [32] were employed. To study the behavior of PPE polymers in explicit solvents, water molecules were described using the SPC/E [33] model for water and the OPLS force field for toluene. Geometric mixing rules [σ ij = (σ ii σ jj ) 1 2 and ǫ ij = (ǫ ii ǫ jj ) 1 2 ] for Lennard-Jones (LJ) diameters (σ) and LJ energies (ǫ) were used for atoms of different species according to the OPLS conventions [27,28,29]. The reader is referred to Ref. [34] for more in-depth discussion on mixing rules. Non-bonded interactions between atom pairs within a molecule separated by one or two bonds were excluded. Interaction was reduced by a factor of 1/2 for atoms separated by three bonds and more. Simulations were run using GROMACS version 5 [35]. A 1.2 nm cutoff was employed for the real space part of electrostatics and Lennard-Jones interactions. The long-range electrostatics was calculated using particle-mesh Ewald (PME) [36,37] with the reciprocal-space interactions evaluated on a 0.16 grid with cubic interpolation of order 4. The importance of proper treatment of electrostatics in MM/MD simulations is discussed in detail in Ref. [38]. The velocity-Verlet algorithm [39] was employed to integrate the equations of motions with 1 fs time step. A Langevin thermostat [40] with a 100 fs damping was used to keep the temperature of the system at 300 K. The systems were energy minimized using the steepest descents algorithm. 100 ps simulations in constant particle number, volume and temperature (NVT) ensemble at 300 K were performed on the energy minimized systems. The simulation box size was (15 × 13 × 13) nm 3 for dinonyl-10-PPE. Simulations were continued in constant particle number, pressure and temperature (NPT) ensemble at 300 K and 1 bar controlled by Parrinello-Rahman [41] barostat with a coupling time constant of 2.0 ps. Molecular visualizations were done using Visual Molecular Dynamics (VMD) software [42].
The excited state electronic structure is evaluated using many-body Green's functions theory within the GW approximation and the Bethe-Salpeter equation (BSE). Green's functions equations of motion are the basis for this approach which contains both the nonlocal energy-dependent electronic self-energy Σ and the electron-hole interaction leading to the formation of excitons, described by the BSE. The first of three steps of the procedure is the determination of molecular orbitals and energies on the level of density-functional theory (DFT) by solving the Kohn-Sham equations Here, V PP is a pseudo-potential (or effective-core potential), V H the Hartree potential, and V xc the exchange-correlation potential. Single-particle excitations are then obtained within the GW approximation of many-body Green's functions theory, as introduced by Hedin and Lundqvist [43], by substitution of the energy-dependent self-energy operator Σ(r, r ′ , E) for the DFT exchange-correlation potential, giving rise to the quasi-particle equations (2) The self-energy operator is evaluated as where is the one-body Green's function in quasiparticle (QP) approximation and W = ǫ −1 v is the dynamically screened Coulomb interaction, comprising the dielectric function ǫ, computed within the random-phase approximation, and the bare Coulomb interaction v. The ground state Kohn-Sham wave functions and energies are used to determine both G and W . Since DFT underestimates the fundamental HOMO-LUMO gap, the self-energy and the resulting QP energies may deviate from self-consistent results. To avoid such deviations, we employ an iterative procedure in which W is calculated from a scissor-shifted Kohn-Sham spectrum. From the resulting QP gap, a new value for the scissor-shift is determined and this procedure is repeated until convergence is reached. For each step, the QP energy levels are iterated and the Green's function of eq. 4 and thus the self-energy are updated. A one-shot G 0 W 0 calculation from Kohn-Sham energies may differ from iterated results by up to several 0.1 eV. Note that this (limited) self-consistency treatment does not change the QP structure of eq. 4 (due to satellite structures or other consequences of a self-consistent spectral shape of G(ω)). An electron-hole excitation, e.g., resulting from photoexcitation, can not be described in an effective single-particle picture but requires explicit treatment of a coupled two-particle system. Within the Tamm-Dancoff approximation (TDA) [44], the electron-hole wavefunction is given by Φ S (r e , r h ) = occ α virt β A S αβ ψ β (r e )ψ * α (r h ) where α and β denote the single-particle occupied and virtual orbitals, respectively, and A αβ represent the electron-hole amplitudes. These amplitudes and the associated transition energies Ω S can be obtained by solving the Bethe-Salpeter equation in which K eh = ηK x + K d (η = 2 for singlets, η = 0 for triplets) is the electron-hole interaction kernel comprised of bare exchange (K x ) and screened direct terms (K d ), respectively. For practical calculations according to the GW -BSE method, first single-point Kohn-Sham calculations are performed using ORCA [45], the B3LYP functional [46,47,48,49], effective core potentials of the Stuttgart/Dresden type [50], and the associated basis sets that are augmented by additional polarization functions [51] of d symmetry [52]. All steps involving the actual GW -BSE calculations are performed using the implementation for isolated systems [53,54,25,21], available in the VOTCA software Table 1. Electronic structure data for n-PPE oligomers with n = 1, . . . , 10 based on QM and MM optimized geometries: HOMO-LUMO gap from Kohn-Sham (E KS g ), quasi-particle (E QP g ) energies, optical excitation energy (Ω), and the contributions to it from free interlevel transitions ( D ) and electron-hole interaction ( K eh = K d + 2K x ). All energies in eV. package [55,56]. In VOTCA, the quantities in the GW self-energy operator (dielectric matrix, exchange and correlation terms) and the electron-hole interaction in the BSE are expressed in terms of auxiliary atom-centered Gaussian basis functions. We include orbitals of s, p, d symmetry with the decay constants α (in a.u.) 0.25, 0.90, 3.0 for C and 0.4 and 1.5 for H atoms, yielding converged excitation energies. It was also confirmed that the addition of diffuse functions with decay constants smaller than 0.06 a.u. to the wave function basis set does not affect the low-lying excitations. For all systems considered in this paper, polarizability is calculated using the full manifold of occupied and virtual states in the random-phase approximation. Quasiparticle corrections are calculated for the 2n occ lowest-energy states, and n occ occupied and n occ virtual states are considered in the Bethe-Salpeter equation. Further technical details can be found in Refs [53,54,25].

Optical absorption energies of isolated oligomers
Conformations of solvated PPE oligomers as obtained from classical MD simulations form the basis for the determination of optical excitations in a mixed QM/MM setup. The underlying assumption is that the details of molecular geometries resulting from the use of a classical force field are consistent with the chosen quantum mechanical description of the ground state. To confirm this the geometries of single n-PPE oligomers with n = 1, . . . , 10 were optimized in vacuum using both DFT (def2-TZVP [57] basis set and B3LYP functional with Grimme's D3 dispersion corrections [58]) and MD on the basis of the modified PCFF force field. Both approaches yield planar configurations of the PPE backbone for all values of n.
While qualitatively identical, quantitative differences can be observed. Most notably, the bonds in the phenyl rings and the C − C bond connecting the ring to the ethyne are elongated by about 2% in MD compared to DFT. In contrast, the length of the C ≡ C triple bond results 1.21Å in both cases. In the next step, GW -BSE calculations were performed to gauge the effect of these differences on the electronic and optical properties of the oligomers. The resulting electronic structure data, summarized in Tab. 1, illustrates the effects of GW -BSE with respect to the calculation of excitation energies: Taking the quasi-particle corrections to the Kohn-Sham energies into account increases the HOMO-LUMO gap for, e.g., the QM optimized 10-PPE from 2.88 eV to 5.53 eV, which reflects the well-known underestimation of the fundamental gap by DFT. The energy for the lowest optically active coupled electron-hole excitation gives 3.11 eV. Due to the fact that the excitation is not a pure HOMO-LUMO transition but has additional contributions from lower occupied and higher virtual single-particle orbitals, the contribution of the independent transitions, D , is with 6.01 eV slightly larger than E QP g . The associated effectively attractive electronhole interaction, K eh = K d + 2K x , in this structure amounts to 2.90 eV. The obtained excitation energy is in good agreement with the experimental values of 3.0-3.2 eV obtained from absorption peaks of dilute solutions of PPE in good solvents. As a function of the number of repeat units n, a monotonous decrease of all energies is found. This quantum-size effect is anticipated for strongly conjugated systems such as PPEs. From the particle-in-a-box model one can estimate, e.g., the optical excitation energy of an infinitely long chain via Ω(n) = Ω ∞ − a/n. By fitting the data for n > 3 to this model, one obtains a value of Ω ∞ = 3.08 eV indicating that for further studies of solvent effects on the excitations in PPE, it is reasonable to consider 10-PPE in the QM/MM setup.
Qualitatively identical results were obtained for the GW -BSE calculations based on the MM optimized geometries. However, there are some noticeable quantitative differences, most importantly with respect to the excitation energy Ω which is consistently larger as compared to the QM structures. For large n, the deviation amounts to 0.25 eV and is a cumulative result of the geometric differences discussed above. Overall, the agreement is satisfying enough to conclude that the use of MD simulated conformations of 10-PPE in the following is capturing the relevant physico-chemical details.

Structural properties of solvated 2,5-dinonyl-10-PPE
Conformations of 2,5-dinonyl-10-PPE were studied in explicit water and toluene. Water is a poor solvent for both the backbone and the side chains while toluene is a good solvent for the backbone and a poor solvent for the side chains [30,31]. Figure 2 shows the structure of 10-PPE (a) in toluene and (b) in water at 7.7 ns in the NpT ensemble. For clarity, water and toluene molecules are not shown. In toluene, the backbone remains extended and the side chains are dispersed and separated from each other as well as the backbone. This is in agreement with the results of Ref. [30]. Structural studies using small angle neutron scattering (SANS) have shown that dialkyl PPE forms a molecular solution with an extended backbone at high temperature and low concentrations [30,59]. In water ( Fig. 2(b)), the side chains start to aggregate toward each other and the backbone. This is in agreement with Ref. [30,31]. Another important parameter is the correlation of aromatic rings along the backbone of PPE. The interplay between the arrangements of aromatic rings in PPE polymers and their electro-optical properties has been studied by several groups (see e.g. [60,61,62]). The orientational order parameter [63], given by is a measure to quantify how aromatic rings within PPE polymer backbone are correlated. θ is the angle between the normal vectors to the planes of two aromatic rings which are apart from each other by a distance ∆n. P θ describes the average alignment of aromatic rings. Since each vector normal to the planes of the aromatic rings can be considered as a reference direction to calculate θ and then P θ , one needs to consider the vector normal to each plane as a reference direction and take an average: there are two averages in the calculation of P θ , one of which is the time average over time frames and the other one is over the selection of a vector normal to the plane of the rings. P θ can have values [− 1 2 ,1] [63]. P θ > 0 describes a co-planar alignment of aromatic rings, while P θ < 0 indicates perpendicular alignments. P θ = 0 and P θ = 1 refer to completely random and fully co-planar alignment of the rings, respectively [30,31]. Figure 3 shows the order parameter versus ∆n for 10-PPE with nonyl side chains in toluene and water after 7.7 ns MD simulations in the NpT ensemble. The time average is taken over the frames of the last 1 ns (0.6 ns) of the MM/MD trajectory with 100 ps time step between the frames for 10-PPE in toluene (water). Having a value of around 0.4 (0.2) in toluene (water) indicates a correlation between the aromatic rings. This refers to an angle of around 39 • and 47 • for 10-PPE in toluene and water, respectively. In Ref. [60], the authors discussed optical properties of dialkyl and dialkoxy-PPEs in chloroform and dichloromethane, and the average angle of aromatic rings using a configuration-coordinate model in which they concluded the angle to be around 40 degrees.

Optical absorption of solvated 2,5-dinonyl-10-PPE
Excitation energies in complex molecular environments, such as molecular aggregates or solute-solvent mixtures, can accurately be obtained by treating the subpart of interest quantum-mechanically and embedding it into an environment at molecular mechanics resolution [64,65,66,67]. Here, we realize such a QM/MM scheme based on GW -BSE by representing the molecules in the MM region by a set of atomic properties such as static partial charges and polarizabilities, which then interact among each other and the GW -BSE part via classical electrostatic potentials. Polarizable interactions are taken into account via Thole's model [68,69]. For both neutral and excited complexes, total energies of the combined QM/MM system are obtained selfconsistently and their difference defines the excitation energy in the polarizable environment. This procedure assumes that the states of interest and, in particular, their localization characteristics on the QM cluster are easily identifiable. Typically, this can be expected to be the case for the lowest optically active excitations in the PPEs studied here.
A two-layer scheme is employed. Within a cutoff R 1 around the QM part, atomic partial and polarizable interactions are taken into account. In the more extended buffer region with R 2 ≥ R 1 , only electrostatics are active. An example for such a partitioning is depicted in Figure 4. For a snapshot taken from the MD simulated morphology of 10-PPE with nonyl side chains in toluene, this approach is adopted using cutoffs of R 1 = 2.5 nm and R 2 = 4.0 nm. The 10-PPE backbone is treated quantum-mechanically while the side chains and solvent molecules belong to the MM region. To split the functionalized oligomer into backbone and side chains, a linkatom approach [70] with hydrogen saturation of the backbone-side chain bridge are employed. Partial charges for the solvent molecules and side chain fragments are determined from CHELPG [71] fits to the electrostatic potentials, while the atomic Thole polarizabilities are parametrized to match the molecular polarizability tensors obtained from DFT. Figure 5 shows a typical evolution of total energies of the coupled QM/MM system E QMMM during the self-consistency procedure [72]. The zero of the energy scale is defined to correspond to the total energy of the neutral 10-PPE in vacuum (iteration m = 0). It is apparent that the most significant change to the total energy in both neutral and excited states occurs during the very first step of the calculation. This is further corroborated by considering the change of total energy at iteration m compared to the previous iteration as shown in the inset of Figure 5.  iterations the respective changes are of the order of 0.01 eV and, more importantly, no significant differences are observed for the two states. Overall, the effect of polarization is small for the solvated 10-PPE and consequently, the excitation energy is nearly unaffected by the environment. More precisely, the excitation energy is 3.47 eV from the polarized QM/MM system, while a calculation using pure static MM yields 3.44 eV. Omitting the environment altogether, i.e., performing a GW -BSE calculation on the isolated oligomer conformation, yields 3.45 eV. The fact that environment effects only have a negligible impact on the calculated excitation energies can be attributed to a combination of the diluteness of the solution and the associated randomness of local electric fields and the small change of dipole moment between neutral and excited states. Similar observations have been made, e.g., for the optically excited states of push-pull oligomers embedded in a polarizable lattice [26]. Based on the above results, it is justified to limit the QM/MM setup to only electrostatic interactions in the following. Having realized that the direct electronic effects of solvent molecules on the excitations in 10-PPE are small, the focus is now on indirect effects that originate from the influence on the backbone conformations. To this end, 10-PPE in both toluene and water is considered and sample the conformations at different time intervals (∆t = 10 fs, 100 fs, 1 ps, and 10 ps), all starting from an identical starting point t 0 = 7.7 ns of our MD simulations are taken. For each of these snapshots the absorption spectrum is calculated in a static MM environment defined by R 2 = 4 nm (implies R 1 = 0 nm). The obtained discrete spectra of excitation energies and associated oscillator strengths are broadened by Gaussian functions with a FWHM (full width at half maximum) of 0.3 eV. It is found that the absorption properties are insensitive to structural dynamics of the backbone at time scales of 100 fs. Only for times exceeding about 500 fs fluctuations in the peak positions and heights of the spectra can be observed both in toluene and water. Figure 6 shows the evolution of the absorption spectrum for the time steps ∆t = 1 ps, as well as the average over the eleven respective snapshots. While the dynamics of the backbone is comparatively slow since it is to a significant extent constrained by the nonyl side chain dynamics in both poor solvents, one can observe stronger fluctuations of the absorption spectra in water as compared to toluene.
To understand the origin of these fluctuations with respect to backbone conformations in more detail, the electron-hole wave functions of the excitations is analyzed at times t = 4 ps and t = 5 ps, c.f. Figure 6(a). In the top row of Figure 7 isosurfaces for the hole (red) and electron (blue) density distributions are shown. The overall conformation of the 10-PPE exhibits a characteristic bend as a result of the stress caused by side chain interactions. At both times, the excitation appears to be localized at the apex of the bend, more pronounced for the structure at 4 ps which is lower in energy by 0.13 eV. The different characteristics can be attributed to a slightly stronger out-of-plane bent angle between the phenylene and etynlene. Over all, co-planarity of the phenyl rings along the backbone (or the lack thereof) does not appear to affect the excitations significantly.
Analysis of the composition of the electron-hole wavefunction reveals striking differences between the two snapshots. The excitation shown in Figure 7(a) is formed to 60% by a transition between the two frontier orbitals. The isosurfaces of these orbitals show that both HOMO and LUMO are practically extended along the full length of the backbone. Slight intensity variations can be noted with the HOMO being more attracted to the apex while the LUMO is thinning out at the same spot. These variations give rise to the coupled excitation being localized. At t = 5 ps, in contrast, there is not a single dominant contribution to the electron-hole wavefunction. Rather, a superposition of several transitions is found, with HOMO-1 to LUMO and HOMO to LUMO+1 transitions being most significant. As can be seen in Figure 7(b) conformational changes result in a different localization characteristics of the underlying single-particle orbitals not at the apex but left and right from it, respectively. As a pure transition between two localized states, such as the one from HOMO-1 to LUMO, is energetically penalized by stronger exchange interactions. By mixing in transitions between lower lying occupied and higher unoccupied levels, an effectively more delocalized excitation is formed. An analogous analysis of the respective excitations of 10-PPE in water, i.e, at t = 9 ps and t = 10 ps, reveals qualitatively similar behaviour.

Summary
Electronic excitations of PPE were computed using a QM/MM approach combining many-body Green's functions theory within GW approximation and the Bethe-Salpeter equation. Conformations of solvated PPE as obtained from atomistic MD simulations were used in the mixed QM/MM setup in order to determine optical excitations of solvated PPE. The reliability of optical excitations based on MM/MD conformations were investigated by comparing optical excitations of n-PPE (n = 1, 2, . . . , 10) using both optimized DFT geometries and MD geometries in vacuum. The results show that the excitation energies Ω calculated based on MM/MD conformations are larger than the ones calculated based on QM optimized geometries. For large n, the deviation amounts to 0.25 eV and is the cumulative result of geometric differences between MM/MD geometries and QM geometries. Overall agreement between the excitation energies based on MM/MD conformations and QM geometries is good enough to conclude that the use of MM/MD conformations for 10-PPE captures the relevant physico-chemical properties.
Conformations of 2,5-dinonyl-10-PPE with nonyl side chains were studied in toluene and water. The side chains were found to be dispersed from each other and from the backbone in toluene. In water, the side chains tend to aggregate. Optical excitations were calculated for 10-PPE in the QM/MM setup. The results show that the electronic environment contributions are negligible compared to the conformation dynamics of the conjugated PPE. From the analysis of the electron-hole wave function, sensitivity of energy and localization characteristics of the excited states to bends in global conformation of PPE polymer was observed.