Analytic energy gradient of excited electronic state within TDDFT/MMpol framework: Benchmark tests and parallel implementation

The OpenMMPol library for polarizable QM/MM calculations of properties and dynamics

We present a new library designed to provide a simple and straightforward way to implement QM/AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Applications) and other polarizable QM/MM (Molecular Mechanics) methods based on induced point dipoles. The library, herein referred to as OpenMMPol, is free and open-sourced and is engineered to address the increasing demand for accurate and efficient QM/MM simulations. OpenMMPol is specifically designed to allow polarizable QM/MM calculations of ground state energies and gradients and excitation properties. Key features of OpenMMPol include a modular architecture facilitating extensibility, parallel computing capabilities for enhanced performance on modern cluster architectures, a user-friendly interface for intuitive implementation, and a simple and flexible structure for providing input data. To show the capabilities offered by the library, we present an interface with PySCF to perform QM/AMOEBA molecular dynamics, geometry optimization, and excited-state calculation based on (time-dependent) density functional theory.

Multistate, Polarizable QM/MM Embedding Scheme Based on the Direct Reaction Field Method: Solvatochromic Shifts, Analytical Gradients and Optimizations of Conical Intersections in Solution

We recently introduced a polarizable embedding scheme based on an integral-exact reformulation of the direct reaction field method (IEDRF) that accounts for the differential solvation of ground and excited states in QM/MM simulations. The polarization and dispersion interactions between the quantum-mechanical (QM) and molecular-mechanical (MM) regions are described by the DRF Hamiltonian, while the Pauli repulsion between explicitly treated QM electrons and the implicit electron density around MM atoms is modeled with effective core potentials. A single Hamiltonian is used for all electronic states so that Born-Oppenheimer states belonging to the same geometry are orthogonal and state crossings are well-defined. In this work, we describe the implementation of the method using graphical processing unit acceleration in TeraChem, where it is combined with multiple electronic structure methods, including Hartree-Fock, time-dependent density functional theory, and complete active space self-consistent field. In contrast with older implementations of the DRF method, integrals of the polarization operators are evaluated exactly. Expressions for ingredients needed to construct analytical gradients and nonadiabatic coupling vectors are derived and tested by optimizing a conical intersection between two excited states in the presence of a polarizable solvent shell. The method is applied to estimate the solvent shifts of absorption energies of a series of donor-acceptor dyes having low-lying charge-transfer states. Even for a nonpolar solvent such as n-hexane, the inclusion of its static polarizability leads to non-negligible shifts that improve the agreement to essentially quantitative levels (0.03 eV) with full-system calculations. Good agreement with the positions of the experimental absorption maxima measured in solution is also observed.

Transferable Potential Function for Flexible H2O Molecules Based on the Single-Center Multipole Expansion

A potential function is presented for describing a system of flexible H₂O molecules based on the single-center multipole expansion (SCME) of the electrostatic interaction. The model, referred to as SCME/f, includes the variation of the molecular quadrupole moment as well as the dipole moment with changes in bond length and angle so as to reproduce results of high-level electronic structure calculations. The multipole expansion also includes fixed octupole and hexadecapole moments, as well as anisotropic dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole polarizability tensors. The model contains five adjustable parameters related to the repulsive interaction and damping functions in the electrostatic and dispersion interactions. Their values are adjusted to reproduce the lowest energy isomers of small clusters, (H₂O)_n with n = 2-6, as well as measured properties of the ice Ih crystal. Subsequent calculations of the energy difference between the various isomer configurations of the clusters show that SCME/f gives good agreement with results of electronic structure calculations and represents a significant improvement over the previously presented rigid SCME potential function. Analysis of the vibrational frequencies of the clusters and structural properties of ice Ih crystal show the importance of accurately describing the variation of the quadrupole moment with molecular structures.

Conical Intersections in Solution with Polarizable Embedding: Integral-Exact Direct Reaction Field

A common strategy to exploring the properties and reactivity of complex systems is to use quantum mechanics/molecular mechanics (QM/MM) embedding, wherein a QM region is defined and treated with electronic structure theory, and the remainder of the system is treated with a force field. Important to the description of electronic excited states, especially those of charge-transfer character, is the treatment of the coupling between the QM and MM subsystems. The state of the art is to use a polarizable force field for the MM region and mutually couple the QM wavefunction and MM induced dipoles, in addition to the usual electrostatic embedding, yielding a polarizable embedding (QM/MM-Pol) approach. However, we showed previously that current popular QM/MM-Pol approaches exhibit issues of root flipping and/or incorrect descriptions of electronic crossings in multistate calculations [J. Chem. Theory Comput. 14, 2137 (2018)]. Here, we demonstrate a solution to these problems with an integral-exact reformulation of the direct reaction field approach of Thole and Van Duijnen (QM/MM-IEDRF). The resulting embedding potential includes one- and two-electron operators and many-body dipole-induced dipole interactions and thus includes a natural description of the screening of electron-electron interactions by the MM induced dipoles. Pauli repulsion from the environment is mimicked by effective core potentials on the MM atoms. Inherent to the DRF approach is the assumption that MM dipoles respond instantaneously to the positions of the QM electrons; therefore, dispersion interactions are captured approximately. All electronic states are eigenfunctions of the same Hamiltonian, while the polarization induced in the environment and the associated energetic stabilization are unique to each state. This allows for a consistent definition of transition properties and state crossings. We demonstrate QM/MM-IEDRF by exploring the influence of a (polarizable) inert xenon matrix environment on the conical intersection underlying the photoisomerization of ethylene.

Evaluation of molecular photophysical and photochemical properties using linear response time-dependent density functional theory with classical embedding: Successes and challenges

Time-dependent density functional theory (TDDFT) based approaches have been developed in recent years to model the excited-state properties and transition processes of the molecules in the gas-phase and in a condensed medium, such as in a solution and protein microenvironment or near semiconductor and metal surfaces. In the latter case, usually, classical embedding models have been adopted to account for the molecular environmental effects, leading to the multi-scale approaches of TDDFT/polarizable continuum model (PCM) and TDDFT/molecular mechanics (MM), where a molecular system of interest is designated as the quantum mechanical region and treated with TDDFT, while the environment is usually described using either a PCM or (non-polarizable or polarizable) MM force fields. In this Perspective, we briefly review these TDDFT-related multi-scale models with a specific emphasis on the implementation of analytical energy derivatives, such as the energy gradient and Hessian, the nonadiabatic coupling, the spin–orbit coupling, and the transition dipole moment as well as their nuclear derivatives for various radiative and radiativeless transition processes among electronic states. Three variations of the TDDFT method, the Tamm–Dancoff approximation to TDDFT, spin–flip DFT, and spin-adiabatic TDDFT, are discussed. Moreover, using a model system (pyridine–Ag20 complex), we emphasize that caution is needed to properly account for system–environment interactions within the TDDFT/MM models. Specifically, one should appropriately damp the electrostatic embedding potential from MM atoms and carefully tune the van der Waals interaction potential between the system and the environment. We also highlight the lack of proper treatment of charge transfer between the quantum mechanics and MM regions as well as the need for accelerated TDDFT modelings and interpretability, which calls for new method developments.

Polarizable Embedding with a Transferable H2O Potential Function I: Formulation and Tests on Dimer

The incorporation of mutual polarization in multiscale simulations where different regions of the system are treated at different level of theory is important in studies of, for example, electronic excitations and charge transfer processes. We present here an energy functional for describing a quantum mechanics/molecular mechanics (QM/MM) scheme that includes reciprocal polarization between the two subsystems. The inclusion of polarization alleviates shortcomings inherent in electrostatic embedding QM/MM models based on point-charge force fields. A density functional theory (DFT) description of the QM subsystem is coupled to a single center multipole expansion (SCME) description of H₂O molecules in the MM subsystem that includes anisotropic dipole and quadrupole polarizability as well as static multipoles up to and including the hexadecapole. The energy functional and the coupling scheme is general and can be extended to arbitrary order in terms of both the static and induced moments. Tests of the energy surface for the H₂O dimer show that the QM/MM results lie in between the pure DFT and pure SCME values. The consistency of the many-body contributions to the energy and analytical forces is demonstrated for an H₂O pentamer.

Electronic transitions for a fully polarizable QM/MM approach based on fluctuating charges and fluctuating dipoles: Linear and corrected linear response regimes

The fully polarizable Quantum Mechanics/Molecular Mechanics (QM/MM) approach based on fluctuating charges and fluctuating dipoles, named QM/FQFμ [T. Giovannini et al., J. Chem. Theory Comput. 15, 2233 (2019)], is extended to the calculation of vertical excitation energies of solvated molecular systems. Excitation energies are defined within two different solvation regimes, i.e., linear response (LR), where the response of the MM portion is adjusted to the QM transition density, and corrected-Linear Response (cLR) in which the MM response is adjusted to the relaxed QM density, thus being able to account for charge equilibration in the excited state. The model, which is specified in terms of three physical parameters (electronegativity, chemical hardness, and polarizability) is applied to vacuo-to-water solvatochromic shifts of aqueous solutions of para-nitroaniline, pyridine, and pyrimidine. The results show a good agreement with their experimental counterparts, thus highlighting the potentialities of this approach.

Coupled Cluster Theory with Induced Dipole Polarizable Embedding for Ground and Excited States

In this work, we present the theory and implementation of the coupled cluster single and double excitations (CCSD) method combined with a classical polarizable molecular mechanics force field (MMPol) based on the induced dipole model. The method is developed to compute electronic excitation energies within the state specific (SS) and linear response (LR) formalisms for the interaction of the quantum mechanical and classical regions. Furthermore, we consider an approximate expression of the correlation energy, originally developed for CCSD with implicit solvation models, where the interaction term is linear in the coupled cluster density. This approximation allows us to include the explicit contribution of the environment to the CC equations without increasing the computational effort. The test calculations on microsolvated systems, where the CCSD/MMPol method is compared to full CCSD calculations, demonstrates the reliability of this computational protocol for all interaction schemes (errors < 2%). We also show that it is important to include induced dipoles on all atom centers of the classical region and that too diffuse functions in the basis set may be problematic due to too strong interaction with the environment.

Polarizable Force Fields for Biomolecular Simulations: Recent Advances and Applications

Realistic modeling of biomolecular systems requires an accurate treatment of electrostatics, including electronic polarization. Due to recent advances in physical models, simulation algorithms, and computing hardware, biomolecular simulations with advanced force fields at biologically relevant timescales are becoming increasingly promising. These advancements have not only led to new biophysical insights but also afforded opportunities to advance our understanding of fundamental intermolecular forces. This article describes the recent advances and applications, as well as future directions, of polarizable force fields in biomolecular simulations.

Analytic Excited State Gradients for the QM/MM Polarizable Embedded Second-Order Algebraic Diagrammatic Construction for the Polarization Propagator PE-ADC(2)

An implementation of a QM/MM embedding in a polarizable environment is presented for second-order Møller-Plesset perturbation theory, MP2, for ground state energies and molecular gradients and for the second-order Algebraic Diagrammatic Construction, ADC(2), for excitation energies and excited state molecular gradients. In this implementation of PE-MP2 and PE-ADC(2), the polarizable embedded Hartree-Fock wave function is used as uncorrelated reference state. The polarization-correlation cross terms for the ground and excited states are included in this model via an approximate coupling density. A Lagrangian formulation is used to derive the relaxed electron densities and molecular gradients. The resolution-of-the-identity approximation speeds up the calculation of four-index electron repulsion integrals in the molecular orbital basis. As a first application, the method is used to study the photophysical properties of host-guest complexes where the accuracy and weaknesses of the model are also critically examined. It is demonstrated that the ground state geometries of the full quantum mechanical calculation for the supermolecule can be well reproduced. For excited state geometries, the deviations from the supermolecular calculation are slightly larger, but still the environment effects are captured qualitatively correctly, and energy gaps between the ground and excited states are obtained with sufficient accuracy.

Polarizable Embedding for Excited-State Reactions: Dynamically Weighted Polarizable QM/MM

Recent interest in polarizable embedding methods for electronic excited states has so far been focused on optical absorption and emission spectra calculations. To explore the suitability of these methods for excited-state reactions, we constructed a simple molecular system with an electronic crossing coupled to a polarizable species: the triatomic LiFBe. We found that current polarizable QM/MM methods inadequately describe the potential energy surfaces in this system, particularly close to the electronic crossing, so we developed a new polarizable embedding method called dynamically weighted polarizable QM/MM. The new method reproduces the potential energy surfaces of LiFBe from full-system multireference configuration interaction singles and doubles calculations with near-quantitative accuracy.

Characterizing the Structures, Spectra, and Energy Landscapes Involved in the Excited-State Proton Transfer Process of Red Fluorescent Protein LSSmKate1

By applying molecular dynamics (MD) simulations and quantum chemical calculations, we have characterized the states and processes involved in the excited-state proton transfer (ESPT) of LSSmKate1. MD simulations identify two stable structures in the electronic ground state of LSSmKate1, one with a protonated chromophore and the other with a deprotonated chromophore, thus leading to two separate low-energy absorption maxima with a large energy spacing, as observed in the calculated and experimentally measured absorption spectra. Proton transfer is induced by electronic excitation. When LSSmKate1 is excited, the electrons in the chromophore are transferred from the phenol ring to the N-acylimine moiety; the acidity of a phenolic hydroxyl group is thus enhanced. The calculated potential energy curves (PECs) exhibit energetic feasibility in the generation of the fluorescent species in LSSmKate1, and the exact agreement between the calculated and experimentally measured values of the large Stokes shift further provides solid theoretical evidence for the ESPT process taking place in photoexcited LSSmKate1. The molecular environments play a significant role in the geometries and absorption/emission energies of the chromophores. Overall, TD-ωB97X-D/molecular mechanics (MM) provides a better description of the optical properties of LSSmKate1 than TD-B3LYP/MM, although it always overestimates the excitation energies.