Electronic Absorption Spectra from MM and ab initio QM/MM Molecular Dynamics: Environmental Effects on the Absorption Spectrum of Photoactive Yellow Protein

Effect of Leaving Centrosymmetric Character on Spectral Properties in Mono-, Bi-, and Triphotonic Absorption Spectroscopies

Numerical simulations of the absorption bands of photoswitch E-o-tetrafluoroazobenzene in DMSO solution under one-, two-, and three-photon absorption conditions combined with the analysis of the behavior of transition probability under distortion of planarity reveal many similarities between the mono- and triphoton spectroscopic behaviors with a two-photon spectrum being set apart. The position of the absorption peak for the studied nπ* and ππ* transitions appears shifted to lower energies (longer wavelengths) than the conventional estimate based on vertical excitation from the ground-state potential energy minimum.

Modeling UV/Vis Absorption Spectra of Food Colorants in Solution: Anthocyanins and Curcumin as Case Studies

We present a comprehensive computational study of UV/Vis absorption spectra of significant food colorants, specifically anthocyanins and curcumin tautomers, dissolved in polar protic solvents, namely water and ethanol. The absorption spectra are simulated using two fully polarizable quantum mechanical (QM)/molecular mechanics (MM) models based on the fluctuating charge (FQ) and fluctuating charge and dipoles (FQFμ) force fields. To accurately capture the dynamical aspects of the solvation phenomenon, atomistic approaches are combined with configurational sampling obtained through classical molecular dynamics (MD) simulations. The calculated QM/FQ and QM/FQFμ spectra are then compared with experiments. Our findings demonstrate that a precise reproduction of the UV/Vis spectra of the studied pigments can be achieved by adequately accounting for configurational sampling, polarization effects, and hydrogen bonding interactions.

GPU-accelerated on-the-fly nonadiabatic semiclassical dynamics

GPU-accelerated on-the-fly nonadiabatic dynamics is enabled by interfacing the linearized semiclassical dynamics approach with the TeraChem electronic structure program. We describe the computational workflow of the "PySCES" code interface, a Python code for semiclassical dynamics with on-the-fly electronic structure, including parallelization over multiple GPU nodes. We showcase the abilities of this code and present timings for two benchmark systems: fulvene solvated in acetonitrile and a charge transfer system in which a photoexcited zinc-phthalocyanine donor transfers charge to a fullerene acceptor through multiple electronic states on an ultrafast timescale. Our implementation paves the way for an efficient semiclassical approach to model the nonadiabatic excited state dynamics of complex molecules, materials, and condensed phase systems.

Calculating absorption and fluorescence spectra for chromophores in solution with ensemble Franck-Condon methods

Accurately modeling absorption and fluorescence spectra for molecules in solution poses a challenge due to the need to incorporate both vibronic and environmental effects, as well as the necessity of accurate excited state electronic structure calculations. Nuclear ensemble approaches capture explicit environmental effects, Franck-Condon methods capture vibronic effects, and recently introduced ensemble-Franck-Condon approaches combine the advantages of both methods. In this study, we present and analyze simulated absorption and fluorescence spectra generated with combined ensemble-Franck-Condon approaches for three chromophore-solvent systems and compare them to standard ensemble and Franck-Condon spectra, as well as to the experiment. Employing configurations obtained from ground and excited state ab initio molecular dynamics, three combined ensemble-Franck-Condon approaches are directly compared to each other to assess the accuracy and relative computational time. We find that the approach employing an average finite-temperature Franck-Condon line shape generates spectra nearly identical to the direct summation of an ensemble of Franck-Condon spectra at one-fourth of the computational cost. We analyze how the spectral simulation method, as well as the level of electronic structure theory, affects spectral line shapes and associated Stokes shifts for 7-nitrobenz-2-oxa-1,3-diazol-4-yl and Nile red in dimethyl sulfoxide and 7-methoxy coumarin-4-acetic acid in methanol. For the first time, our studies show the capability of combined ensemble-Franck-Condon methods for both absorption and fluorescence spectroscopy and provide a powerful tool for simulating linear optical spectra.

MiMiC: A high-performance framework for multiscale molecular dynamics simulations

MiMiC is a framework for performing multiscale simulations in which loosely coupled external programs describe individual subsystems at different resolutions and levels of theory. To make it highly efficient and flexible, we adopt an interoperable approach based on a multiple-program multiple-data (MPMD) paradigm, serving as an intermediary responsible for fast data exchange and interactions between the subsystems. The main goal of MiMiC is to avoid interfering with the underlying parallelization of the external programs, including the operability on hybrid architectures (e.g., CPU/GPU), and keep their setup and execution as close as possible to the original. At the moment, MiMiC offers an efficient implementation of electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) that has demonstrated unprecedented parallel scaling in simulations of large biomolecules using CPMD and GROMACS as QM and MM engines, respectively. However, as it is designed for high flexibility with general multiscale models in mind, it can be straightforwardly extended beyond QM/MM. In this article, we illustrate the software design and the features of the framework, which make it a compelling choice for multiscale simulations in the upcoming era of exascale high-performance computing.

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.

GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

The development of the QM/MM interface and its application for the on-the-fly QM/MM nonadiabatic dynamics in JADE package: Theory, implementation, and applications

Understanding the nonadiabatic dynamics of complex systems is a challenging task in computational photochemistry. Herein, we present an efficient and user-friendly quantum mechanics/molecular mechanics (QM/MM) interface to run on-the-fly nonadiabatic dynamics. Currently, this interface consists of an independent set of codes designed for general-purpose use. Herein, we demonstrate the ability and feasibility of the QM/MM interface by integrating it with our long-term developed JADE package. Tailored to handle nonadiabatic processes in various complex systems, especially condensed phases and protein environments, we delve into the theories, implementations, and applications of on-the-fly QM/MM nonadiabatic dynamics. The QM/MM approach is established within the framework of the additive QM/MM scheme, employing electrostatic embedding, link-atom inclusion, and charge-redistribution schemes to treat the QM/MM boundary. Trajectory surface-hopping dynamics are facilitated using the fewest switches algorithm, encompassing classical and quantum treatments for nuclear and electronic motions, respectively. Finally, we report simulations of nonadiabatic dynamics for two typical systems: azomethane in water and the retinal chromophore PSB3 in a protein environment. Our results not only illustrate the power of the QM/MM program but also reveal the important roles of environmental factors in nonadiabatic processes.

Axial H-Bonding Solvent Controls Inhomogeneous Spectral Broadening, While Peripheral H-Bonding Solvent Controls Vibronic Broadening: Cresyl Violet in Methanol

The dynamics of the nuclei of both a chromophore and its condensed-phase environment control many spectral features, including the vibronic and inhomogeneous broadening present in spectral line shapes. For the cresyl violet chromophore in methanol, we here analyze and isolate the effect of specific chromophore-solvent interactions on simulated spectral densities, reorganization energies, and linear absorption spectra. Employing both chromophore and its condensed-phase environment control many spectral features, including the vibronic and inhomogeneous broadening present in spectral line shapes. For the cresyl violet chromophore in methanol, we here analyze and isolate the effect of specific chromophore-solvent interactions on simulated spectral densities, reorganization energies, and linear absorption spectra. Employing both force field and ab initio molecular dynamics trajectories along with the inclusion of only certain solvent molecules in the excited-state calculations, we determine that the methanol molecules axial to the chromophore are responsible for the majority of inhomogeneous broadening, with a single methanol molecule that forms an axial hydrogen bond dominating the response. The strong peripheral hydrogen bonds do not contribute to spectral broadening, as they are very stable throughout the dynamics and do not lead to increased energy-gap fluctuations. We also find that treating the strong peripheral hydrogen bonds as molecular mechanical point charges during the molecular dynamics simulation underestimates the vibronic coupling. Including these peripheral hydrogen bonding methanol molecules in the quantum-mechanical region in a geometry optimization increases the vibronic coupling, suggesting that a more advanced treatment of these strongly interacting solvent molecules during the molecular dynamics trajectory may be necessary to capture the full vibronic spectral broadening.

Multiscale biomolecular simulations in the exascale era

The complexity of biological systems and processes, spanning molecular to macroscopic scales, necessitates the use of multiscale simulations to get a comprehensive understanding. Quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations are crucial for capturing processes beyond the reach of classical MD simulations. The advent of exascale computing offers unprecedented opportunities for scientific exploration, not least within life sciences, where simulations are essential to unravel intricate molecular mechanisms underlying biological processes. However, leveraging the immense computational power of exascale computing requires innovative algorithms and software designs. In this context, we discuss the current status and future prospects of multiscale biomolecular simulations on exascale supercomputers with a focus on QM/MM MD. We highlight our own efforts in developing a versatile and high-performance multiscale simulation framework with the aim of efficient utilization of state-of-the-art supercomputers. We showcase its application in uncovering complex biological mechanisms and its potential for leveraging exascale computing.

Computational Spectroscopy of Aqueous Solutions: The Underlying Role of Conformational Sampling

Fully atomistic multiscale polarizable quantum mechanics (QM)/molecular mechanics (MM) approaches, combined with techniques to sample the solute-solvent phase space, constitute the most accurate method to compute spectral signals in aqueous solution. Conventional sampling strategies, such as classical molecular dynamics (MD), may encounter drawbacks when the conformational space is particularly complex, and transition barriers between conformers are high. This can lead to inaccurate sampling, which can potentially impact the accuracy of spectral calculations. For this reason, in this work, we compare classical MD with enhanced sampling techniques, i.e., replica exchange MD and metadynamics. In particular, we show how the different sampling techniques affect computed UV, electronic circular dichroism, nuclear magnetic resonance shielding, and optical rotatory dispersion of N-acetylproline-amide in aqueous solution. Such a system is a model peptide characterized by complex conformational variability. Calculated values suggest that spectral properties are influenced by solute conformers, relative population, and solvent effects; therefore, particular care needs to be paid for when choosing the sampling technique.

Exploring the Mechanisms behind Non-aromatic Fluorescence with the Density Functional Tight Binding Method

Recent experimental findings reveal nonconventional fluorescence emission in biological systems devoid of conjugated bonds or aromatic compounds, termed non-aromatic fluorescence (NAF). This phenomenon is exclusive to aggregated or solid states and remains absent in monomeric solutions. Previous studies focused on small model systems in vacuum show that the carbonyl stretching mode along with strong interaction of short hydrogen bonds (SHBs) remains the primary vibrational mode explaining NAF in these systems. In order to simulate larger model systems taking into account the effects of the surrounding environment, in this work we propose using the density functional tight-binding (DFTB) method in combination with non-adiabatic molecular dynamics (NAMD) and the mixed quantum/molecular mechanics (QM/MM) approach. We investigate the mechanism behind NAF in the crystal structure of l-pyroglutamine-ammonium, comparing it with the related nonfluorescent amino acid l-glutamine. Our results extend our previous findings to more realistic systems, demonstrating the efficiency and robustness of the proposed DFTB method in the context of NAMD in biological systems. Furthermore, due to its inherent low computational cost, this method allows for a better sampling of the nonradiative events at the conical intersection which is crucial for a complete understanding of this phenomenon. Beyond contributing to the ongoing exploration of NAF, this work paves the way for future application of this method in more complex biological systems such as amyloid aggregates, biomaterials, and non-aromatic proteins.

The Role of Aqueous Solvation on the Intersystem Crossing of Nitrophenols

The photochemistry of nitrophenols is a source of smog as nitrous acid is formed from their photolysis. Nevertheless, computational studies of the photochemistry of these widespread toxic molecules are scarce. In this work, the initial photodeactivation of ortho-nitrophenol and para-nitrophenol is modeled, both in gas phase and in aqueous solution to simulate atmospheric and aerosol environments. A large number of excited states, six for ortho-nitrophenol and 11 for para-nitrophenol, have been included and were all populated during the decay. Moreover, periodic time-dependent density functional theory (TDDFT) is used for both the explicitly included solvent and the solute. A comparison to periodic QM/MM (TDDFT/MM), with electrostatic embedding, is made, showing notable differences between the decays of solvated nitrophenols simulated with QM/MM and full (TD)DFT. A reduced intersystem crossing in aqueous solution could be observed thanks to the surface hopping approach using explicit, periodic TDDFT solvation including spin-orbit couplings.

Taming the third order cumulant approximation to linear optical spectroscopy

The second order cumulant method offers a promising pathway to predicting optical properties in condensed phase systems. It allows for the computation of linear absorption spectra from excitation energy fluctuations sampled along molecular dynamics (MD) trajectories, fully accounting for vibronic effects, direct solute-solvent interactions, and environmental polarization effects. However, the second order cumulant approximation only guarantees accurate line shapes for energy gap fluctuations obeying Gaussian statistics. A third order correction has recently been derived but often yields unphysical spectra or divergent line shapes for moderately non-Gaussian fluctuations due to the neglect of higher order terms in the cumulant expansion. In this work, we develop a corrected cumulant approach, where the collective effect of neglected higher order contributions is approximately accounted for through a dampening factor applied to the third order cumulant term. We show that this dampening factor can be expressed as a function of the skewness and kurtosis of energy gap fluctuations and can be parameterized from a large set of randomly sampled model Hamiltonians for which exact spectral line shapes are known. This approach is shown to systematically remove unphysical contributions in the form of negative absorbances from cumulant spectra in both model Hamiltonians and condensed phase systems sampled from MD and dramatically improves over the second order cumulant method in describing systems exhibiting Duschinsky mode mixing effects. We successfully apply the approach to the coumarin-153 dye in toluene, obtaining excellent agreement with experiment.

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.

Embedding Nonrigid Solutes in an Averaged Environment: A Case Study on Rhodopsins

Many simulation methods concerning solvated molecules are based on the assumption that the solvated species and the solvent can be characterized by some representative structures of the solute and some embedding potential corresponding to this structure. While the averaging of the solvent configurations to obtain an embedding potential has been studied in great detail, this hinges on a single solute structure representation. This assumption is re-examined and generalized for conformationally flexible solutes and tested on 4 nonrigid systems. In this generalized approach, the solute is characterized by a set of representative structures and the corresponding embedding potentials. The representative structures are identified by means of subdividing the statistical ensemble, which in this work is generated by a constant-temperature molecular dynamics simulation. The embedding potential defined in the Frozen-Density Embedding Theory is used to characterize the average effect of the solvent in each subensemble. The numerical examples concern the vertical excitation energies of protonated retinal Schiff bases in protein environments. It is comprehensively shown that subensemble averaging leads to huge computational savings compared with explicit averaging of the excitation energies in the whole ensemble while introducing only minor errors in the case of the systems examined.

An Expedited Route to Optical and Electronic Properties at Finite Temperature via Unsupervised Learning

Electronic properties and absorption spectra are the grounds to investigate molecular electronic states and their interactions with the environment. Modeling and computations are required for the molecular understanding and design strategies of photo-active materials and sensors. However, the interpretation of such properties demands expensive computations and dealing with the interplay of electronic excited states with the conformational freedom of the chromophores in complex matrices (i.e., solvents, biomolecules, crystals) at finite temperature. Computational protocols combining time dependent density functional theory and ab initio molecular dynamics (MD) have become very powerful in this field, although they require still a large number of computations for a detailed reproduction of electronic properties, such as band shapes. Besides the ongoing research in more traditional computational chemistry fields, data analysis and machine learning methods have been increasingly employed as complementary approaches for efficient data exploration, prediction and model development, starting from the data resulting from MD simulations and electronic structure calculations. In this work, dataset reduction capabilities by unsupervised clustering techniques applied to MD trajectories are proposed and tested for the ab initio modeling of electronic absorption spectra of two challenging case studies: a non-covalent charge-transfer dimer and a ruthenium complex in solution at room temperature. The K-medoids clustering technique is applied and is proven to be able to reduce by ∼100 times the total cost of excited state calculations on an MD sampling with no loss in the accuracy and it also provides an easier understanding of the representative structures (medoids) to be analyzed on the molecular scale.

Cost-Effective Simulations of Vibrationally-Resolved Absorption Spectra of Fluorophores with Machine-Learning-Based Inhomogeneous Broadening

The results of electronic and vibrational structure simulations are an invaluable support for interpreting experimental absorption/emission spectra, which stimulates the development of reliable and cost-effective computational protocols. In this work, we contribute to these efforts and propose an efficient first-principle protocol for simulating vibrationally-resolved absorption spectra, including nonempirical estimations of the inhomogeneous broadening. To this end, we analyze three key aspects: (i) a metric-based selection of density functional approximation (DFA) so to benefit from the computational efficiency of time-dependent density function theory (TD-DFT) while safeguarding the accuracy of the vibrationally-resolved spectra, (ii) an assessment of two vibrational structure schemes (vertical gradient and adiabatic Hessian) to compute the Franck-Condon factors, and (iii) the use of machine learning to speed up nonempirical estimations of the inhomogeneous broadening. In more detail, we predict the absorption band shapes for a set of 20 medium-sized fluorescent dyes, focusing on the bright ππ^★ S₀ → S₁ transition and using experimental results as references. We demonstrate that, for the studied 20-dye set which includes structures with large structural variability, the preselection of DFAs based on an easily accessible metric ensures accurate band shapes with respect to the reference approach and that range-separated functionals show the best performance when combined with the vertical gradient model. As far as band widths are concerned, we propose a new machine-learning-based approach for determining the inhomogeneous broadening induced by the solvent microenvironment. This approach is shown to be very robust offering inhomogeneous broadenings with errors as small as 2 cm^-1 with respect to genuine electronic-structure calculations, with a total CPU time reduced by 98%.

Spatial Decay and Limits of Quantum Solute-Solvent Interactions

Molecular excitations in the liquid-phase environment are renormalized by the surrounding solvent molecules. Herein, we employ the GW approximation to investigate the solvation effects on the ionization energy of phenol in various solvent environments. The electronic effects differ by up to 0.4 eV among the five investigated solvents. This difference depends on both the macroscopic solvent polarizability and the spatial decay of the solvation effects. The latter is probed by separating the electronic subspace and the GW correlation self-energy into fragments. The fragment correlation energy decays with increasing intermolecular distance and vanishes at ∼9 Å, and this pattern is independent of the type of solvent environment. The 9 Å cutoff defines an effective interacting volume within which the ionization energy shift per solvent molecule is proportional to the macroscopic solvent polarizability. Finally, we propose a simple model for computing the ionization energies of molecules in an arbitrary solvent environment.

Quantum Mechanics/Molecular Mechanics Simulations on NVIDIA and AMD Graphics Processing Units

We have ported and optimized the graphics processing unit (GPU)-accelerated QUICK and AMBER-based ab initio quantum mechanics/molecular mechanics (QM/MM) implementation on AMD GPUs. This encompasses the entire Fock matrix build and force calculation in QUICK including one-electron integrals, two-electron repulsion integrals, exchange-correlation quadrature, and linear algebra operations. General performance improvements to the QUICK GPU code are also presented. Benchmarks carried out on NVIDIA V100 and AMD MI100 cards display similar performance on both hardware for standalone HF/DFT calculations with QUICK and QM/MM molecular dynamics simulations with QUICK/AMBER. Furthermore, with respect to the QUICK/AMBER release version 21, significant speedups are observed for QM/MM molecular dynamics simulations. This significantly increases the range of scientific problems that can be addressed with open-source QM/MM software on state-of-the-art computer hardware.

TeraChem protocol buffers (TCPB): Accelerating QM and QM/MM simulations with a client-server model

The routine use of electronic structures in many chemical simulation applications calls for efficient and easy ways to access electronic structure programs. We describe how the graphics processing unit (GPU) accelerated electronic structure program TeraChem can be set up as an electronic structure server, to be easily accessed by third-party client programs. We exploit Google's protocol buffer framework for data serialization and communication. The client interface, called TeraChem protocol buffers (TCPB), has been designed for ease of use and compatibility with multiple programming languages, such as C++, Fortran, and Python. To demonstrate the ease of coupling third-party programs with electronic structures using TCPB, we have incorporated the TCPB client into Amber for quantum mechanics/molecular mechanics (QM/MM) simulations. The TCPB interface saves time with GPU initialization and I/O operations, achieving a speedup of more than 2× compared to a prior file-based implementation for a QM region with ∼250 basis functions. We demonstrate the practical application of TCPB by computing the free energy profile of p-hydroxybenzylidene-2,3-dimethylimidazolinone (p-HBDI^-)-a model chromophore in green fluorescent proteins-on the first excited singlet state using Hamiltonian replica exchange for enhanced sampling. All calculations in this work have been performed with the non-commercial freely-available version of TeraChem, which is sufficient for many QM region sizes in common use.

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.

Multiple Facets of Modeling Electronic Absorption Spectra of Systems in Solution

In this Perspective, we outline the essential physicochemical aspects that need to be considered when building a reliable approach to describe absorption properties of solvated systems. In particular, we focus on how to properly model the complexity of the solvation phenomenon, arising from dynamical aspects and specific, strong solute-solvent interactions. To this end, conformational and configurational sampling techniques, such as Molecular Dynamics, have to be coupled to accurate fully atomistic Quantum Mechanical/Molecular Mechanics (QM/MM) methodologies. By exploiting different illustrative applications, we show that an effective reproduction of experimental spectral signals can be achieved by delicately balancing exhaustive sampling, hydrogen bonding, mutual polarization, and nonelectrostatic effects.

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.

Spectroscopic and QM/MM studies of the Cu(I) binding site of the plant ethylene receptor ETR1

Herein, we present, to our knowledge, the first spectroscopic characterization of the Cu(I) active site of the plant ethylene receptor ETR1. The x-ray absorption (XAS) and extended x-ray absorption fine structure (EXAFS) spectroscopies presented here establish that ETR1 has a low-coordinate Cu(I) site. The EXAFS resolves a mixed first coordination sphere of N/O and S scatterers at distances consistent with potential histidine and cysteine residues. This finding agrees with the coordination of residues C65 and H69 to the Cu(I) site, which are critical for ethylene activity and well conserved. Furthermore, the Cu K-edge XAS and EXAFS of ETR1 exhibit spectroscopic changes upon addition of ethylene that are attributed to modifications in the Cu(I) coordination environment, suggestive of ethylene binding. Results from umbrella sampling simulations of the proposed ethylene binding helix of ETR1 at a mixed quantum mechanics/molecular mechanics level agree with the EXAFS fit distance changes upon ethylene binding, particularly in the increase of the distance between H69 and Cu(I), and yield binding energetics comparable with experimental dissociation constants. The observed changes in the copper coordination environment might be the triggering signal for the transmission of the ethylene response.

Excited-State Barrier Controls E → Z Photoisomerization in p-Hydroxycinnamate Biochromophores

Molecules based on the deprotonated p-hydroxycinnamate moiety are widespread in nature, including serving as UV filters in the leaves of plants and as the biochromophore in photoactive yellow protein. The photophysical behavior of these chromophores is centered around a rapid E → Z photoisomerization by passage through a conical intersection seam. Here, we use photoisomerization and photodissociation action spectroscopies with deprotonated 4-hydroxybenzal acetone (pCK^-) to characterize a wavelength-dependent bifurcation between electron autodetachment (spontaneous ejection of an electron from the S₁ state because it is situated in the detachment continuum) and E → Z photoisomerization. While autodetachment occurs across the entire S₁(ππ*) band (370-480 nm), E → Z photoisomerization occurs only over a blue portion of the band (370-430 nm). No E → Z photoisomerization is observed when the ketone functional group in pCK^- is replaced with an ester or carboxylic acid. The wavelength-dependent bifurcation is consistent with potential energy surface calculations showing that a barrier separates the Franck-Condon region from the E → Z isomerizing conical intersection. The barrier height, which is substantially higher in the gas phase than in solution, depends on the functional group and governs whether E → Z photoisomerization occurs more rapidly than autodetachment.

Quantum chemistry study of the multiphoton absorption in enhanced green fluorescent protein at the single amino acid residue level

The chromophore (CRO) of fluorescent proteins (FPs) is embedded in a complex environment that is a source of specific interactions with the CRO. Understanding how these interactions influence FPs spectral properties is important for a directed design of novel markers with desired characteristics. In this work, we apply computational chemistry methods to gain insight into one-, two- and three-photon absorption (1PA, 2PA, 3PA) tuning in enhanced green fluorescent protein (EGFP). To achieve this goal, we built EGFP models differing in: i) number and position of hydrogen-bonds (h-bonds) donors to the CRO and ii) the electric field, as approximated by polarizable force  field, acting on the CRO. We  find that h-bonding to the CRO's phenolate oxygen results in stronger one- and multiphoton absorption. The brighter absorption can be also achieved by creating more positive electric field near the CRO's phenolate moiety. Interestingly, while individual CRO-environment h-bonds usually enhance 1PA and 2PA, it takes a few h-bond donors to enhance 3PA. Clearly, response of the absorption intensity to many-body effects depends on the excitation mechanism.  We further employ symmetry-adapted perturbation theory (SAPT) to reveal excellent (2PA) and good (3PA) correlation of multiphoton intensity with electrostatic and induction interaction energies. This points to importance of accounting for mutual CRO-environment polarization in quantitative calculations of absorption spectra in FPs.

Multistructural microiteration combined with QM/MM-ONIOM electrostatic embedding

Multistructural microiteration (MSM) is a method to take account of contributions of multiple surrounding structures in a geometrical optimization or reaction path calculation using the quantum mechanics/molecular mechanics (QM/MM) ONIOM method. In this study, we combined MSM with the electrostatic embedding (EE) scheme of the QM/MM-ONIOM method by extending its original formulation for mechanical embedding (ME). MSM-EE takes account of the polarization in the QM region induced by point charges assigned to atoms in the multiple surrounding structures, where the point charges are scaled by the weight factor of each surrounding structure determined through MSM. The performance of MSM-EE was compared with that of the other methods, i.e., ONIOM-ME, ONIOM-EE, and MSM-ME, by applying them to three chemical processes: (1) chorismate-to-prephenate transformation in aqueous solution, (2) the same transformation as (1) in an enzyme, and (3) hydroxylation in p-hydroxybenzoate hydroxylase. These numerical tests of MSM-EE yielded barriers and reaction energies close to experimental values with computational costs comparable to those of the other three methods.

Transition-Based Constrained DFT for the Robust and Reliable Treatment of Excitations in Supramolecular Systems

Despite the variety of available computational approaches, state-of-the-art methods for calculating excitation energies, such as time-dependent density functional theory (TDDFT), are computationally demanding and thus limited to moderate system sizes. Here, we introduce a new variation of constrained DFT (CDFT), wherein the constraint corresponds to a particular transition (T), or a combination of transitions, between occupied and virtual orbitals, rather than a region of the simulation space as in traditional CDFT. We compare T-CDFT with TDDFT and ΔSCF results for the low-lying excited states (S1 and T1) of a set of gas-phase acene molecules and OLED emitters and with reference results from the literature. At the PBE level of theory, T-CDFT outperforms ΔSCF for both classes of molecules, while also proving to be more robust. For the local excitations seen in the acenes, T-CDFT and TDDFT perform equally well. For the charge transfer (CT)-like excitations seen in the OLED molecules, T-CDFT also performs well, in contrast to the severe energy underestimation seen with TDDFT. In other words, T-CDFT is equally applicable to both local excitations and CT states, providing more reliable excitation energies at a much lower computational cost than TDDFT cost. T-CDFT is designed for large systems and has been implemented in the linear-scaling BigDFT code. It is therefore ideally suited for exploring the effects of explicit environments on excitation energies, paving the way for future simulations of excited states in complex realistic morphologies, such as those which occur in OLED materials.

Factors That Determine the Variation of Equilibrium and Kinetic Properties of QM/MM Enzyme Simulations: QM Region, Conformation, and Boundary Condition

To analyze the impact of various technical details on the results of quantum mechanical (QM)/molecular mechanical (MM) enzyme simulations, including the QM region size, catechol-O-methyltransferase (COMT) is studied as a model system using an approximate QM/MM method (DFTB3/CHARMM). The results show that key equilibrium and kinetic properties for methyl transfer in COMT exhibit limited variations with respect to the size of the QM region, which ranges from ∼100 to ∼500 atoms in this study. With extensive sampling, local and global structural characteristics of the enzyme are largely conserved across the studied QM regions, while the nature of the transition state (e.g., secondary kinetic isotope effect) and reaction exergonicity are largely maintained. Deviations in the free energy profile with different QM region sizes are similar in magnitude to those observed with changes in other simulation protocols, such as different initial enzyme conformations and boundary conditions. Electronic structural properties, such as the covariance matrix of residual charge fluctuations, appear to exhibit rather long-range correlations, especially when the peptide backbone is included in the QM region; this observation holds when a range-separated DFT approach is used as the QM region, suggesting that delocalization error is unlikely the origin. Overall, the analyses suggest that multiple simulation details determine the results of QM/MM enzyme simulations with comparable contributions.

Efficiently Computing Excitations of Complex Systems: Linear-Scaling Time-Dependent Embedded Mean-Field Theory in Implicit Solvent

Quantum embedding schemes have the potential to significantly reduce the computational cost of first-principles calculations while maintaining accuracy, particularly for calculations of electronic excitations in complex systems. In this work, I combine time-dependent embedded mean field theory (TD-EMFT) with linear-scaling density functional theory and implicit solvation models, extending previous work within the ONETEP code. This provides a way to perform multilevel calculations of electronic excitations on very large systems, where long-range environmental effects, both quantum and classical in nature, are important. I demonstrate the power of this method by performing simulations on a variety of systems, including a molecular dimer, a chromophore in solution, and a doped molecular crystal. This work paves the way for high accuracy calculations to be performed on large-scale systems that were previously beyond the reach of quantum embedding schemes.

The role of hydrogen bonds and electrostatic interactions in enhancing two-photon absorption in green and yellow fluorescent proteins

The spectral properties of fluorescent proteins (FPs) depend on the protein environment of the chromophore (CRO). A deeper understanding of the CRO - environment interactions in terms of FPs spectral characteristics will allow for a rational design of novel markers with desired properties. Here, we are taking a step towards achieving this important goal. With the time-dependent density functional theory (TDDFT), we calculate one- and two-photon absorption (OPA and TPA) spectra for 5 green FPs (GFPs) and 3 yellow FPs (YFPs) differing in amino acid sequence. The goal is to reveal a role of: (i) electrostatic interactions, (ii) hydrogen-bonds (h-bonds), and (iii) h-bonds together with distant electrostatic field in absorption spectra tuning. Our results point to design hypothesis towards FPs optimised for TPA-based applications. Both h-bonds and electrostatic interactions co-operate in enhancing TPA cross-section (σ TPA ) for the S 0 ->S 1   transition in GFPs. Furthermore, it seems that details of h-bonds network in the CRO's vicinity influences σ TPA   response to CRO - environment electrostatic interactions in YFPs. We postulate that engineering FPs with more hydrophilic CRO's environment can lead to greater σ TPA . We also find that removing h-bonds formed with the CRO's phenolate leads to TPA enhancement for transition to higher excited states than S 1 . Particularly Y145 and T203 residues are important in this regard.

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.

Bridging the Experiment-Calculation Divide: Machine Learning Corrections to Redox Potential Calculations in Implicit and Explicit Solvent Models

Prediction of redox potentials is essential for catalysis and energy storage. Although density functional theory (DFT) calculations have enabled rapid redox potential predictions for numerous compounds, prominent errors persist compared to experimental measurements. In this work, we develop machine learning (ML) models to reduce the errors of redox potential calculations in both implicit and explicit solvent models. Training and testing of the ML correction models are based on the diverse ROP313 data set with experimental redox potentials measured for organic and organometallic compounds in a variety of solvents. For the implicit solvent approach, our ML models can reduce both the systematic bias and the number of outliers. ML corrected redox potentials also demonstrate less sensitivity to DFT functional choice. For the explicit solvent approach, we significantly reduce the computational costs by embedding the microsolvated cluster in implicit bulk solvent, obtaining converged redox potential results with a smaller solvation shell. This combined implicit-explicit solvent model, together with GPU-accelerated quantum chemistry methods, enabled rapid generation of a large data set of explicit-solvent-calculated redox potentials for 165 organic compounds, allowing detailed investigation of the error sources in explicit solvent redox potential calculations.

Quantum-Mechanical/Molecular-Mechanical (QM/MM) Simulations for Understanding Enzyme Dynamics

Quantum mechanics/molecular mechanics (QM/MM) methods have become widely used for computational modeling of enzyme structure and mechanism. In these approaches, a portion of the enzyme of great interest (e.g., where a chemical reaction is occurring) is treated with QM, whereas the surrounding region is treated with MM. A critical challenge with these methods is the choice of the region to partition into QM and which to treat with MM along with numerous practical choices that must be made at each step of the modeling procedure. Here, we attempt to simplify this process by describing the steps involved in preparing protein structures, choosing the appropriate QM region size and electronic structure methods, preparing all necessary input files, and troubleshooting common errors for QM/MM simulations of enzymes.

Flexible boundary layer using exchange for embedding theories. I. Theory and implementation

Embedding theory is a powerful computational chemistry approach to exploring the electronic structure and dynamics of complex systems, with Quantum Mechanical/Molecular Mechanics (QM/MM) being the prime example. A challenge arises when trying to apply embedding methodology to systems with diffusible particles, e.g., solvents, if some of them must be included in the QM region, for example, in the description of solvent-supported electronic states or reactions involving proton transfer or charge-transfer-to-solvent: without a special treatment, inter-diffusion of QM and MM particles will eventually lead to a loss of QM/MM separation. We have developed a new method called Flexible Boundary Layer using Exchange (FlexiBLE) that solves the problem by adding a biasing potential to the system that closely maintains QM/MM separation. The method rigorously preserves ensemble averages by leveraging their invariance to an exchange of identical particles. With a careful choice of the biasing potential and the use of a tree algorithm to include only important QM and MM exchanges, we find that the method has an MM-forcefield-like computational cost and thus adds negligible overhead to a QM/MM simulation. Furthermore, we show that molecular dynamics with the FlexiBLE bias conserves total energy, and remarkably, sub-diffusional dynamical quantities in the inner QM region are unaffected by the applied bias. FlexiBLE thus widens the range of chemistry that can be studied with embedding theory.

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.

Removing artifacts in polarizable embedding calculations of one- and two-photon absorption spectra of fluorescent proteins

The multiscale calculations involving excited states may suffer from the electron spill-out (ESO) problem. This seems to be especially the case when the environment of the core region, described with the electronic structure method, is approximated by a polarizable force field. The ESO effect often leads to incorrect physical character of electronic excitations, spreading outside the quantum region, which, in turn, results in erroneous absorption spectra. In this work, we investigate means to remove the artifacts in one-photon absorption (OPA) and two-photon absorption (TPA) spectra of green and yellow fluorescent protein representatives. This includes (i) using different basis sets, (ii) extending the core subsystem beyond the chromophore, (iii) modification of polarization interaction between the core region and its environment, and (iv) including the Pauli repulsion through effective core potentials (ECPs). Our results clearly show that ESO is observed when diffuse functions are used to assemble the multielectron wave function regardless of the exchange-correlation functional used. Furthermore, extending the core region, thus accounting for exchange interactions between the chromophore and its environment, leads to even more spurious excited states. Also, damping the interactions between the core subsystem and the polarizable force field is hardly helpful. In contrast, placing ECPs in the position of sites creating the embedding potential leads to the removal of artificious excited states that presumably should not be observed in the OPA and TPA spectra. We prove that it is a reliable and cost-effective approach for systems where the covalent bond(s) between the core region and its environment must be cut.

Open-Source Multi-GPU-Accelerated QM/MM Simulations with AMBER and QUICK

The quantum mechanics/molecular mechanics (QM/MM) approach is an essential and well-established tool in computational chemistry that has been widely applied in a myriad of biomolecular problems in the literature. In this publication, we report the integration of the QUantum Interaction Computational Kernel (QUICK) program as an engine to perform electronic structure calculations in QM/MM simulations with AMBER. This integration is available through either a file-based interface (FBI) or an application programming interface (API). Since QUICK is an open-source GPU-accelerated code with multi-GPU parallelization, users can take advantage of "free of charge" GPU-acceleration in their QM/MM simulations. In this work, we discuss implementation details and give usage examples. We also investigate energy conservation in typical QM/MM simulations performed at the microcanonical ensemble. Finally, benchmark results for two representative systems in bulk water, the N-methylacetamide (NMA) molecule and the photoactive yellow protein (PYP), show the performance of QM/MM simulations with QUICK and AMBER using a varying number of CPU cores and GPUs. Our results highlight the acceleration obtained from a single or multiple GPUs; we observed speedups of up to 53× between a single GPU vs a single CPU core and of up to 2.6× when comparing four GPUs to a single GPU. Results also reveal speedups of up to 3.5× when the API is used instead of FBI.

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.

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.

ReaxFF/AMBER-A Framework for Hybrid Reactive/Nonreactive Force Field Molecular Dynamics Simulations

Combined quantum mechanical/molecular mechanical (QM/MM) models using semiempirical and ab initio methods have been extensively reported on over the past few decades. These methods have been shown to be capable of providing unique insights into a range of problems, but they are still limited to relatively short time scales, especially QM/MM models using ab initio methods. An intermediate approach between a QM based model and classical mechanics could help fill this time-scale gap and facilitate the study of a range of interesting problems. Reactive force fields represent the intermediate approach explored in this paper. A widely used reactive model is ReaxFF, which has largely been applied to materials science problems and is generally used as a stand-alone (i.e., the full system is modeled using ReaxFF). We report a hybrid ReaxFF/AMBER molecular dynamics (MD) tool, which introduces ReaxFF capabilities to capture bond breaking and formation within the AMBER MD software package. This tool enables us to study local reactive events in large systems at a fraction of the computational costs of QM/MM models. We describe the implementation of ReaxFF/AMBER, validate this implementation using a benzene molecule solvated in water, and compare its performance against a range of similar approaches. To illustrate the predictive capabilities of ReaxFF/AMBER, we carried out a Claisen rearrangement study in aqueous solution. In a first for ReaxFF, we were able to use AMBER's potential of mean force (PMF) capabilities to perform a PMF study on this organic reaction. The ability to capture local reaction events in large systems using combined ReaxFF/AMBER opens up a range of problems that can be tackled using this model to address both chemical and biological processes.

The Mechanics of the Bicycle Pedal Photoisomerization in Crystalline cis,cis-1,4-Diphenyl-1,3-butadiene

Direct irradiation of crystalline cis,cis-1,4-diphenyl-1,3-butadiene (cc-DPB) forms trans,trans-1,4-diphenyl-1,3,-butadiene via a concerted two-bond isomerization called the bicycle pedal (BP) mechanism. However, little is known about photoisomerization pathways in the solid state and there has been much debate surrounding the interpretation of volume-conserving isomerization mechanisms. The bicycle pedal photoisomerization is investigated using the quantum mechanics/molecular mechanics complete active space self-consistent field/Amber force-field method. Important details about how the steric environment influences isomerization mechanisms are revealed including how the one-bond flip and hula-twist mechanisms are suppressed by the crystal cavity, the nature of the seam space in steric environments, and the features of the bicycle pedal mechanism. Specifically, in the bicycle pedal, the phenyl rings of cc-DPB are locked in place and the intermolecular packing allows a passageway for rotation of the central diene in a volume-conserving manner. In contrast, the bicycle pedal rotation in the gas phase is not a stable pathway, so single-bond rotation mechanisms become operative instead. Furthermore, the crystal BP mechanism is an activated process that occurs completely on the excited state; the photoproduct can decay to the ground state through radiative and non-radiative pathways. The present models, however, do not capture the quantitative activation barriers, and more work is needed to better model reactions in crystals. Last, the reaction barriers of the different crystalline conformations within the unit cell of cc-DPB are compared to investigate the possibility for conformation-dependent isomerization. Although some difference in reaction barriers is observed, the difference is most likely not responsible for the experimentally observed periods of fast and slow conversion.

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.

Exchange Repulsion in Quantum Mechanical/Effective Fragment Potential Excitation Energies: Beyond Polarizable Embedding

Hybrid quantum mechanical and molecular mechanical (QM/MM) approaches facilitate computational modeling of large biological and materials systems. Typically, in QM/MM, a small region of the system is modeled with an accurate quantum mechanical method and its surroundings with a more efficient alternative, such as a classical force field or the effective fragment potential (EFP). The reliability of QM/MM calculations depends largely on the treatment of interactions between the two subregions, also known as embedding. The polarizable embedding, which allows mutual polarization between solvent and solute, is considered to be essential for describing electronic excitations in polar solvents. In this work, we employ the QM/EFP model and extend the polarizable embedding by incorporating two short-range terms-a charge penetration correction to the electrostatic term and the exchange-repulsion term-both of which are modeled with one-electron contributions to the quantum Hamiltonian. We evaluate the accuracy of these terms by computing excitation energies across 37 molecular clusters consisting of biologically relevant chromophores surrounded by polar solvent molecules. QM/EFP excitation energies are compared to the fully quantum mechanical calculations with the configuration interaction singles (CIS) method. We find that the charge penetration correction diminishes the accuracy of the QM/EFP calculations. On the other hand, while the effect of exchange-repulsion is negligible for most ππ* transitions, the exchange-repulsion significantly improves description of nπ* transitions with blue solvatochromic shifts. As a result, addition of the exchange-repulsion term improves the overall accuracy of QM/EFP. Performances of QM/EFP models remain similar when excitation energies are modeled with cc-pVDZ and aug-cc-pVDZ basis sets.