A Case Study of Flugi-2 - American Chemical Society

The classical force field parameters were obtained and validated specifically for the purpose of the present work. The calculated gas phase. 1. Page 1...
1 downloads 0 Views 3MB Size
Subscriber access provided by UNIV AUTONOMA DE COAHUILA UADEC

A: Spectroscopy, Molecular Structure, and Quantum Chemistry

The Origin of the Solvatochromism in Organic Fluorophores With Flexible Side Chains: A Case Study of Flugi-2 Franziska Elisabeth Wolff, Sebastian Hoefener, Marcus Elstner, and Tomasz A. Wesolowski J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.9b02474 • Publication Date (Web): 02 May 2019 Downloaded from http://pubs.acs.org on May 3, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The Origin of the Solvatochromism in Organic Fluorophores with Flexible Side Chains: A Case Study of Flugi-2 Franziska E. Wolff,† Sebastian Höfener,‡ Marcus Elstner,¶ and Tomasz A. Wesołowski∗,§ †Department of Theoretical Chemical Biology, Institute of Physical Chemistry, Karlsruhe Institute of Technology, Kaiserstraße 12, 76131 Karlsruhe, Germany ‡Theoretical Chemistry Group, Institute of Physical Chemistry, Karlsruhe Institute of Technology, Kaiserstraße 12, 76131 Karlsruhe, Germany ¶Department of Theoretical Chemical Biology, Institute of Physical Chemistry and Institute of Biological Interfaces (IBG-2), Karlsruhe Institute of Technology, Kaiserstraße 12, 76131 Karlsruhe, Germany §Université de Genève, Département de Chimie Physique, 30 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland E-mail: [email protected]

Abstract The emission band for Flugi-2 solvated in DMSO is obtained from the combined quantum-classical simulations in which the QM/MM excitation energies are evaluated at the equilibrated segment of the classical molecular dynamics trajectory on the lowest excites state potential energy surface. The classical force field parameters were obtained and validated specifically for the purpose of the present work. The calculated gas phase

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to DMSO solvatochromic shifts amounts -0.21 eV which is in line with experimentally available difference between the maxima of the emission bands for Flugi-2 in decane and in DMSO (-0.26 eV). The used model describes rather well the effect on DMSO on broadening of the emission band. The solvatochromic shift in DMSO shift originates from two competing effects. The structure deformation of Flugi-2 due to the interaction with DMSO, which results in a positive contribution and the negative contribution of larger magnitude due to favorable specific interactions with the solvent. The later is dominated by a single hydrogen bond between in which the proton of N3 acts is the donor.

2

ACS Paragon Plus Environment

Page 2 of 23

Page 3 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction Fluorescence imaging is an important tool in clinical diagnostics and in biomedical research. For this method efficient fluorophores are required, which are sometimes called fluorescent pharmacophores. Burchak et al. discovered a class of efficient fluorophores using combinatorial multi-component reactions 1 . Among the discovered chromophores, the Flugi-2 molecule (IUPAC name: Methyl 3-(cyclohexylamino)-2-(4-methoxyphenyl)imidazo[1,2-a]pyridine-7carboxylate) features the largest solvent effect on emission spectra. Depending on the flexible side-chain attached to the rigid conjugated core, excitation spectra of such fluorophores in polar solvents can be strongly affected 1 . The combinatorial approach is a strategy which often results in molecules with optimized properties, but it can be cumbersome to synthesize a large number of trail compounds. It is, therefore, desirable to use numerical simulations to rationalize and ultimately predict molecules with tailored spectroscopic properties. Understanding the origin of the solvatochromic shifts is indispensable to establish an adequate computational protocol. QM/MM methods are commonly used to study solvent effect on properties of embedded molecules. Whether explicit quantum mechanical treatment should be limited to the chromophore only or include also selected solvent molecules depends strongly on the nature of the solute-solvent interactions. Even if the quantum treatment is limited to the chromophore, setting up an adequate QM/MM simulation involves choices concerning: (a) the quantum mechanical method for describing the excited state of the chromophore, (b) the embedding potential, (c) the fluctuations of the structure of the solvent structure and the chromophore. The latter one is particularly important in certain cases such as solvated Flugi molecules due to the presence of a flexible side chain(s). In particular, specific solvent-chromophore interactions and the dielectric response of the solvent must be properly accounted for in QM/MM modeling. Accurate simulation techniques have been developed to quantify the importance of such effects. For electronic excited states, specific solute-solvent interactions are accounted for explicitly in QM/MM simulations either through including the nearest solvent molecules 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

into the QM part or through the embedding potential if only the chromophore is described at the QM level. In the literature, several strategies to construct the embedding potential for modelling UV/Vis absorption spectra were proposed (see Refs. 2–8 for instance). In this work, we compute fluorescence spectra by explicitly sampling the configurational space of the system, i.e. of chromphore and solvent, on the potential energy surface of the lowest excited state. Upon excitation, the dipole moment of the chromophore changes inducing solvent reorganization. To capture this, sampling has to be extended into the nanosecondregime, which is not possible using any first-principle based quantum method like small basis set Hartree-Fock Configuration Interaction Singles (HF-CIS) or time-dependent Density Functional Theory (TD-DFT). Force field approaches, on the other hand, are usually well parametrized for ground state properties, but not able to describe excited energies and forces. We therefore derive force field parameters for the excited states surface, which is a challenging task since there are no standard procedures for the determination of atomic charges and force constants 9,10 . The excited states trajectories are then post-processed computing (de-)excitation energies. The large numbers of configurations to be computed require an efficient excited states method to be applied. To this end, the semi-empirical method - Orthogonalization Model 2 with Multireference Configurational Interaction (OM2/MRCI) 11,12 , was used. OM2/MRC method was carefully tested in previous work on biological chromomophores 11,13–18 .

Methods The used simulation protocol for solvated (and isolated) Flugi-2 consists of performing classical force-field MD simulation on the lowest excited potential energy surface followed by quantum-mechanical calculations of the vertical excitation energies along the equilibrated segment of the trajectory. Similarly as it is made commonly in simulation of the UV/Vis absorption bands, the band shape is obtained from a histogram of vertical (de-)excitation en-

4

ACS Paragon Plus Environment

Page 4 of 23

Page 5 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ergies. Radiation-less transitions are thus not accounted for because the dynamics takes place on one potential energy surface. Turning back to the used force-field parameters describing the internal degrees of freedom of the chromophore, they are derived based on reference data from CC2 calculations. The intramolecular degrees of freedom on the other hand are taken care by a non-polarizable force field. The details of the parameterization and validation of the excited-state force-field are given below.

Parameterization and validation of the force-field The reference structural data for the isolated chromophore in its lowest excited state (S1 ) were obtained from CC2 calculations 19–21 . The program package Turbomole Version 7.0 22 and the def2-TZVP 23 basis set was used for all CC2 calculations. The results obtained applying also other methods for comparison purposes, further details are given in Supporting Information. For the CC2 structure as a reference target, AmberTools15 24 was used to fit the force-field parameters describing bonds. The standard General Amber force field (GAFF) 25 parameters were refined in the iterative procedure leading to a gradual convergence of the force-field minimized structure to the reference CC2 data. The parameters are provided in the Supporting Information. The atomic charges for the first excited state of Flugi-2 at CC2 optimized geometry were obtained from the procedure based on the fit to the electrostatic potential 26 . The fit was made on 14 points per atom at the distance from the atom equal to two van der Waals radii. The potential was evaluated using the charge density obtained from LR-TDDFT calculation for the first excited state using the CAM-B3LYP 27 functional and the KOALA program package 28,29 . The obtained charges are provided in the Supporting Information. The parameters for the dihedral angle around the C8-C10 atoms are not available in General Amber Force Field. The coefficients for the fitting function F = a0 + a1 cos(2θ) + a2 cos(θ), were obtained from fitting to the reference CC2 energies obtained from a scan with 10 degree step. All other degrees of freedom were optimized for each torsional angle. In the 5

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

reference CC2 data, the double counting of the contribution to the energy due to the C8-C10 bond turning was eliminated by means of subtracting the corresponding energy contribution evaluated with the GROMACS program package 30 . The potential energy profile of the CC2 calculations and the potential profile of the parameterized dihedral are shown in Fig. 1 in the Supporting Information. All parameters of the Flugi-2 molecule can be found in the Supporting Information. The final validation of the parameters used for intramolecular degrees of freedom for the chromophore was made by comparing the OM2/MRCI 11,12 excitation energies at force-fields optimized and at the reference (CC2) geometries (see the Results section). DMSO has been parameterized by the standard procedure provided in the AmberTools15 24 and the General Amber force field (GAFF) 25 on a structure optimized with the Hartree-Fock method by using a 6-31G* basis set in Gaussian09 program package 31 . The charges were determined by the standard approach with a Hartree-Fock calculation and a Merz-Singh-Kollman-population analysis and the RESP-Fit 26 .

Classical molecular dynamics simulations For both cases - Flugi-2 in gas phase and in DMSO - force-field based molecular dynamics trajectories with the time step of 2 fs were obtained using the GROMACS package version 5.0.4 30 using the excited-state parameterized force field described in the previous sections. For solvated Flugi-2, cubic simulation box (5x5x5 nm3 ) was used. The snapshots used in statistical analysis were taken from 20 ns long equilibrated MD trajectory. The trajectory was generated according to a standard protocol consisting of the following steps. The whole box was subjected to standard MM equilibration. The first step was a minimization of the potential energy of the system with the steepest descent algorithm. 1000 steps have been performed and the maximum force tolerance was set to 100 kcal/mol. It was followed by a subsequent equilibration of the velocities of the atoms in the simulation box to a temperature of 300 K in a NVT ensemble over a simulation time of 6

ACS Paragon Plus Environment

Page 6 of 23

Page 7 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1 ns. Initial velocities have been generated by the Maxwell-Boltzmann distribution at 300K. The positions of the atoms of the chromophore have been restrained (force constant 1000 kJ/mol/nm) over the whole equilibration time. In all equilibration steps, the group cut-off scheme has been chosen for the non-bonded interaction. In the next equilibration step, the pressure of the solvent box has been adjusted to 1 bar with a Parinello-Rahman barostat in a NPT ensemble over 10 ns (τp =2.0ps, β=4.5e−5 ). The temperature was set to 300K and was scaled by the Nosé-Hoover thermostat.

Calculation of emission spectra To gain sufficient statistics, the emission spectra was computed from the equilibrated trajectories of the length of 20 ns corresponding to the NPT ensemble for 300 K, where snapshots were selected every 10 ps from the trajectory. Since this amounts to a significant number of excited states calculations, we used the computationally efficient semi-empirical OM2 Hamiltonian, as implemented within the MRCI framework in the MNDO program package 32 . This method has been extensively tested for excitation energies of small molecules 12,17 and biological chromophores 15,18 . The point of particular importance is, that OM2/MRCI is able to accurately reproduce the environment effect on excitation energies 13 . OM2/MRCI calculations were performed for 20 electrons in 20 orbitals, including single and double excitations. The obtained energies were used to generate the spectra in form of histograms. For each histogram, the is shape was fitted with a Gaussian function g(x) = a0 exp(a2 (a1 − x)2 ). In addition to the two histograms one for the gas phase trajectory and one for the condensed phase, an additional histogram was also constructed. This Gaussian fit was obtained from OM2/MRCI excitation energies which were computed for the snapshots of the solvated trajectory from which all solvent molecules have been removed. The isotropic long-range electrostatic contribution to the solvatochromism of Flugi-2 in excited state was estimated using single-point DFT/ωB97X calculations with solvent 7

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Structure of the Flugi-2 molecule. The dihedral angles used throughout the present work are defined as follows: α = ^(N2 − C8 − C10 − C14); β = ^(C9 − N3 − C16 − C17); γ = ^(C8 − C9 − N3 − H11. described by means of the conductor-like screening model (COSMO) 33 model implemented in the Gaussian09 program package 31 . The 6-31G* basis set has been used. OM2/MRCI calculations indicate that the electron density distribution changes significantly upon electronic excitation. The magnitude of the dipole moment increases from 2.1 D in the ground state to 7.7 D in the excited state. Also the DFT-derived charges collected in the Supporting Information are in line with this result. The largest change between the ground state and excited state occurs at the cyclohexane moiety. While this part is almost uncharged in the electronic ground state resulting in a weak interaction with polar solvents, it becomes slightly positively charged in the excited state resulting in an increased interaction with polar solvents. For the sake of verification, the Spectroscopic Oriented Configuration Interaction (SORCI) method 34 implemented in the ORCA program package, was also used. For the SORCI single point calculations, the split-valence basis set def2-SV(P) was used whereas the Complete Active Space included 12 electrons in 12 orbitals.

8

ACS Paragon Plus Environment

Page 8 of 23

Page 9 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1: The key structural parameters of Flugi-2 in excited state obtained form geometry optimization (either CC2 or force-field) and the average structure along the classical force filed trajectory. Vertical OM2/MRCI excitation energies obtained for these structures are given in the bottom line. For the definition of the bond length alternation parameter dBLA , see Eq. (1). hdDB i and hdSB i denote the average length of double and single bonds, respectively. Geometry optimization method Force field DFT ωB97X hdDB i [Å] 1.428 1.406 hdSB i [Å] 1.361 1.367 dBLA 0.067 0.039 α [degree] -30 -27 β [degree] -90 -87 VEEa [eV] 3.01 3.01 a Vertical fluorescence energy.

CC2 1.423 1.364 0.060 -26 -90 2.88

Results and Discussion Properties of the isolated Flugi-2 molecule and validation of the force-field parameterization We start with the analysis of the structural changes in the isolated Flugi-2 molecule upon excitation. The complete list of structural parameters of Flugi-2 in ground- and excited states is given in the Supporting Information. Upon excitation no significant changes in the bond lengths in the phenyl ring of the Flugi-2 molecule are observed, whereas the bond lengths in the imidazopyridine rings change sizeably. For example, the N1-C3 and C9-N3 single bonds decrease by about 0.04 Å, while the C3-C5, C5-C4 and C9-N1 bonds increase by 0.03 to 0.05 Å. Such a behaviour can be expected as polyenes with conjugated π-systems, e.g. retinal, invert the bond lengths after excitation, that is the double bonds are elongated and increase the single-bond character, while the single bonds are shortened and increase the double-bond character. The bond length alternation has been found to be directly correlated with the excitation energy of the chromophore 35,36 . The higher is the bond length alternation the further blue shifted is the excitation energy. For linear molecules, a good measure of bond length

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 23

alternation could be just the average of the difference of neighboring single and double bond lengths. In a more general case, like the conjugated π system in Flugi-2, the definition of a good measure involves some arbitrary choices. In the case of Flugi-2, not all bonds in the conjugate ring systems change sizeably. Only the bonds, for which the excitation affects their lengths by at least 0.02 Å were considered for quantification of the bond length alternation (see Table 3 in the Supporting Information). The bonds, for which their bond length inverts in the excited state, were used to define the bond lengths alternation parameter dBLA defined as:  dBLA =

L(N 1 − C3) + L(C9 − N 3) L(C3 − C5) + L(C5 − 4) + L(C9 − N 1) − 2 3

 . (1)

where, L(N 1 − C3) is the length of the N1-C3 bond, etc., cf. Figure 1. The geometrical parameters collected in Table 1 reveal that the dBLA decreases significantly in the excited state, in line with trends found for other chromophores 35,36 . The chosen measure of bond length alternation dBLA , being obviously arbitrary, provides a meaningful quantity to characterise structural changes in the geometry of the chromophore. Our parameterization of the force-field yields geometries which reproduce the structural parameters of the reference CC2 structure with sufficient accuracy. The force-field optimized structure yields these structural parameters which are closer or even superior to the ones obtained from alternative quantum mechanical methods based on DFT. However, the excitation energy at the force-field optimized structure is the same as the one at the DFT-ωB97X optimized geometry. They both differ from the reference value by about 0.13 eV corresponding to the CC2 optimized geometry. The difference in the excitation energy at force-fieldand CC2 optimized geometries reflect the fact that the simplicity of the force-field precludes the accurate description of small changes in the internal degrees of freedom which accompany the electronic excitation. Note that DFT-ωB97X optimized geometry is not better in this respect. We can expect, however, that this systematic geometrical effect (and conse-

10

ACS Paragon Plus Environment

Page 11 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 2: The histogram of OM2/MRCI fluorescence energies calculated for MM-MDgeometries in excited state of the gas phase MM-MD (blue) and of the MM-MD in DMSO with the electrostatic interaction of the DMSO by point charges (orange) and without (green). quently the effect on the excitation energy) remains constant and negligibly influences the solvatochromic shifts which is the primary target of the present simulations.

Classical molecular dynamics simulation of Flugi-2 in the excited state Molecular dynamics simulations were used to generate two ensembles of geometries, namely one for isolated Flugi-2 in vacuum and one for the solvated Flugi-2 molecule in DMSO. They were subsequently used to obtain the histograms shown in Figure 2. The figure shows clearly that while the maximum fluorescence emission occurs at 2.96 eV in vacuum, the band maximum is red shifted by –0.21 eV to 2.75 eV in DMSO. Although the emission spectrum of DMSO solvated Flugi-2 is available - the maximum at

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3: OM2/MRCI vertical excitation energies (S1 -S0 ) along the MM-MD trajectory on the S1 surface in the gas phase (blue) and in DMSO represented as dielectric continuum (pink), point charges (orange), or with the embedding potential switched off (green). For the sake of clarity, the graph has been smoothed along the simulation time, the running average of 100 data points is plotted here. 528 nm (2.35 eV) - a direct comparison of the solvatochromic shift obtained in our simulations with experimental data is not possible because of missing data for the gas phase. In decane, the emission band maximum measured is situated at 275 nm (2.61 eV) 1 . Taking this value as the experimental counterpart of the simulated vacuum spectrum and experimental data for DMSO solvated Flugi-2 yields the experimental solvatochromic shift equal to -0.26 eV. Taking into account a number approximations made in the applied simulation a better agreement between experimental and simulated shifts (-0.21 eV vs -0.26 eV) could hardly be expected. Besides the relatively good agreement between measured and simulated shifts, the histograms shown in Figure 2 reproduce quite fairly the measured band widths. The half-width of the emission band for DMSO solvated Flugi-2 estimated from the figures given in Ref. 1 equals 0.31 eV whereas that from the simulated spectrum is 0.42 eV. The analysis of the 12

ACS Paragon Plus Environment

Page 12 of 23

Page 13 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

structural and electrostatic contributions to the solvatochromic shifts indicates that the electrostatic interactions are responsible for the solvatochromic shifts whereas the structural effects lead to broadening of the band in DMSO. It is also worthwhile to notice the difference between the excitation energy evaluated at the force field minimum equals to 3.01 eV (see Table 1) and the center of the emission band in gas-phase simulations (2.96 eV) (see Figure 2). This small difference might can be most likely attributed to the asymmetry of the potential energy surface revealed in the finite-temperature simulations (see also Fig. 3 in the Supporting Information). The relatively good agreement between the experimentally and theoretically obtained shifts (-0.26 vs -0.21 eV) justify a more detailed analysis of the simulation in order to identify the factors determining the overall solvatochromic shift.

Geometric and electronic contributions In order to decompose the total solvatochromic shift into its electronic and geometric components, the vertical (de-)excitation energies were calculated along the previously analysed trajectory made for solvated Flugi-2 but with the embedding potential switched off. Obviously such computational experiment does not correspond to any experimental situation. Forcing the isolated Flugi-2 molecule to sample the conformational space in the same way as in the solvent, was made to determine the relative importance of the effect of the solvent on the dynamics of the chromophore. Figure 2 shows that the effect on changing the dynamics of isolated chromophore on the position of the maximum of the emission band is small but not negligible. The maximum shifts from 2.96 eV (dynamics as in the gas phase) to 3.05 eV (dynamics as in the solvent), i.e., in the opposite direction then the whole solvatochromic shift. The contribution to the shift due to the electric field generated by the solvent is thus quite significant. It results in red-shift of the histogram by about 0.30 eV. Not only it compensates the geometric effect, but its magnitude is about three times larger. In the following part, we consider simpler treatments of the environment to verify if

13

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4: Iso-surface of the transition density upon the lowest singlet excitation computed at the excited-state geometry using RICC2. Red indicates the density increase. The lowest singlet excitation has 86% HOMO-LUMO character. For the picture of these orbitals, see Supporting Information. they could effectively describe the electrostatic effects. We start with noticing that most significant changes in the electron density distribution due to (de-)excitation are localized on the cyclohexane moiety (decrease) and on the N3 atom (increase) (see Figure 4). This observation, together with the significant contribution of the electrostatic field, suggests that even more approximate models such COSMO could be applied in this case. Such continuum models are not able to describe anisotropic interactions but they take into account the solvent polarization by means of the changes in the effective surface charges on the dielectric boundary. Figure 3 shows how the implicit solvation model COSMO performs in this case. The excitation energy for isolated Flugi-2 at conformations taken the solvatedphase trajectory evaluated with (orange) and without (green) COSMO treatment of the solvent are almost identical. This indicates that explicit treatment of specific solute-solvent interactions is indispensable in explaining large solvatochromism of emission spectra of Flugi2. The results obtained using COSMO treatment of the solvent and excitation energies from SORCI (see Supporting Information) confirm this conclusion.

14

ACS Paragon Plus Environment

Page 14 of 23

Page 15 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) Dihedral angles: α = −29 degree, β = −23 degree, γ = 149 degree

(b) Dihedral angles: α = 19, β = −131 degree, γ = 31 degree

Figure 5: Two representative snapshots taken from the 20 ns simulation featuring the key hydrogen bond between Flugi-2 and a DMSO molecule. Side chains and dihedral angles The final part of the analysis concerns possible relation between the structure of the flexible sidechains and the excitation energy. A detailed analysis of the trajectory reveals a hydrogen bond between the oxygen atom of a DMSO molecule and the N3 hydrogen atom of the Flugi2 molecule. This hydrogen bond is stable over the 20 ns simulation time and influences the C8-C10 and N3-C16 dihedral angles quite significantly, see Figure 5 and Figure 6. This bond, in turn, affects the C8-C10 and N3-C16 dihedral angles. A detailed cluster analysis shows the correlation between the dihedral angles and the excitation energy in DMSO, see Figures 7c and 7d. The centroid points for each cluster - black points in Figure 7 - mark the mean value of the excitation energy for each cluster. The dihedral angle distribution in the trajectory is significantly affected by the interactions with DMSO. The sampled part of the conformational space corresponding to variation of this dihedral angle is different in isolated and in solvated Flugi-2. The distribution of the dihedral angles in gas phase trajectory is less spread out. The dihedral angles in gas phase

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6: α and γ dihedral angles over 20 ns simulation time. do not fluctuate as strongly as in the DMSO trajectory, because of the lacking interaction with DMSO molecules. Furthermore, the hydrogen bond to the DMSO enables the Flugi-2 to sample more dihedral angles than in gas phase. Comparing the gas phase and DMSO trajectory, the distribution of the dihedral angle β is shifted to lower dihedral angles in the presence of DMSO, while the region between 0 to 135 degree is not sampled at all, see Figure 7d. For the α dihedral angle, different areas of the phase-space can be sampled by adding DMSO as a solvent, see Figure 7a and Figure 7c. The different distribution of the dihedral angles results in a more distributed excitation energy in the individual clusters of the α and β dihedral angles in the DMSO trajectory, see Figure 7. While in gas phase trajectory, the cluster center vary for the dihedral angles

16

ACS Paragon Plus Environment

Page 16 of 23

Page 17 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

only up to 0.1 eV (for the β Figure 7b). In DMSO trajectory the shift between the cluster centroids is up to 0.2 eV (for the β Figure 7d). Comparing the centroid points of the clusters with the same dihedral angle between the gas phase Figure 7b and the DMSO trajectory, Figure 7d, shows that the red shift can not be exclusively attributed to the dihedral angles. The centroid points of the same clusters in gas phase and DMSO, for example the cluster around –45 degree in Figure 7b and Figure 7d, do not correspond to the same excitation energy.

Conclusions In the present work, the solvatochromic shift upon fluorescence emission for the Flugi-2 compound solvated in DMSO is investigated using a force field that has been parameterized exclusively for the purpose to model Flugi-2 in the excited state. The parametrization of the force-field is a key element of the simulations reported in this work. The CC2 calculations revealed the inversion of the bond length alternation in the excited state of Flugi-2, which must be taken into account in the description of the internal degrees of freedom of the chromophore. The solvent-chromophore interactions result in a red-shift of the magnitude of 0.21 eV which agrees quite well with its experimental counterpart - the red-shift of the maximum of the emission band between decane and DMSO amounting to 0.26 eV. The red-shift results mainly from direct electronic interactions which are compensated only partially by the effect of the solvent on the dynamics of the chromophore. We demonstrated that quasi-equilibrium simulations in which the chromophore is forced to remain in its lowest excited state in thermal equilibrium with the environment can describe quite accurately not only the solvatochromic shifts but also the widths of the bands. The dynamics of the relaxation processes can obviously not be described using the proposed protocol. The excitation involves a large reorganization of electron density. The cyclohexane moi-

17

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ety becomes slightly positively charged in the excited state. This suggests that continuum dielectric models could be adequate. Our comparative simulations using explicit and implicit models of the solvent indicate, however, that the solvatochromism of the emission spectra cannot be properly accounted with the latter. Taking into account specific solventchromophore interactions is indispensable. Among these interactions, the formation of the N3-H–O hydrogen bond in the excited state was identified as the one responsible for the observed red-shift. Further improvements, leading agreement in absolute positions as well as the shapes of the emission bands would require significant refinement of the applied combined classicalquantum simulation protocol. For instance, better description of the potential energy surface of the isolated chromophore at the excited state would require significantly more expensive quantum mechanical methods then the ones used in the present work, better description of the embedding potential would require taking into account mutual polarisation of the chromophore and the solvent and the Pauli repulsion. Such refinements although possible lie outside of the scope of the present work which aimed at the determination of the most likely molecular mechanism responsible for large solvatochromic shift in emission of Flugi-2.

Acknowledgement The authors thank the Deutsche Forschungsgemeinschaft for funding via GRK 2039 and acknowledge support by the state of Baden-Württemberg through bwHPC and the bwHPC project provided through associated compute services of the JUSTUS HPC facility at the University of Ulm and the German Research Foundation (DFG) through grant no INST 40/467-1 FUGG. This research was supported by grant from Swiss National Science Foundation (Grant No. 200020-172532).

18

ACS Paragon Plus Environment

Page 18 of 23

Page 19 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Supporting Information Available Optimized geometries of the chromophore using various quantum mechanical methods together with the corresponding vertical excitation energies, atomic charges, parameters of the force field, orbitals involved in electronic transitions, histograms of the vertical excitation energies along the molecular dynamics trajectory.

References (1) Burchak, O. N.; Mugherli, L.; Ostuni, M.; Lacapere, J. J.; Balakirev, M. Y. Journal of the American Chemical Society 2011, 133, 10058–10061. (2) Wesolowski, T. A. J. Am. Chem. Soc. 2004, 126, 11444–11445. (3) Nielsen, C. B.; Christiansen, O.; Mikkelsen, K. V.; Kongsted, J. Journal of Chemical Physics 2007, 126, 154112. (4) Defusco, A.; Minezawa, N.; Slipchenko, L. V.; Zahariev, F.; Gordon, M. S. Journal of Physical Chemistry Letters 2011, 2, 2184–2192. (5) Daday, C.; Koenig, C.; Valsson, O.; Neugebauer, J.; Filippi, C. J. Chem. Theory Comput. 2013, 9, 2355–2367. (6) Wesolowski, T. A. J. Chem. Phys. 2014, 140, 18A530. (7) Mennucci, B. International Journal of Quantum Chemistry 2015, 115, 1202–1208. (8) Gurunathan, P. K.; Acharya, A.; Ghosh, D.; Kosenkov, D.; Kaliman, I.; Shao, Y.; Krylov, A. I.; Slipchenko, L. V. Journal of Physical Chemistry B 2016, 120, 6562– 6574. (9) Graen, T.; Hoefling, M.; Grubmüller, H. Journal of Chemical Theory and Computation 2014, 10, 5505–5512. 19

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(10) Vaiana, A. C.; Schulz, A.; Wolfrum, J.; Sauer, M.; Smith, J. C. Journal of Computational Chemistry 2003, 24, 632–639. (11) Dral, P. O.; Wu, X.; Spörkel, L.; Koslowski, A.; Weber, W.; Steiger, R.; Scholten, M.; Thiel, W. Journal of Chemical Theory and Computation 2016, 12, 1082–1096. (12) Weber, W.; Thiel, W. Theoretical Chemistry Accounts 2000, 103, 495–506. (13) Guo, Y.; Beyle, F. E.; Bold, B. M.; Watanabe, H. C.; Koslowski, A.; Thiel, W.; Hegemann, P.; Marazzi, M.; Elstner, M. Chem. Sci. 2016, 7, 3879–3891. (14) Wanko, M.; Hoffmann, M.; Frauenheim, T.; Elstner, M. Journal of Computer-Aided Molecular Design 2006, 20, 511–518. (15) Wanko, M.; Hoffmann, M.; Strodel, P.; Koslowski, A.; Thiel, W.; Neese, F.; Frauenheim, T.; Elstner, M. Journal of Physical Chemistry B 2005, 109, 3606–3615. (16) Hoffmann, M.; Wanko, M.; Strodel, P.; König, P. H.; Frauenheim, T.; Schulten, K.; Thiel, W.; Tajkhorshid, E.; Elstner, M. Journal of the American Chemical Society 2006, 128, 10808–10818. (17) Silva-Junior, M. R.; Thiel, W. Journal of Chemical Theory and Computation 2010, 6, 1546–64. (18) Sun, Q.; Li, Z.; Lan, Z.; Pfisterer, C.; Doerr, M.; Fischer, S.; Smith, S. C.; Thiel, W. Physical chemistry chemical physics : PCCP 2012, 14, 11413–24. (19) Köhn, A.; Hättig, C. The Journal of Chemical Physics 2003, 119, 5021–5036. (20) Hättig, C.; Köhn, A. The Journal of Chemical Physics 2002, 117, 6939–6951. (21) Hättig, C. The Journal of Chemical Physics 2003, 118, 7751–7761. (22) TURBOMOLE V6.2 2010,

a development of University of Karlsruhe and

Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; 20

ACS Paragon Plus Environment

Page 20 of 23

Page 21 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

available from http://www.turbomole.com. (23) Yousaf, K. E.; Peterson, K. A. The Journal of Chemical Physics 2008, 129, 184108. (24) Case, D. et al. AMBER 2015, University of California, San Francisco. 2015. (25) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Journal of Computational Chemistry 2004, 25, 1157–1174. (26) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. The Journal of Physical Chemistry 1993, 97, 10269–10280. (27) Yanai, T.; Tew, D. P.; Handy, N. C. Chemical Physics Letters 2004, 393, 51–57. (28) KOALA, an ab-initio electronic structure program, written by S. Höfener, with contributions from A.-S. Hehn, J. Heuser and N. Schieschke. (29) Höfener, S. J. Comput. Chem. 2014, 35, 1716–1724. (30) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435–447. (31) Gaussian 09, Revision A.02, Gaussian, Inc., Wallingford CT, 2016. (32) Thiel, W. MNDO99 Program Package, Version 6, Max–Planck–Institute für Kohlenforschung Mülheim an der Ruhr, Germany. 2006. (33) Klamt, A.; Schüürmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 799–805. (34) Neese, F. The Journal of Chemical Physics 2003, 119, 9428–9443. (35) Hoffmann, M.; Wanko, M.; Strodel, P.; König, P. H.; Frauenheim, T.; Schulten, K.; Thiel, W.; Tajkhorshid, E.; Elstner, M. Journal of the American Chemical Society 2006, 128, 10808–10818. 21

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(36) Kochendoerfer, G. G.; Wang, Z.; Oprian, D. D.; Mathies, R. A. Biochemistry 1997, 36, 6577–6587.

22

ACS Paragon Plus Environment

Page 22 of 23

Page 23 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) Cluster analysis α dihedral angle in gas phase.

(b) Cluster analysis β dihedral angle in gas phase.

(c) Cluster analysis α dihedral angle in DMSO.

(d) Cluster analysis β dihedral angle in DMSO.

Figure 7: Vertical fluorescence energies plotted against the dihedral angle of α (a) and (c) and β (b) and (d) in gas phase (top) and DMSO (bottom). The data points were classified and the centroid points were obtained using a k-means cluster algorithm. Each color represents one determined cluster and the black dots are the calculated cluster centers, which represent the average of each cluster.

23

ACS Paragon Plus Environment