Domain decomposition for implicit solvation models

Bridging Eulerian and Lagrangian Poisson-Boltzmann solvers by ESES

The Poisson-Boltzmann (PB) model is a widely used electrostatic model for biomolecular solvation analysis. Formulated as an elliptic interface problem, the PB model can be numerically solved on either Eulerian meshes using finite difference/finite element methods or Lagrangian meshes using boundary element methods. Molecular surface generators, which produce the discretized dielectric interfaces between solutes and solvents, are critical factors in determining the accuracy and efficiency of the PB solvers. In this work, we investigate the utility of the Eulerian Solvent Excluded Surface (ESES) software for rendering conjugated Eulerian and Lagrangian surface representations, which enables us to numerically validate and compare the quality of Eulerian PB solvers, such as the MIBPB solver, and the Lagrangian PB solvers, such as the TABI-PB solver. Furthermore, with the ESES software and its associated PB solvers, we are able to numerically validate an interesting and useful but often neglected source-target symmetric property associated with the linearized PB model.

Multireference Protonation Energetics of a Dimeric Model of Nitrogenase Iron-Sulfur Clusters

Characterizing the electronic structure of the iron-sulfur clusters in nitrogenase is necessary to understand their role in the nitrogen fixation process. One challenging task is to determine the protonation state of the intermediates in the nitrogen fixing cycle. Here, we use a dimeric iron-sulfur model to study relative energies of protonation at C, S, or Fe. Using a composite method based on coupled cluster and density matrix renormalization group energetics, we converge the relative energies of four protonated configurations with respect to basis set and correlation level. We find that accurate relative energies require large basis sets as well as a proper treatment of multireference and relativistic effects. We have also tested ten density functional approximations for these systems. Most of them give large errors in their relative energies. The best performing functional in this system is B3LYP, which gives mean absolute and maximum deviations of only 10 and 13 kJ/mol with respect to our correlated wave function estimates, respectively, comparable to the uncertainty in our correlated estimates. Our work provides benchmark results for the calibration of new approximate electronic structure methods and density functionals for these problems.

KineticNet: Deep learning a transferable kinetic energy functional for orbital-free density functional theory

Orbital-free density functional theory (OF-DFT) holds promise to compute ground state molecular properties at minimal cost. However, it has been held back by our inability to compute the kinetic energy as a functional of electron density alone. Here, we set out to learn the kinetic energy functional from ground truth provided by the more expensive Kohn-Sham density functional theory. Such learning is confronted with two key challenges: Giving the model sufficient expressivity and spatial context while limiting the memory footprint to afford computations on a GPU and creating a sufficiently broad distribution of training data to enable iterative density optimization even when starting from a poor initial guess. In response, we introduce KineticNet, an equivariant deep neural network architecture based on point convolutions adapted to the prediction of quantities on molecular quadrature grids. Important contributions include convolution filters with sufficient spatial resolution in the vicinity of nuclear cusp, an atom-centric sparse but expressive architecture that relays information across multiple bond lengths, and a new strategy to generate varied training data by finding ground state densities in the face of perturbations by a random external potential. KineticNet achieves, for the first time, chemical accuracy of the learned functionals across input densities and geometries of tiny molecules. For two-electron systems, we additionally demonstrate OF-DFT density optimization with chemical accuracy.

Extended Conductor-like Polarizable Continuum Solvation Model (CPCM-X) for Semiempirical Methods

We have developed a new method to accurately account for solvation effects in semiempirical quantum mechanics based on a polarizable continuum model (PCM). The extended conductor-like polarizable continuum model (CPCM-X) incorporates a computationally efficient domain decomposition conductor-like screening model (ddCOSMO) for extended tight binding (xTB) methods and uses a post-processing approach based on established solvation models, like the conductor-like screening model for real solvents (COSMO-RS) and the universal solvent model based on solute electron density (SMD). According to various benchmarks, the approach performs well across a broad range of systems and applications, including hydration free energies, non-aqueous solvation free energies, and large supramolecular association reactions of neutral and charged species. Our method for computing solvation free energies is much more accurate than the current methods in the xtb program package. It improves the accuracy of solvation free energies by up to 40% for larger supramolecular association reactions to match even the accuracy of higher-level DFT-based solvation models like COSMO-RS and SMD while being computationally more than 2 orders of magnitude faster. The proposed method and the underlying ddCOSMO model are readily available for a wide variety of solvents and are accessible in xtb for use in various computational applications.

Constrained Nuclear-Electronic Orbital Density Functional Theory with a Dielectric Continuum Solvent Model

Solvent effects are crucial for simulating chemical and biological processes in solutions. The continuum solvation model is widely used for incorporating solvent effects with different levels of theoretical descriptions of solutes. For solutes and solutions containing hydrogen atoms, nuclear quantum effects can also be nonnegligible for reliable simulations. In this work, we couple our recently developed constrained nuclear-electronic orbital density functional theory with a dielectric continuum solvation model to cover nuclear quantum effects and solvent effects simultaneously. This approach is applied to the formate ion, where an anomalous solvatochromic shift in C-H stretch frequency was reported in experiments. By using this new approach to account for nuclear quantum effects and solvent effects, we show that the vibrational frequency of the C-H stretch and the solvatochromic shift are accurately described.

Core-Projected Hybrids Fix Systematic Errors in Time-Dependent Density Functional Theory Predicted Core-Electron Excitations

Linear response time-dependent density functional theory (TDDFT) is widely applied to valence, Rydberg, and charge-transfer excitations but, in its current form, makes large errors for core-electron excitations. This work demonstrates that the admixture of nonlocal exact exchange in atomic core regions significantly improves TDDFT-predicted core excitations. Exact exchange admixture is accomplished using projected hybrid density functional theory [ J. Chem. Theory Comput. 2023, 19, 837-847]. Scalar relativistic TDDFT calculations using core-projected B3LYP accurately model core excitations of second-period elements C-F and third-period elements Si-Cl, without sacrificing performance for the relative shifts of core excitation energies. Predicted K-edge X-ray near absorption edge structure (XANES) of a series of sulfur standards highlight the value of this approach. Core-projected hybrids appear to be a practical solution to TDDFT's limitations for core excitations, in the way that long-range-corrected hybrids are a practical solution to TDDFT's limitations for Rydberg and charge-transfer excitations.

Solv: An Alternative Continuum Model Implementation Based on Fixed Atomic Charges, Scaled Particle Theory, and the Atom-Atom Potential Method

An alternative continuum model implementation is reported. The electrostatic contribution to the solvation Gibbs free energy utilizes the noniterative conductor-like screening model of Vyboishchikov and Voityuk (DOI: 10.1002/jcc.26531) based on the fixed partial atomic charges. The nonelectrostatic solute-solvent dispersion-repulsion energy is calculated through the Caillet-Claverie atom-atom potential method employing the grid-based approach. The nonelectrostatic cavitation energy is computed within the scaled particle theory (SPT) formalism with the solute hard-sphere radius obtained via the Pierotti-Claverie (PC) scheme, from the solute molecular surface (SPT-S) or volume (SPT-V). The solvent hard-sphere radius is derived through the fitting to the experimental total solvation free energies of 2530 neutral species in 92 solvents. Application of the model to reproduce both the absolute and relative (reaction net) solvation free energies indicates that the SPT-V approach based on the CM5 charges is the best performer. The method is suggested for the solvation free energy calculation in the nonaqueous solvents.

Linear scaling computation of forces for the domain-decomposition linear Poisson-Boltzmann method

The Linearized Poisson-Boltzmann (LPB) equation is a popular and widely accepted model for accounting solvent effects in computational (bio-) chemistry. In the present article, we derive the analytical forces using the domain-decomposition-based LPB-method with a van-der Waals or solvent-accessible surface. We present an efficient strategy to compute the forces and its implementation, allowing linear scaling of the method with respect to the number of atoms using the fast multipole method. Numerical tests illustrate the accuracy of the computation of the analytical forces and compare the efficiency with other available methods.

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.

Linearly scaling computation of ddPCM solvation energy and forces using the fast multipole method

This paper proposes the first linear scaling implementation for the domain decomposition approach of the polarizable continuum model (ddPCM) for the computation of the solvation energy and forces. The ddPCM-equation consists of a (non-local) integral equation on the van der Waals or solvent accessible surface of the solute's cavity resulting in a dense solution matrix, and, in turn, one matrix-vector multiplication has a quadratic arithmetic complexity with respect to the number of atoms of the solute molecule. The use of spherical harmonics as basis functions makes it natural to employ the fast multipole method (FMM) in order to provide an asymptotically linear scaling method. In this paper, we employ the FMM in a non-uniform manner with a clusterization based on a recursive inertial bisection. We present some numerical tests illustrating the accuracy and scaling of our implementation.

On the Convergence of the Damped Additive Schwarz Methods and the Subdomain Coloring

In this paper, we consider that the subdomains of the domain decomposition are colored such that the subdomains with the same color do not intersect and introduce and analyze the convergence of a damped additive Schwarz method related to such a subdomain coloring for the resolution of variational inequalities and equations. In this damped method, a single damping value is associated with all the subdomains having the same color. We first make this analysis both for variational inequalities and, as a special case, for equations in an abstract framework. By introducing an assumption on the decomposition of the convex set of the variational inequality, we theoretically analyze in a reflexive Banach space the convergence of the damped additive Schwarz method. The introduced assumption contains a constant C0, and we explicitly write the expression of the convergence rates, depending on the number of colors and the constant C0, and find the values of the damping constants which minimize them. For problems in the finite element spaces, we write the constant C0 as a function of the overlap parameter of the domain decomposition and the number of colors of the subdomains. We show that, for a fixed overlap parameter, the convergence rate, as a function of the number of subdomains has an upper limit which depends only on the number of the colors of the subdomains. Obviously, this limit is independent of the total number of subdomains. Numerical results are in agreement with the theoretical ones. They have been performed for an elasto-plastic problem to verify the theoretical predictions concerning the choice of the damping parameter, the dependence of the convergence on the overlap parameter and on the number of subdomains.

Molecular Dynamics Simulations Enforcing Nonperiodic Boundary Conditions: New Developments and Application to the Solvent Shifts of Nitroxide Magnetic Parameters

Multiscale methods combining quantum mechanics and molecular mechanics (QM/MM) have become the most suitable and effective strategies for the investigation of the spectroscopic properties of medium-to-large size chromophores in condensed phases. In this context, we are developing a novel workflow aimed at improving the generality, reliability, and ease of use of the available computational tools. In this paper, we report our latest developments with specific reference to a general protocol based on atomistic simulations, carried out under nonperiodic boundary conditions (NPBC). In particular, we add to our in house MD engine a new efficient treatment of mean field electrostatic contributions to energy and forces, together with the capability of performing the simulations either in the canonical (NVT) or in the isothermal–isobaric (NPT) ensemble. Next, we provide convincing evidence that the NBPC approach enhanced by specific tweaks for rigid body propagation, allows for the simulation of solute–solvent systems with a minimum number of degrees of freedom and large integration time step. After its validation, this new approach is applied to the challenging case of solvatochromic effects on the electron paramagnetic resonance (EPR) spectrum of a prototypical nitroxide radical. To this end, we propose and validate also an automated protocol to extract and weight simulation snapshots, making use of a continuous description of the strength of solute–solvent hydrogen bridges. While further developments are being worked on, the effectiveness of our approach, even in its present form, is demonstrated by the accuracy of the results obtained through an unsupervised approach characterized by a strongly reduced computational cost as compared to that of conventional QM/MM models, without any appreciable deterioration of accuracy.

Electronic excitations through the prism of mean-field decomposition techniques

The potential of mean-field decomposition techniques in interpreting electronic transitions in molecules is explored, in particular, the usefulness of these for offering computational signatures of different classes of such excitations. When viewed as a conceptual lens for this purpose, decomposed results are presented for ground- and excited-state energies and dipole moments of selected prototypical organic dyes, and the discrete nature of these properties as well as how they change upon transitioning from one state to another is analyzed without recourse to a discussion based on the involved molecular orbitals. On the basis of results obtained both with and without an account of continuum solvation, our work is further intended to shed new light on practical and pathological differences in between various functional approximations in orbital-optimized Kohn-Sham density functional theory for excited states, equipping practitioners and developers in the field with new probes and possible validation tools.

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.

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

We present the implementation of a fully coupled polarizable QM/MM/continuum model based on the AMOEBA polarizable force field and the domain decomposition implementation of the conductor-like screening model. Energies, response properties, and analytical gradients with respect to both QM and MM nuclear positions are available, and a generic, atomistic cavity can be employed. The model is linear scaling in memory requirements and computational cost with respect to the number of classical atoms and is therefore suited to model large, complex systems. Using three variants of the green-fluorescent protein, we investigate the overall computational cost of such calculations and the effect of the continuum model on the convergence of the computed properties with respect to the size of the embedding. We also demonstrate the fundamental role of polarization effects by comparing polarizable and nonpolarizable embeddings to fully QM ones.

Quantum chemistry for molecules at extreme pressure on graphical processing units: Implementation of extreme-pressure polarizable continuum model

Pressure plays essential roles in chemistry by altering structures and controlling chemical reactions. The extreme-pressure polarizable continuum model (XP-PCM) is an emerging method with an efficient quantum mechanical description of small- and medium-sized molecules at high pressure (on the order of GPa). However, its application to large molecular systems was previously hampered by a CPU computation bottleneck: the Pauli repulsion potential unique to XP-PCM requires the evaluation of a large number of electric field integrals, resulting in significant computational overhead compared to the gas-phase or standard-pressure polarizable continuum model calculations. Here, we exploit advances in graphical processing units (GPUs) to accelerate the XP-PCM-integral evaluations. This enables high-pressure quantum chemistry simulation of proteins that used to be computationally intractable. We benchmarked the performance using 18 small proteins in aqueous solutions. Using a single GPU, our method evaluates the XP-PCM free energy of a protein with over 500 atoms and 4000 basis functions within half an hour. The time taken by the XP-PCM-integral evaluation is typically 1% of the time taken for a gas-phase density functional theory (DFT) on the same system. The overall XP-PCM calculations require less computational effort than that for their gas-phase counterpart due to the improved convergence of self-consistent field iterations. Therefore, the description of the high-pressure effects with our GPU-accelerated XP-PCM is feasible for any molecule tractable for gas-phase DFT calculation. We have also validated the accuracy of our method on small molecules whose properties under high pressure are known from experiments or previous theoretical studies.

Computing Protein pKas Using the TABI Poisson–Boltzmann Solver

A common approach to computing protein pKas uses a continuum dielectric model in which the protein is a low dielectric medium with embedded atomic point charges, the solvent is a high dielectric medium with a Boltzmann distribution of ionic charges, and the pKa is related to the electrostatic free energy which is obtained by solving the Poisson–Boltzmann equation. Starting from the model pKa for a titrating residue, the method obtains the intrinsic pKa and then computes the protonation probability for a given pH including site–site interactions. This approach assumes that acid dissociation does not affect protein conformation aside from adding or deleting charges at titratable sites. In this work, we demonstrate our treecode-accelerated boundary integral (TABI) solver for the relevant electrostatic calculations. The pKa computing procedure is enclosed in a convenient Python wrapper which is publicly available at the corresponding author’s website. Predicted results are compared with experimental pKas for several proteins. Among ongoing efforts to improve protein pKa calculations, the advantage of TABI is that it reduces the numerical errors in the electrostatic calculations so that attention can be focused on modeling assumptions.

Implicit Solvation Using a Generalized Finite-Difference Approach in CRYSTAL: Implementation and Results for Molecules, Polymers, and Surfaces

We present the implementation of an implicit solvation model in the CRYSTAL code. The solvation energy is separated into two components: the electrostatic contribution arising from a self-consistent reaction field treatment obtained within a generalized finite-difference Poisson model, augmented by a nonelectrostatic contribution proportional to the solvent-accessible surface area of the solute. A discontinuous dielectric boundary is used, along with a solvent-excluded surface built from interlocking atom-centered spheres on which apparent surface point charges are mapped. The procedure is general and can be performed at both the Hartree-Fock and density functional theory levels, with pure or hybrid functionals, for systems periodic in 0, 1, and 2 directions, that is, for isolated molecules and extended polymers and surfaces. The Poisson equation resolution and apparent surface charge formalism is first validated on model analytical test cases. The good agreement obtained on solvation free energies is further confirmed by calculations performed on a large test set of 501 neutral molecules, for which a mean unsigned error of 1.3 kcal/mol is obtained when compared to the available experimental data. Importantly, the self-consistent reaction field procedure converges well for all molecules tested. This is further verified for all polymers and surfaces considered. In particular, for periodic systems, results obtained on an infinite glycine chain and on the wettability parameters of SiO₂ surfaces are in good agreement with previously published data. The size extensivity of the energetic terms involved in the electrostatic contribution to the solvation energy is also well verified. These encouraging results constitute a first step to take into account complex environments in the CRYSTAL code, potentially allowing for a more accurate modeling of complex processes for both periodic and nonperiodic systems.

Polarizable molecular dynamics in a polarizable continuum solvent

We present, for the first time, scalable polarizable molecular dynamics (MD) simulations within a polarizable continuum solvent with molecular shape cavities and exact solution of the mutual polarization. The key ingredients are a very efficient algorithm for solving the equations associated with the polarizable continuum, in particular, the domain decomposition Conductor-like Screening Model (ddCOSMO), which involves a rigorous coupling of the continuum with the polarizable force field achieved through a robust variational formulation and an effective strategy to solve the coupled equations. The coupling of ddCOSMO with nonvariational force fields, including AMOEBA, is also addressed. The MD simulations are feasible, for real-life systems, on standard cluster nodes; a scalable parallel implementation allows for further acceleration in the context of a newly developed module in Tinker, named Tinker-HP. NVE simulations are stable, and long-term energy conservation can be achieved. This paper is focused on the methodological developments, the analysis of the algorithm, and the stability of the simulations; a proof-of-concept application is also presented to attest to the possibilities of this newly developed technique.

Quantum, classical, and hybrid QM/MM calculations in solution: general implementation of the ddCOSMO linear scaling strategy

We present the general theory and implementation of the Conductor-like Screening Model according to the recently developed ddCOSMO paradigm. The various quantities needed to apply ddCOSMO at different levels of theory, including quantum mechanical descriptions, are discussed in detail, with a particular focus on how to compute the integrals needed to evaluate the ddCOSMO solvation energy and its derivatives. The overall computational cost of a ddCOSMO computation is then analyzed and decomposed in the various steps: the different relative weights of such contributions are then discussed for both ddCOSMO and the fastest available alternative discretization to the COSMO equations. Finally, the scaling of the cost of the various steps with respect to the size of the solute is analyzed and discussed, showing how ddCOSMO opens significantly new possibilities when cheap or hybrid molecular mechanics/quantum mechanics methods are used to describe the solute.

Wavelet formulation of the polarizable continuum model. II. Use of piecewise bilinear boundary elements

The sparsity pattern obtained in the wavelet Galerkin boundary element representation of the PCM boundary integral operators. By employing a wavelet basis on a smooth solvent-excluded molecular surface the method is able to guarantee solutions with high accuracy at a linear cost in memory and computational time.

Quantum Calculations in Solution for Large to Very Large Molecules: A New Linear Scaling QM/Continuum Approach

We present a new implementation of continuum solvation models for semiempirical Hamiltonians that allows the description of environmental effects on very large molecular systems. In this approach based on a domain decomposition strategy of the COSMO model (ddCOSMO), the solution to the COSMO equations is no longer the computational bottleneck but becomes a negligible part of the overall computation time. In this Letter, we analyze the computational impact of COSMO on the solution of the SCF equations for large to very large molecules, using semiempirical Hamiltonians, for both the new ddCOSMO implementation and the most recent, linear scaling one, based on the fast multipole method. A further analysis is on the simulation of the UV/visible spectrum of a light-harvesting pigment-protein complex. All of the results show how the new ddCOSMO algorithm paves the way to routine computations for large molecular systems in the condensed phase.

Scalable Evaluation of Polarization Energy and Associated Forces in Polarizable Molecular Dynamics: II.Towards Massively Parallel Computations using Smooth Particle Mesh Ewald

In this paper, we present a scalable and efficient implementation of point dipole-based polarizable force fields for molecular dynamics (MD) simulations with periodic boundary conditions (PBC). The Smooth Particle-Mesh Ewald technique is combined with two optimal iterative strategies, namely, a preconditioned conjugate gradient solver and a Jacobi solver in conjunction with the Direct Inversion in the Iterative Subspace for convergence acceleration, to solve the polarization equations. We show that both solvers exhibit very good parallel performances and overall very competitive timings in an energy-force computation needed to perform a MD step. Various tests on large systems are provided in the context of the polarizable AMOEBA force field as implemented in the newly developed Tinker-HP package which is the first implementation for a polarizable model making large scale experiments for massively parallel PBC point dipole models possible. We show that using a large number of cores offers a significant acceleration of the overall process involving the iterative methods within the context of spme and a noticeable improvement of the memory management giving access to very large systems (hundreds of thousands of atoms) as the algorithm naturally distributes the data on different cores. Coupled with advanced MD techniques, gains ranging from 2 to 3 orders of magnitude in time are now possible compared to non-optimized, sequential implementations giving new directions for polarizable molecular dynamics in periodic boundary conditions using massively parallel implementations.