Energy, Structures, and Response Properties with a Fully Coupled QM/AMOEBA/ddCOSMO Implementation

Analytical First Derivatives of the SCF Energy for the Conductor-Like Polarizable Continuum Model With Non-Static Radii

Within this work, we present the derivation and implementation of analytical gradients for the Gaussian-switching (SwiG) Conductor-like Polarizable Continuum Model (CPCM) with general nuclear coordinate-dependent non-static radii used for the creation of van der Waals-type cavities. This is done using the recently presented dynamic radii adjustment for continuum solvation (Draco) scheme. This allows for efficient geometry optimization and reasonable numerical Hessian calculations. The derived gradient is implemented in ORCA, and therefore is easily applicable. The derivation and implementation is validated by comparing analytical and numerical gradients and testing geometry optimizations on a diverse test set, including small organic compounds, metal-organic complexes, and highly charged species. We additionally test the continuity of the potential energy surface using an example where very strong changes in the radii occur. The computational efficiency of the derived gradient is investigated.

Energetic Landscape and Terminal Emitters of Phycobilisome Cores from Quantum Chemical Modeling

Phycobilisomes (PBs) are giant antenna supercomplexes of cyanobacteria that use phycobilin pigments to capture sunlight and transfer the collected energy to membrane-bound photosystems. In the PB core, phycobilins are bound to particular allophycocyanin (APC) proteins. Some phycobilins are thought to be terminal emitters (TEs) with red-shifted fluorescence. However, the precise identification of TEs is still under debate. In this work, we employ multiscale quantum-mechanical calculations to disentangle the excitation energy landscape of PB cores. Using the recent atomistic PB structures from Synechoccoccus PCC 7002 and Synechocystis PCC 6803, we compute the spectral properties of different APC trimers and assign the low-energy pigments. We show that the excitation energy of APC phycobilins is determined by geometric and electrostatic factors and is tuned by the specific protein-protein interactions within the core. Our findings challenge the simple picture of a few red-shifted bilins in the PB core and instead suggest that the red-shifts are established by the entire TE-containing APC trimers. Our work provides a theoretical microscopic basis for the interpretation of energy migration and time-resolved spectroscopy in phycobilisomes.

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.

Polarizable Embedding without Artificial Boundary Polarization

We present a fully self-consistent polarizable embedding (PE) model that does not suffer from unphysical boundary polarization. This is achieved through the use of the minimum-image convention (MIC) in the induced electrostatics. It is a simple yet effective approach that includes a more physically accurate description of the polarization throughout the molecular system. Using PE with MIC (PE-MIC), we shed new light on the limitations of commonly employed cutoff models, such as the droplet model, when used in PE calculations. Specifically, we investigate the effects of the unphysical polarization at the outer boundary by comparing induced dipoles and the associated electrostatic potentials, as well as some optical properties of solute-solvent and biomolecular systems. We show that the magnitude of the inaccuracies caused by the unphysical polarization depends on multiple parameters: the nature of the quantum subsystem and of the environment, the cutoff model and distance, and the calculated property.

Challenges in Protein QM/MM Simulations with Intra-Backbone Link Atoms

Hybrid quantum mechanical/molecular mechanical (QM/MM) simulations fuel discoveries in many fields of science including computational biochemistry and enzymology. Development of more convenient tools leads to an increase in the number of works in which mechanical insights into enzymes' mode of operation are obtained. Most commonly, these tools feature hydrogen-capping (link atom) approach to provide coupling between QM and MM subsystems across a covalent bond. Extensive studies were conducted to provide a solid foundation for the correctness of such an approach when a bond to a nonpolar MM atom is considered. However, not every task may be accomplished this way. Certain scenarios of using QM/MM in computational enzymology encourage or even necessitate the incorporation of backbone atoms into the QM region. Two out of three backbone atoms are polar, and in QM/MM with electrostatic embedding, a neighboring link atom will be hyperpolarized. Several schemes to mitigate this effect were previously proposed alongside a rigorous assessment of quantitative effects on model systems. However, it was not clear whether they may translate into qualitatively different results and how link atom hyperpolarization may manifest itself in a real-life enzymological scenario. Here, we show that the consequences of such an artifact may be severe and may completely overturn the conclusions drawn from the simulations. Our case advocates for the use of charge redistribution schemes whenever intra-backbone QM/MM boundaries are considered. Moreover, we addressed how different boundary types and charge redistribution schemes influence backbone dynamics. We showed that the results are heavily dependent on which boundary MM terms are retained, with charge alteration being of secondary importance. In the worst case, only three intra-backbone boundaries may be used with relative confidence in the adequacy of resulting simulations, irrespective of the hyperpolarization mitigation scheme. Thus, advances in the field are certainly needed to fuel new discoveries. As of now, we believe that issues raised in this work might encourage authors in the field to report what boundaries, boundary MM terms, and charge redistribution schemes they are using, so their results may be correctly interpreted.

State averaged CASSCF in AMOEBA polarizable water model for simulating nonadiabatic molecular dynamics with nonequilibrium solvation effects

This paper presents a state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular application (AMOEBA) polarizable water model, which enables rigorous simulation of non-adiabatic molecular dynamics with nonequilibrium solvation effects. The molecular orbital and configuration interaction coefficients of the solute wavefunction, and the induced dipoles on solvent atoms, are solved by minimizing the state averaged energy variationally. In particular, by formulating AMOEBA water models and the polarizable continuum model (PCM) in a unified way, the algorithms developed for computing SA-CASSCF/PCM energies, analytical gradients, and non-adiabatic couplings in our previous work can be generalized to SA-CASSCF/AMOEBA by properly substituting a specific list of variables. Implementation of this method will be discussed with the emphasis on how the calculations of different terms are partitioned between the quantum chemistry and molecular mechanics codes. We will present and discuss results that demonstrate the accuracy and performance of the implementation. Next, we will discuss results that compare three solvent models that work with SA-CASSCF, i.e., PCM, fixed-charge force fields, and the newly implemented AMOEBA. Finally, the new SA-CASSCF/AMOEBA method has been interfaced with the ab initio multiple spawning method to carry out non-adiabatic molecular dynamics simulations. This method is demonstrated by simulating the photodynamics of the model retinal protonated Schiff base molecule in water.

Coarse-Graining ddCOSMO through an Interface between Tinker and the ddX Library

The domain decomposition conductor-like screening model is an efficient way to compute the solvation energy of solutes within a polarizable continuum medium in a linear scaling computational time. Despite its efficiency, the application to very large systems is still challenging. A possibility to further accelerate the algorithm is resorting to coarse-graining strategies. In this paper we present a preliminary interface between the molecular dynamics package Tinker and the ddX library. The interface was used to test a united atom coarse-graining strategy that allowed us to push ddCOSMO to its limits by computing solvation energies on systems with up to 7 million atoms. We first present benchmarks to find an optimal discretization, and then, we discuss the performance and results obtained with fine- and coarse-grained solvation energy calculations.

The atomistic modeling of light-harvesting complexes from the physical models to the computational protocol

The function of light-harvesting complexes is determined by a complex network of dynamic interactions among all the different components: the aggregate of pigments, the protein, and the surrounding environment. Complete and reliable predictions on these types of composite systems can be only achieved with an atomistic description. In the last few decades, there have been important advances in the atomistic modeling of light-harvesting complexes. These advances have involved both the completeness of the physical models and the accuracy and effectiveness of the computational protocols. In this Perspective, we present an overview of the main theoretical and computational breakthroughs attained so far in the field, with particular focus on the important role played by the protein and its dynamics. We then discuss the open problems in their accurate modeling that still need to be addressed. To illustrate an effective computational workflow for the modeling of light harvesting complexes, we take as an example the plant antenna complex CP29 and its H111N mutant.

Piecewise Multipole-Expansion Implicit Solvation for Arbitrarily Shaped Molecular Solutes

The multipole-expansion (MPE) model is an implicit solvation model used to efficiently incorporate solvent effects in quantum chemistry. Even within the recent direct approach, the multipole basis used in MPE to express the dielectric response still solves the electrostatic problem inefficiently or not at all for solutes larger than approximately ten non-hydrogen atoms. In existing MPE parametrizations, the resulting systematic underestimation of the electrostatic solute-solvent interaction is presently compensated for by a systematic overestimation of nonelectrostatic attractive interactions. Even though the MPE model can thus reproduce experimental free energies of solvation of small molecules remarkably well, the inherent error cancellation makes it hard to assign physical meaning to the individual free-energy terms in the model, raising concerns about transferability. Here we resolve this issue by solving the electrostatic problem piecewise in 3D regions centered around all non-hydrogen nuclei of the solute, ensuring reliable convergence of the multipole series. The resulting method thus allows for a much improved reproduction of the dielectric response of a medium to a solute. Employing a reduced nonelectrostatic model with a single free parameter, in addition to the density isovalue defining the solvation cavity, our method yields free energies of solvation of neutral, anionic, and cationic solutes in water in good agreement with experiment.