Convergence of Computed Aqueous Absorption Spectra with Explicit Quantum Mechanical Solvent

Two-Dimensional Electronic Spectroscopy in the Condensed Phase Using Equivariant Transformer Accelerated Molecular Dynamics Simulations

Two-dimensional electronic spectroscopy (2DES) provides rich information about how the electronic states of molecules, proteins, and solid-state materials interact with each other and their surrounding environment. Atomistic molecular dynamics simulations offer an appealing route to uncover how nuclear motions mediate electronic energy relaxation and their manifestation in electronic spectroscopies but are computationally expensive. Here we show that by using an equivariant transformer-based machine learning architecture trained with only 2500 ground state and 100 excited state electronic structure calculations, one can construct accurate machine-learned potential energy surfaces for both the ground-state electronic surface and excited-state energy gap. We demonstrate the utility of this approach for simulating the dynamics of Nile blue in ethanol, where we experimentally validate and decompose the simulated 2DES to establish the nuclear motions of the chromophore and the solvent that couple to the excited state, connecting the spectroscopic signals to their molecular origin.

Modeling electronic absorption spectra with nuclear quantum effects in constrained nuclear-electronic orbital framework

Electronic absorption spectra serve as versatile and powerful tools in experiments. Accurate theoretical simulation of electronic absorption spectra is challenging because multiple factors such as environmental effects and nuclear quantum effects contribute to spectrum lineshapes. This work proposes a protocol to model electronic absorption spectra in the constrained nuclear-electronic orbital framework. Solvent effects, temperature effects, and particularly nuclear quantum effects can be taken into consideration in this unified framework. This protocol is applied to investigate the electronic absorption spectrum of the pyridine molecule in water. Nuclear quantum effects are found to induce a broadening and red shift of the absorption spectrum of pyridine.

How Accurate Are QM/MM Models?

Despite the success and widespread use of QM/MM methods in modeling (bio)chemically important processes, their accuracy is still not well understood. A key reason is because these methods are ultimately approximations to direct QM calculations of very large systems, which are impractical to perform in most cases. We highlight recent progress toward the development of realistic model systems where it is possible to obtain full QM reference data to directly and systematically evaluate the effectiveness of different QM/MM generation schemes. These model systems are highly flexible and can be tailored to probe the sensitivity of a QM/MM model to different reaction types and simulation parameters such as pairing of QM and MM potentials, QM region size, and composition. It is envisaged that this strategy could be used to directly validate different QM/MM generation schemes and spur the development of more robust models in the future.

Toward a numerically efficient description of bulk-solvated anionic states

We investigate the vertical electron attachment energy (VAE) of 1-methyl-4-nitroimidazole, a model radiosensitizer, employing quantum mechanics/molecular mechanics (QM/MM) and QM/polarized continuum (QM/PCM) solvation models. We considered the solvent-excluded surface (QM/PCM-SES) and Van der Waals (QM/PCM-VDW) cavities within the PCM framework, the electrostatic embedding QM/MM (EE-QM/MM) model, and the self-consistent sequential QM/MM polarizable electrostatic embedding (scPEE-S-QM/MM) model. Due to slow VAE convergence concerning the number of QM solvent molecules, full QM calculations prove inefficient. Ensemble averages in these calculations do not align with VAEs computed for the representative solute-solvent configuration. QM/MM and QM/PCM calculations show agreement with each other for sufficiently large QM regions, although the QM/PCM-VDW model exhibits artifacts linked to the cavity. QM/MM models demonstrate good agreement between ensemble averages and VAEs calculated with the representative configuration. Notably, the VAE computed with the scPEE-S-QM/MM model achieves faster convergence concerning the number of QM water molecules compared to the EE-QM/MM model, attributed to enhanced efficiency from MM charge polarization in the scPEE-S-QM/MM approach. This emphasizes the importance of QM/classical models with accurate solute-solvent and solvent-solvent mutual polarization for obtaining converged VAEs at a reasonable computational cost. The full-QM approach is very inefficient, while the microsolvation model is inaccurate. Computational savings in QM/MM models result from electrostatic embedding and the representative configuration, with the scPEE-S-QM/MM approach emerging as an efficient tool for describing bulk-solvated anions within the QM/MM framework. Its potential extends to improving transient anion state descriptions in biomolecules and radiosensitizers, especially given the frequent employment of microsolvation models.

First-Principles Modeling of the Absorption Spectrum of Crystal Violet in Solution: The Importance of Environmentally Driven Symmetry Breaking

Theoretical spectroscopy plays a crucial role in understanding the properties of the materials and molecules. One of the most promising methods for computing optical spectra of chromophores embedded in complex environments from the first principles is the cumulant approach, where both (generally anharmonic) vibrational degrees of freedom and environmental interactions are explicitly accounted for. In this work, we verify the capabilities of the cumulant approach in describing the effect of complex environmental interactions on linear absorption spectra by studying Crystal Violet (CV) in different solvents. The experimental absorption spectrum of CV strongly depends on the nature of the solvent, indicating strong coupling to the condensed-phase environment. We demonstrate that these changes in absorption line shape are driven by an increased splitting between absorption bands of two low-lying excited states that is caused by a breaking of the D3 symmetry of the molecule and that in polar solvents, this symmetry breaking is mainly driven by electrostatic interactions with the condensed-phase environment rather than distortion of the structure of the molecule, in contrast with conclusions reached in a number of previous studies. Our results reveal the importance of explicitly including a counterion in the calculations in nonpolar solvents due to electrostatic interactions between CV and the ion. In polar solvents, these interactions are strongly reduced due to solvent screening effects, thus minimizing the symmetry breaking. Computed spectra in methanol are found to be in reasonable agreement with the experiment, demonstrating the strengths of the outlined approach in modeling strong environmental interactions.

Photophysics of a nucleic acid-protein crosslinking model strongly depends on solvation dynamics: an experimental and theoretical study

We present a combined experimental and theoretical study of the photophysics of 5-benzyluracil (5BU) in methanol, which is a model system for interactions between nucleic acids and proteins. A molecular dynamics study of 5BU in solution through efficient DFT-based hybrid ab initio potentials revealed a remarkable conformational flexibility - allowing the population of two main conformers - as well as specific solute-solvent interactions, which both appear as relevant factors for the observed 5BU optical absorption properties. The simulated absorption spectrum, calculated on such an ensemble, enabled a molecular interpretation of the experimental UV-Vis lowest energy band, which is also involved in the induced photo-reactivity upon irradiation. In particular, the first two excited states (mainly involving the uracil moiety) both contribute to the 5BU lowest energy absorption. Moreover, as a key finding, the nature and brightness of such electronic transitions are strongly influenced by 5BU conformation and the microsolvation of its heteroatoms.

Beyond the Condon limit: Condensed phase optical spectra from atomistic simulations

While dark transitions made bright by molecular motions determine the optoelectronic properties of many materials, simulating such non-Condon effects in condensed phase spectroscopy remains a fundamental challenge. We derive a Gaussian theory to predict and analyze condensed phase optical spectra beyond the Condon limit. Our theory introduces novel quantities that encode how nuclear motions modulate the energy gap and transition dipole of electronic transitions in the form of spectral densities. By formulating the theory through a statistical framework of thermal averages and fluctuations, we circumvent the limitations of widely used microscopically harmonic theories, allowing us to tackle systems with generally anharmonic atomistic interactions and non-Condon fluctuations of arbitrary strength. We show how to calculate these spectral densities using first-principles simulations, capturing realistic molecular interactions and incorporating finite-temperature, disorder, and dynamical effects. Our theory accurately predicts the spectra of systems known to exhibit strong non-Condon effects (phenolate in various solvents) and reveals distinct mechanisms for electronic peak splitting: timescale separation of modes that tune non-Condon effects and spectral interference from correlated energy gap and transition dipole fluctuations. We further introduce analysis tools to identify how intramolecular vibrations, solute-solvent interactions, and environmental polarization effects impact dark transitions. Moreover, we prove an upper bound on the strength of cross correlated energy gap and transition dipole fluctuations, thereby elucidating a simple condition that a system must follow for our theory to accurately predict its spectrum.

Time-Dependent Density Functional Theory for X-ray Absorption Spectra: Comparing the Real-Time Approach to Linear Response

We simulate X-ray absorption spectra at elemental K-edges using time-dependent density functional theory (TDDFT) in both its conventional linear-response implementation and its explicitly time-dependent or "real-time" formulation. Real-time TDDFT simulations enable broadband spectra calculations without the need to invoke frozen occupied orbitals ("core/valence separation"), but we find that these spectra are often contaminated by transitions to the continuum that originate from lower-energy core and semicore orbitals. This problem becomes acute in triple-ζ basis sets, although it is sometimes sidestepped in double-ζ basis sets. Transitions to the continuum acquire surprisingly large dipole oscillator strengths, leading to spectra that are difficult to interpret. Meaningful spectra can be recovered by means of a filtering technique that decomposes the spectrum into contributions from individual occupied orbitals, and the same procedure can be used to separate L- and K-edge spectra arising from different elements within a given molecule. In contrast, conventional linear-response TDDFT requires core/valence separation but is free of these artifacts. It is also significantly more efficient than the real-time approach, even when hundreds of individual states are needed to reproduce near-edge absorption features and even when Padé approximants are used to reduce the real-time simulations to just 2-4 fs of time propagation. Despite the cost, the real-time approach may be useful to examine the validity of the core/valence separation approximation.

Effect of the QM Size, Basis Set, and Polarization on QM/MM Interaction Energy Decomposition Analysis

Herein, an Energy Decomposition Analysis (EDA) scheme extended to the framework of QM/MM calculations in the context of electrostatic embeddings (QM/MM-EDA) including atomic charges and dipoles is applied to assess the effect of the QM region size on the convergence of the different interaction energy components, namely, electrostatic, Pauli, and polarization, for cationic, anionic, and neutral systems interacting with a strong polar environment (water). Significant improvements are found when the bulk solvent environment is described by a MM potential in the EDA scheme as compared to pure QM calculations that neglect bulk solvation. The predominant electrostatic interaction requires sizable QM regions. The results reported here show that it is necessary to include a surprisingly large number of water molecules in the QM region to obtain converged values for this energy term, contrary to most cluster models often employed in the literature. Both the improvement of the QM wave function by means of a larger basis set and the introduction of polarization into the MM region through a polarizable force field do not translate to a faster convergence with the QM region size, but they lead to better results for the different interaction energy components. The results obtained in this work provide insight into the effect of each energy component on the convergence of the solute-solvent interaction energy with the QM region size. This information can be used to improve the MM FFs and embedding schemes employed in QM/MM calculations of solvated systems.

Polarizable Force Field for Acetonitrile Based on the Single-Center Multipole Expansion

A polarizable potential function describing the interaction between acetonitrile molecules is introduced. The molecules are described as rigid and linear, with three mass sites corresponding to the CH₃ group (methyl, Me), the central carbon atom (C), and the nitrogen atom (N). The electrostatic interaction is represented using a single-center multipole expansion as has been done previously for H₂O [Wikfeldt et al., Phys. Chem. Chem. Phys. 15, 16542 (2013)], by including multipole moments from dipole up to and including hexadecapole, as well as anisotropic dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole polarizability tensors. The model is free of point charges. The non-electrostatic part is described in a pair-wise fashion by a Born-Mayer repulsion and damped dispersion attraction. The potential function is parameterized to fit the interaction energy of small (CH₃CN)_n, n = 2-6, clusters calculated using the PBE0 hybrid functional with an additional atomic many-body dispersion contribution. The parameterized potential function is found to compare well with results of the electronic structure calculations of dissociation curves for different dimer orientations and cohesive properties (the equilibrium volume, cohesive energy, and the bulk modulus) of the α-phase of acetonitrile crystal. The average value of the molecular dipole moment obtained in the α-phase is 5.53 D, corresponding to ca. 40% increase as compared to the dipole moment of an isolated acetonitrile molecule, 3.92 D. The calculated densities of solid and liquid acetonitrile turn out to be 8-10% higher than experimental values. This appears to be caused by an overestimate of the atomic many-body dispersion interaction in the density functional calculations used as input in the parametrization of the potential function.

Fragment-Based Quantum Mechanical Calculation of Excited-State Properties of Fluorescent RNAs

Fluorescent RNA aptamers have been successfully applied to track and tag RNA in a biological system. However, it is still challenging to predict the excited-state properties of the RNA aptamer–fluorophore complex with the traditional electronic structure methods due to expensive computational costs. In this study, an accurate and efficient fragmentation quantum mechanical (QM) approach of the electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC) scheme was applied for calculations of excited-state properties of the RNA aptamer–fluorophore complex. In this method, the excited-state properties were first calculated with one-body fragment quantum mechanics/molecular mechanics (QM/MM) calculation (the excited-state properties of the fluorophore) and then corrected with a series of two-body fragment QM calculations for accounting for the QM effects from the RNA on the excited-state properties of the fluorophore. The performance of the EE-GMFCC on prediction of the absolute excitation energies, the corresponding transition electric dipole moment (TEDM), and atomic forces at both the TD-HF and TD-DFT levels was tested using the Mango-II RNA aptamer system as a model system. The results demonstrate that the calculated excited-state properties by EE-GMFCC are in excellent agreement with the traditional full-system time-dependent ab initio calculations. Moreover, the EE-GMFCC method is capable of providing an accurate prediction of the relative conformational excited-state energies for different configurations of the Mango-II RNA aptamer system extracted from the molecular dynamics (MD) simulations. The fragmentation method further provides a straightforward approach to decompose the excitation energy contribution per ribonucleotide around the fluorophore and then reveals the influence of the local chemical environment on the fluorophore. The applications of EE-GMFCC in calculations of excitation energies for other RNA aptamer–fluorophore complexes demonstrate that the EE-GMFCC method is a general approach for accurate and efficient calculations of excited-state properties of fluorescent RNAs.

The Influence of Electronic Polarization on Nonlinear Optical Spectroscopy

The environment surrounding a chromophore can dramatically affect the energy absorption and relaxation process, as manifested in optical spectra. Simulations of nonlinear optical spectroscopy, such as two-dimensional electronic spectroscopy (2DES) and transient absorption (TA), will be influenced by the computational model of the environment. We here compare a fixed point charge molecular mechanics model and a quantum mechanical (QM) model of the environment in computed 2DES and TA spectra of Nile red in water and the chromophore of photoactive yellow protein (PYP) in water and protein environments. In addition to simulating these nonlinear optical spectra, we directly juxtapose the computed excitation energy correlation function to the dynamic Stokes shift function often used to analyze environment dynamics. Overall, we find that for the three systems studied here the mutual electronic polarization provided by the QM environment manifests in broader 2DES signals, as well as a larger reorganization energy and a larger static Stokes shift due to stronger coupling between the chromophore and the environment.

Multiscale Modeling of Electronic Spectra Including Nuclear Quantum Effects

Theoretical prediction of electronic absorption spectra without input from experiments is no easy feat, as it requires addressing all of the factors that affect line shapes. In practice, however, the methodologies are limited to treat these ingredients only to a certain extent. Here, we present a multiscale protocol that addresses the temperature, solvent, and nuclear quantum effects as well as anharmonicity and the reconstruction of the final spectra from individual transitions. First, quantum mechanics/molecular mechanics (QM/MM) molecular dynamics is conducted to obtain trajectories of solute-solvent configurations, from which the corresponding quantum-corrected ensembles are generated through the generalized smoothed trajectory analysis (GSTA). The optical spectra of the ensembles are then produced by calculating vertical transitions using time-dependent density-functional theory (TDDFT) with implicit solvation. To obtain the final spectral shapes, the stick spectra from TDDFT are convoluted with Gaussian kernels where the half-widths are determined by a statistically motivated strategy. We have tested our method by calculating the UV-vis spectra of a recently discovered acridine photocatalyst in two redox states. Vibronic progressions and broadenings due to the finite lifetime of the excited states are not included in the methodology yet. Nuclear quantization affects the relative peak intensities and widths, which is necessary to reproduce the experimental spectrum. We have also found that using only the optimized geometry of each molecule works surprisingly well if a proper empirical broadening factor is applied. This is explained by the rigidity of the conjugated chromophore moieties of the selected molecules, which are mainly responsible for the excitations in the spectra. In contrast, we have also shown that other parts of the molecules are flexible enough to feature anharmonicities that impair the use of other techniques such as Wigner sampling.

Generalized energy-based fragmentation approach for calculations of solvation energies of large systems

A generalized energy-based fragmentation (GEBF) approach has been combined with a universal solvation model based on solute electron density (SMD) to compute the solvation energies of general large systems (such as protein molecules) in solutions. In the GEBF-SMD method, the solvation energy of a target system could be combined by the corresponding solvation energies of various subsystems, each of which is embedded in the background point charges and surface charges on the surface of solute cavity at the positions of its atoms and neighbouring atoms outside of the subsystem. Our results show that the GEBF-SMD model could reproduce the conventional SMD solvation energies quite well for various proteins in solutions, and could significantly reduce the computational costs for the SMD calculations of large proteins. In addition, the GEBF-SMD approach is almost independent of the basis sets and the types of solvents (including protic, polar, and nonpolar ones). Also, the GEBF-SMD approach could reproduce the relative energies of various conformers of large systems in solutions. Therefore, the GEBF-SMD method is expected to be applicable for computing the solvation energies of a broad range of large systems.

Benzene Radical Anion Microsolvated in Ammonia Clusters: Modeling the Transition from an Unbound Resonance to a Bound Species

The benzene radical anion, well-known in organic chemistry as the first intermediate in the Birch reduction of benzene in liquid ammonia, exhibits intriguing properties from the point of view of quantum chemistry. Notably, it has the character of a metastable shape resonance in the gas phase, while measurements in solution find it to be experimentally detectable and stable. In this light, our previous calculations performed in bulk liquid ammonia explicitly reveal that solvation leads to stabilization. Here, we focus on the transition of the benzene radical anion from an unstable gas-phase ion to a fully solvated bound species by explicit ionization calculations of the radical anion solvated in molecular clusters of increasing size. The computational cost of the largest systems is mitigated by combining density functional theory with auxiliary methods including effective fragment potentials or approximating the bulk by polarizable continuum models. Using this methodology, we obtain the cluster size dependence of the vertical binding energy of the benzene radical anion converging to the value of -2.3 eV at a modest computational cost.

Vibronic and Environmental Effects in Simulations of Optical Spectroscopy

Including both environmental and vibronic effects is important for accurate simulation of optical spectra, but combining these effects remains computationally challenging. We outline two approaches that consider both the explicit atomistic environment and the vibronic transitions. Both phenomena are responsible for spectral shapes in linear spectroscopy and the electronic evolution measured in nonlinear spectroscopy. The first approach utilizes snapshots of chromophore-environment configurations for which chromophore normal modes are determined. We outline various approximations for this static approach that assumes harmonic potentials and ignores dynamic system-environment coupling. The second approach obtains excitation energies for a series of time-correlated snapshots. This dynamic approach relies on the accurate truncation of the cumulant expansion but treats the dynamics of the chromophore and the environment on equal footing. Both approaches show significant potential for making strides toward more accurate optical spectroscopy simulations of complex condensed phase systems.

Quantum and semiclassical dynamical studies of nonadiabatic processes in solution: achievements and perspectives

We concisely review the main methodological approaches to model nonadiabatic dynamics in isotropic solutions and their applications. Three general classes of models are identified as the most used to include solvent effects in the simulations. The first model describes the solvent as a set of harmonic collective modes coupled to the solute degrees of freedom, and the second as a continuum, while the third explicitly includes solvent molecules in the calculations. The issues related to the use of these models in semiclassical and quantum dynamical simulations are discussed, as well as the main limitations and perspectives of each approach.

Thermally activated delayed fluorescence: A critical assessment of environmental effects on the singlet-triplet energy gap

The effective design of dyes optimized for thermally activated delayed fluorescence (TADF) requires the precise control of two tiny energies: the singlet-triplet gap, which has to be maintained within thermal energy, and the strength of spin-orbit coupling. A subtle interplay among low-energy excited states having dominant charge-transfer and local character then governs TADF efficiency, making models for environmental effects both crucial and challenging. The main message of this paper is a warning to the community of chemists, physicists, and material scientists working in the field: the adiabatic approximation implicitly imposed to the treatment of fast environmental degrees of freedom in quantum-classical and continuum solvation models leads to uncontrolled results. Several approximation schemes were proposed to mitigate the issue, but we underline that the adiabatic approximation to fast solvation is inadequate and cannot be improved; rather, it must be abandoned in favor of an antiadiabatic approach.

Explicit environmental and vibronic effects in simulations of linear and nonlinear optical spectroscopy

Accurately simulating the linear and nonlinear electronic spectra of condensed phase systems and accounting for all physical phenomena contributing to spectral line shapes presents a significant challenge. Vibronic transitions can be captured through a harmonic model generated from the normal modes of a chromophore, but it is challenging to also include the effects of specific chromophore-environment interactions within such a model. We work to overcome this limitation by combining approaches to account for both explicit environment interactions and vibronic couplings for simulating both linear and nonlinear optical spectra. We present and show results for three approaches of varying computational cost for combining ensemble sampling of chromophore-environment configurations with Franck-Condon line shapes for simulating linear spectra. We present two analogous approaches for nonlinear spectra. Simulated absorption spectra and two-dimensional electronic spectra (2DES) are presented for the Nile red chromophore in different solvent environments. Employing an average Franck-Condon or 2DES line shape appears to be a promising method for simulating linear and nonlinear spectroscopy for a chromophore in the condensed phase.

Quantifying Uncertainties in Solvation Procedures for Modeling Aqueous Phase Reaction Mechanisms

Computational quantum chemistry provides fundamental chemical and physical insights into solvated reaction mechanisms across many areas of chemistry, especially in homogeneous and heterogeneous renewable energy catalysis. Such reactions may depend on explicit interactions with ions and solvent molecules that are nontrivial to characterize. Rigorously modeling explicit solvent effects with molecular dynamics usually brings steep computational costs while the performance of continuum solvent models such as polarizable continuum model (PCM), charge-asymmetric nonlocally determined local-electric (CANDLE), conductor-like screening model for real solvents (COSMO-RS), and effective screening medium method with the reference interaction site model (ESM-RISM) are less well understood for reaction mechanisms. Here, we revisit a fundamental aqueous hydride transfer reaction-carbon dioxide (CO₂) reduction by sodium borohydride (NaBH₄)-as a test case to evaluate how different solvent models perform in aqueous phase charge migrations that would be relevant to renewable energy catalysis mechanisms. For this system, quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations almost exactly reproduced energy profiles from QM simulations, and the Na⁺ counterion in the QM/MM simulations plays an insignificant role over ensemble averaged trajectories that describe the reaction pathway. However, solvent models used on static calculations gave much more variability in data depending on whether the system was modeled using explicit solvent shells and/or the counterion. We pinpoint this variability due to unphysical descriptions of charge-separated states in the gas phase (i.e., self-interaction errors), and we show that using more accurate hybrid functionals and/or explicit solvent shells lessens these errors. This work closes with recommended procedures for treating solvation in future computational efforts in studying renewable energy catalysis mechanisms.

Modeling Excited-State Proton Transfer to Solvent: A Dynamics Study of a Super Photoacid with a Hybrid Implicit/Explicit Solvent Model

The rapid growth of time-resolved spectroscopies and the theoretical advances in ab initio molecular dynamics (AIMD) pave the way to look at the real-time molecular motion following the electronic excitation. Here, we exploited the capabilities of AIMD combined with a hybrid implicit/explicit model of solvation to investigate the ultrafast excited-state proton transfer (ESPT) reaction of a super photoacid, known as QCy9, in water solution. QCy9 transfers a proton to a water solvent molecule within 100 fs upon the electronic excitation in aqueous solution, and it is the strongest photoacid reported in the literature so far. Because of the ultrafast kinetics, it has been experimentally hypothesized that the ESPT escapes the solvent dynamics control (Huppert et al., J. Photochem. Photobiol. A2014,277, 90). The sampling of the solvent configuration space on the ground electronic state is the first key step toward the simulation of the ESPT event. Therefore, several configurations in the Franck–Condon region, describing an average solvation, were chosen as starting points for the excited-state dynamics. In all cases, the excited-state evolution spontaneously leads to the proton transfer event, whose rate is strongly dependent on the hydrogen bond network around the proton acceptor solvent molecule. Our study revealed that the explicit representation at least of three solvation shells around the proton acceptor molecule is necessary to stabilize the excess proton. Furthermore, the analysis of the solvent molecule motions in proximity of the reaction site suggested that even in the case of the strongest photoacid, the ESPT is actually assisted by the solvation dynamics of the first and second solvation shells of the water accepting molecule.

Exploiting Machine Learning to Efficiently Predict Multidimensional Optical Spectra in Complex Environments

The excited-state dynamics of chromophores in complex environments determine a range of vital biological and energy capture processes. Time-resolved, multidimensional optical spectroscopies provide a key tool to investigate these processes. Although theory has the potential to decode these spectra in terms of the electronic and atomistic dynamics, the need for large numbers of excited-state electronic structure calculations severely limits first-principles predictions of multidimensional optical spectra for chromophores in the condensed phase. Here, we leverage the locality of chromophore excitations to develop machine learning models to predict the excited-state energy gap of chromophores in complex environments for efficiently constructing linear and multidimensional optical spectra. By analyzing the performance of these models, which span a hierarchy of physical approximations, across a range of chromophore-environment interaction strengths, we provide strategies for the construction of machine learning models that greatly accelerate the calculation of multidimensional optical spectra from first principles.

Ab-initio molecular dynamics and hybrid explicit-implicit solvation model for aqueous and nonaqueous solvents: GFP chromophore in water and methanol solution as case study

Solute-solvent interactions are proxies for understanding how the electronic density of a chromophore interacts with the environment in a more exhaustive way. The subtle balance between polarization, electrostatic, and non-bonded interactions need to be accurately described to obtain good agreement between simulations and experiments. First principles approaches providing accurate configurational sampling through molecular dynamics may be a suitable choice to describe solvent effects on solute chemical-physical properties and spectroscopic features, such as optical absorption of dyes. In this context, accurate energy potentials, obtained by hybrid implicit/explicit solvation methods along with employing nonperiodic boundary conditions, are required to represent bulk solvent around a large solute-solvent cluster. In this work, a novel strategy to simulate methanol solutions is proposed combining ab initio molecular dynamics, a hybrid implicit/explicit flexible solvent model, nonperiodic boundary conditions, and time dependent density functional theory. As case study, the robustness of the proposed protocol has been gauged by investigating the microsolvation and electronic absorption of the anionic green fluorescent protein chromophore in methanol and aqueous solution. Satisfactory results are obtained, reproducing the microsolvation layout of the chromophore and, as a consequence, the experimental trends shown by the optical absorption in different solvents.

Nonlinear spectroscopy in the condensed phase: The role of Duschinsky rotations and third order cumulant contributions

First-principles modeling of nonlinear optical spectra in the condensed phase is highly challenging because both environment and vibronic interactions can play a large role in determining spectral shapes and excited state dynamics. Here, we compute two dimensional electronic spectroscopy (2DES) signals based on a cumulant expansion of the energy gap fluctuation operator, with specific focus on analyzing mode mixing effects introduced by the Duschinsky rotation and the role of the third order term in the cumulant expansion for both model and realistic condensed phase systems. We show that for a harmonic model system, the third order cumulant correction captures effects introduced by a mismatch in curvatures of ground and excited state potential energy surfaces, as well as effects of mode mixing. We also demonstrate that 2DES signals can be accurately reconstructed from purely classical correlation functions using quantum correction factors. We then compute nonlinear optical spectra for the Nile red and methylene blue chromophores in solution, assessing the third order cumulant contribution for realistic systems. We show that the third order cumulant correction is strongly dependent on the treatment of the solvent environment, revealing the interplay between environmental polarization and the electronic-vibrational coupling.

Fragment Quantum Mechanical Method for Excited States of Proteins: Development and Application to the Green Fluorescent Protein

Understanding the excited-state properties of luminescent biomolecules is of central importance to their biophysical applications. In this study, we develop the Electrostatically Embedded Generalized Molecular Fractionation with Conjugate Caps (EE-GMFCC) method for quantitatively characterizing properties of covalently bonded systems with localized excitations (i.e., involving a single chromophore), such as fluorescent proteins. The excitation energy, transition dipole moment, and oscillator strength of wild-type Green Fluorescent Protein (wt-GFP) calculated by EE-GMFCC are found to be in excellent agreement with full system time-dependent density functional theory results. We also applied the Polarized Protein-Specific Charge model to wt-GFP, and found that electronic polarization of the protein is critical in stabilizing hydrogen bonding interactions in wt-GFP, which influences its absorption spectrum. The predicted absorption spectra of wt-GFP in the A and B states qualitatively agree with experiment. The fragmentation approach further allows a straightforward per residue decomposition of the excitation which reveals the influence of the protein environment on the absorption spectra of wt-GFP A and B states. Our results demonstrate that the EE-GMFCC method is both accurate and efficient for excited-state property calculations on proteins.

Both Configuration and QM Region Size Matter: Zinc Stability in QM/MM Models of DNA Methyltransferase

Quantum-mechanical/molecular-mechanical (QM/MM) methods are essential to the study of metalloproteins, but the relative importance of sampling and degree of QM treatment in achieving quantitative predictions is poorly understood. We study the relative magnitude of configurational and QM-region sensitivity of energetic and electronic properties in a representative Zn²⁺ metal binding site of a DNA methyltransferase. To quantify property variations, we analyze snapshots extracted from 250 ns of molecular dynamics simulation. To understand the degree of QM-region sensitivity, we perform analysis using QM regions ranging from a minimal 49-atom region consisting only of the Zn²⁺ metal and its four coordinating Cys residues up to a 628-atom QM region that includes residues within 12 Å of the metal center. Over the configurations sampled, we observe that illustrative properties (e.g., rigid Zn²⁺ removal energy) exhibit large fluctuations that are well captured with even minimal QM regions. Nevertheless, for both energetic and electronic properties, we observe a slow approach to asymptotic limits with similarly large changes in absolute values that converge only with larger (ca. 300-atom) QM region sizes. For the smaller QM regions, the electronic description of Zn²⁺ binding is incomplete: the metal binds too tightly and is too stabilized by the strong electrostatic potential of MM point charges, and the Zn-S bond covalency is overestimated. Overall, this work suggests that efficient sampling with QM/MM in small QM regions is an effective method to explore the influence of enzyme structure on target properties. At the same time, accurate descriptions of electronic and energetic properties require a larger QM region than the minimal metal-coordinating residues in order to converge treatment of both metal-local bonding and the overall electrostatic environment.

Systematic microsolvation approach with a cluster-continuum scheme and conformational sampling

Solvation is a notoriously difficult and nagging problem for the rigorous theoretical description of chemistry in the liquid phase. Successes and failures of various approaches ranging from implicit solvation modeling through dielectric continuum embedding and microsolvated quantum chemical modeling to explicit molecular dynamics highlight this situation. Here, we focus on quantum chemical microsolvation and discuss an explicit conformational sampling ansatz to make this approach systematic. For this purpose, we introduce an algorithm for rolling and automated microsolvation of solutes. Our protocol takes conformational sampling and rearrangements in the solvent shell into account. Its reliability is assessed by monitoring the evolution of the spread and average of the observables of interest.

Ionization energies in solution with the QM:QM approach

We discuss a fragment-based QM:QM scheme as a practical way to access the energetics of vertical electronic processes in the condensed phase. In the QM:QM scheme, we decompose the large molecular system into small fragments, which interact solely electrostatically. The energies of the fragments are calculated in a self-consistent field generated by the other fragments and the total energy of the system is calculated as a sum of the fragment energies. We show on two test cases (cytosine and a sodium cation) that the method allows one to accurately simulate the shift of vertical ionization energies (VIE) while going from the gas phase to the bulk. For both examples, the predicted solvent shifts and peak widths estimated at the DFT level agree well with the experimental observations. We argue that the QM:QM approach is more suitable than either an electrostatic embedding based QM/MM approach, a full quantum description at the DFT level with a generally used functional or a combination of both. We also discuss the potential scope of the applicability for other electronic processes such as Auger decay.

Optical spectra in the condensed phase: Capturing anharmonic and vibronic features using dynamic and static approaches

Simulating optical spectra in the condensed phase remains a challenge for theory due to the need to capture spectral signatures arising from anharmonicity and dynamical effects, such as vibronic progressions and asymmetry. As such, numerous simulation methods have been developed that invoke different approximations and vary in their ability to capture different physical regimes. Here, we use several models of chromophores in the condensed phase and ab initio molecular dynamics simulations to rigorously assess the applicability of methods to simulate optical absorption spectra. Specifically, we focus on the ensemble scheme, which can address anharmonic potential energy surfaces but relies on the applicability of extreme nuclear-electronic time scale separation; the Franck-Condon method, which includes dynamical effects but generally only at the harmonic level; and the recently introduced ensemble zero-temperature Franck-Condon approach, which straddles these limits. We also devote particular attention to the performance of methods derived from a cumulant expansion of the energy gap fluctuations and test the ability to approximate the requisite time correlation functions using classical dynamics with quantum correction factors. These results provide insights as to when these methods are applicable and able to capture the features of condensed phase spectra qualitatively and, in some cases, quantitatively across a range of regimes.

Revealing quantum mechanical effects in enzyme catalysis with large-scale electronic structure simulation

Enzymes have evolved to facilitate challenging reactions at ambient conditions with specificity seldom matched by other catalysts. Computational modeling provides valuable insight into catalytic mechanism, and the large size of enzymes mandates multi-scale, quantum mechanical-molecular mechanical (QM/MM) simulations. Although QM/MM plays an essential role in balancing simulation cost to enable sampling with full QM treatment needed to understand electronic structure in enzyme active sites, the relative importance of these two strategies for understanding enzyme mechanism is not well known. We explore challenges in QM/MM for studying the reactivity and stability of three diverse enzymes: i) Mg2+-dependent catechol O-methyltransferase (COMT), ii) radical enzyme choline trimethylamine lyase (CutC), and iii) DNA methyltransferase (DNMT1), which has structural Zn2+ binding sites. In COMT, strong non-covalent interactions lead to long range coupling of electronic structure properties across the active site, but the more isolated nature of the metallocofactor in DNMT1 leads to faster convergence of some properties. We quantify these effects in COMT by computing covariance matrices of by-residue electronic structure properties during dynamics and along the reaction coordinate. In CutC, we observe spontaneous bond cleavage following initiation events, highlighting the importance of sampling and dynamics. We use electronic structure analysis to quantify the relative importance of CHO and OHO non-covalent interactions in imparting reactivity. These three diverse cases enable us to provide some general recommendations regarding QM/MM simulation of enzymes.

Computational Understanding of the Selectivities in Metalloenzymes

Metalloenzymes catalyze many different types of biological reactions with high efficiency and remarkable selectivity. The quantum chemical cluster approach and the combined quantum mechanics/molecular mechanics methods have proven very successful in the elucidation of the reaction mechanism and rationalization of selectivities in enzymes. In this review, recent progress in the computational understanding of various selectivities including chemoselectivity, regioselectivity, and stereoselectivity, in metalloenzymes, is discussed.

Assessing Configurational Sampling in the Quantum Mechanics/Molecular Mechanics Calculation of Temoporfin Absorption Spectrum and Triplet Density of States

The absorption properties of Temoporfin, a second-generation photosensitizer employed in photodynamic therapy, are calculated with an electrostatic-embedding quantum mechanics/molecular mechanics (QM/MM) scheme in methanol. The suitability of several ensembles of geometries generated by different sampling techniques, namely classical-molecular-dynamics (MD) and QM/MM-MD thermal sampling, Wigner quantum sampling and a hybrid protocol, which combines the thermal and quantum approaches, is assessed. It is found that a QM description of the chromophore during the sampling is needed in order to achieve a good agreement with respect to the experimental spectrum. Such a good agreement is obtained with both QM/MM-MD and Wigner samplings, demonstrating that a proper description of the anharmonic motions of the chromophore is not relevant in the computation of the absorption properties. In addition, it is also found that solvent organization is a rather fast process and a long sampling is not required. Finally, it is also demonstrated that the same exchange-correlation functional should be employed in the sampling and in the computation of the excited states properties to avoid unphysical triplet states with relative energies close or below 0 eV.

Do Better Quality Embedding Potentials Accelerate the Convergence of QM/MM Models? The Case of Solvated Acid Clusters

This study examines whether the use of more accurate embedding potentials improves the convergence of quantum mechanics/molecular mechanics (QM/MM) models with respect to the size of the QM region. In conjunction with density functional theory calculations using the ωB97X-D functional, various embedding potentials including the TIP3P water model, the effective fragment potential (EFP), and semi-empirical methods (PM6, PM7, and DFTB) were used to simulate the deprotonation energies of solvated acid clusters. The calculations were performed on solvated neutral (HA) and cationic (HB⁺) acids clusters containing 160 and 480 water molecules using configurations sampled from molecular dynamics simulations. Consistently, the ωB97X-D/EFP model performed the best when using a minimal QM region size. The performance for the other potentials appears to be highly sensitive to the charge character of the acid/base pair. Neutral acids display the expected trend that semi-empirical methods generally perform better than TIP3P; however, an opposite trend was observed for the cationic acids. Additionally, electronic embedding provided an improvement over mechanical embedding for the cationic systems, but not the neutral acids. For the best performing ωB97X-D/EFP model, a QM region containing about 6% of the total number of solvent molecules is needed to approach within 10 kJ mol^-1 of the pure QM result if the QM region was chosen based on the distance from the reaction centre.

Large-scale QM/MM free energy simulations of enzyme catalysis reveal the influence of charge transfer

Hybrid quantum mechanical-molecular mechanical (QM/MM) simulations provide key insights into enzyme structure-function relationships. Numerous studies have demonstrated that large QM regions are needed to systematically converge ground state, zero temperature properties with electrostatic embedding QM/MM. However, it is not well known if ab initio QM/MM free energy simulations have this same dependence, in part due to the hundreds of thousands of energy evaluations required for free energy estimations that in turn limit QM region size. Here, we leverage recent advances in electronic structure efficiency and accuracy to carry out range-separated hybrid density functional theory free energy simulations in a representative methyltransferase. By studying 200 ps of ab initio QM/MM dynamics for each of five QM regions from minimal (64 atoms) to one-sixth of the protein (544 atoms), we identify critical differences between large and small QM region QM/MM in charge transfer between substrates and active site residues as well as in geometric structure and dynamics that coincide with differences in predicted free energy barriers. Distinct geometric and electronic structure features in the largest QM region indicate that important aspects of enzymatic rate enhancement in methyltransferases are identified with large-scale electronic structure.

Unraveling electronic absorption spectra using nuclear quantum effects: Photoactive yellow protein and green fluorescent protein chromophores in water

Many physical phenomena must be accounted for to accurately model solution-phase optical spectral line shapes, from the sampling of chromophore-solvent configurations to the electronic-vibrational transitions leading to vibronic fine structure. Here we thoroughly explore the role of nuclear quantum effects, direct and indirect solvent effects, and vibronic effects in the computation of the optical spectrum of the aqueously solvated anionic chromophores of green fluorescent protein and photoactive yellow protein. By analyzing the chromophore and solvent configurations, the distributions of vertical excitation energies, the absorption spectra computed within the ensemble approach, and the absorption spectra computed within the ensemble plus zero-temperature Franck-Condon approach, we show how solvent, nuclear quantum effects, and vibronic transitions alter the optical absorption spectra. We find that including nuclear quantum effects in the sampling of chromophore-solvent configurations using ab initio path integral molecular dynamics simulations leads to improved spectral shapes through three mechanisms. The three mechanisms that lead to line shape broadening and a better description of the high-energy tail are softening of heavy atom bonds in the chromophore that couple to the optically bright state, widening the distribution of vertical excitation energies from more diverse solvation environments, and redistributing spectral weight from the 0-0 vibronic transition to higher energy vibronic transitions when computing the Franck-Condon spectrum in a frozen solvent pocket. The absorption spectra computed using the combined ensemble plus zero-temperature Franck-Condon approach yield significant improvements in spectral shape and width compared to the spectra computed with the ensemble approach. Using the combined approach with configurations sampled from path integral molecular dynamics trajectories presents a significant step forward in accurately modeling the absorption spectra of aqueously solvated chromophores.

Combining the ensemble and Franck-Condon approaches for calculating spectral shapes of molecules in solution

The correct treatment of vibronic effects is vital for the modeling of absorption spectra of many solvated dyes. Vibronic spectra for small dyes in solution can be easily computed within the Franck-Condon approximation using an implicit solvent model. However, implicit solvent models neglect specific solute-solvent interactions on the electronic excited state. On the other hand, a straightforward way to account for solute-solvent interactions and temperature-dependent broadening is by computing vertical excitation energies obtained from an ensemble of solute-solvent conformations. Ensemble approaches usually do not account for vibronic transitions and thus often produce spectral shapes in poor agreement with experiment. We address these shortcomings by combining zero-temperature vibronic fine structure with vertical excitations computed for a room-temperature ensemble of solute-solvent configurations. In this combined approach, all temperature-dependent broadening is treated classically through the sampling of configurations and quantum mechanical vibronic contributions are included as a zero-temperature correction to each vertical transition. In our calculation of the vertical excitations, significant regions of the solvent environment are treated fully quantum mechanically to account for solute-solvent polarization and charge-transfer. For the Franck-Condon calculations, a small amount of frozen explicit solvent is considered in order to capture solvent effects on the vibronic shape function. We test the proposed method by comparing calculated and experimental absorption spectra of Nile red and the green fluorescent protein chromophore in polar and non-polar solvents. For systems with strong solute-solvent interactions, the combined approach yields significant improvements over the ensemble approach. For systems with weak to moderate solute-solvent interactions, both the high-energy vibronic tail and the width of the spectra are in excellent agreement with experiments.