Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies

Alchemically derived free energies are artifacted when the perturbed moiety has a nonzero net charge. The source of the artifacts lies in the effective treatment of the electrostatic interactions within and between the perturbed atoms and remaining (partial) charges in the simulated system. To treat the electrostatic interactions effectively, lattice-summation (LS) methods or cutoff schemes in combination with a reaction-field contribution are usually employed. Both methods render the charging component of the calculated free energies sensitive to essential parameters of the system like the cutoff radius or the box side lengths. Here, we discuss the results of three previously published studies of ligand binding. These studies presented estimates of binding free energies that were artifacted due to the charged nature of the ligands. We show that the size of the artifacts can be efficiently calculated and raw simulation data can be corrected. We compare the corrected results with experimental estimates and nonartifacted estimates from path-sampling methods. Although the employed correction scheme involves computationally demanding continuum-electrostatics calculations, we show that the correction estimate can be deduced from a small sample of configurations rather than from the entire ensemble. This observation makes the calculations of correction terms feasible for complex biological systems. To show the general applicability of the proposed procedure, we also present results where the correction scheme was used to correct independent free energies obtained from simulations employing a cutoff scheme or LS electrostatics. In this work, we give practical guidelines on how to apply the appropriate corrections easily.     

The effect of different cutoff schemes in molecular simulations of proteins

Molecular simulations of nanoscale systems invariably involve assumptions and approximations to describe the electrostatic interactions, which are long-ranged in nature. One approach is the use of cutoff schemes with a reaction-field contribution to account for the medium outside the cutoff scheme. Recent reports show that macroscopic properties may depend on the exact choice of cutoff schemes in modern day simulations. In this work, a systematic analysis of the effects of different cutoff schemes was performed using a set of 52 proteins. We find no statistically significant differences between using a twin-range or a single-range cutoff scheme. Applying the cutoff based on charge groups or based on atomic positions, does lead to significant differences, which is traced to the cutoff noise for energies and forces. While group-based cutoff schemes show increased cutoff noise in the potential energy, applying an atomistic cutoff leads to artificial structure in the solvent at the cutoff distance. Carefully setting the temperature control, or using an atomistic cutoff for the solute and a group-based cutoff for the solvent significantly reduces the effects of the cutoff noise, without introducing structure in the solvent. This study aims to deepen the understanding of the implications different cutoffs have on molecular dynamics simulations.     

Computation of methodology-independent single-ion solvation properties from molecular simulations. III. Correction terms for the solvation free energies, enthalpies, entropies, heat capacities, volumes, compressibilities, and expansivities of solvated ions

The raw single-ion solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system or box size) and treatment of electrostatic interactions (Coulombic, lattice-sum, or cutoff-based) used during these simulations. However, as shown by Kastenholz and Hu?nenberger [J. Chem. Phys. 124, 224501 (2006)], correction terms can be derived for the effects of: (A) an incorrect solvent polarization around the ion and an incomplete or/and inexact interaction of the ion with the polarized solvent due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite-size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site, and the possible presence of a polarized air-liquid interface or of a constraint of vanishing average electrostatic potential in the simulated system; and (D) an inaccurate dielectric permittivity of the employed solvent model. Comparison with standard experimental data also requires the inclusion of appropriate cavity-formation and standard-state correction terms. In the present study, this correction scheme is extended by: (i) providing simple approximate analytical expressions (empirically-fitted) for the correction terms that were evaluated numerically in the above scheme (continuum-electrostatics calculations); (ii) providing correction terms for derivative thermodynamic single-ion solvation properties (and corresponding partial molar variables in solution), namely, the enthalpy, entropy, isobaric heat capacity, volume, isothermal compressibility, and isobaric expansivity (including appropriate standard-state correction terms). The ability of the correction scheme to produce methodology-independent single-ion solvation free energies based on atomistic simulations is tested in the case of Na(+) hydration, and the nature and magnitude of the correction terms for derivative thermodynamic properties is assessed numerically.     

Computation of methodology-independent ionic solvation free energies from molecular simulations. II. The hydration free energy of the sodium cation

The raw ionic solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system shape, and size) and treatment of electrostatic interactions (Coulombic, lattice sum, or cutoff based) used during these simulations. In the present article, it is shown that correction terms can be derived for the effect of (A) an incorrect solvent polarization around the ion due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site and the possible presence of a liquid-vacuum interface in the simulated system. Taking the hydration free energy of the sodium cation as a test case, it is shown that the raw solvation free energies obtained using seven different types of boundary conditions and electrostatic schemes commonly used in explicit-solvent simulations (for a total of 72 simulations differing in the corresponding simulation parameters) can be corrected so as to obtain a consistent value for this quantity.     

Computation of methodology-independent ionic solvation free energies from molecular simulations. I. The electrostatic potential in molecular liquids

The computation of ionic solvation free energies from atomistic simulations is a surprisingly difficult problem that has found no satisfactory solution for more than 15 years. The reason is that the charging free energies evaluated from such simulations are affected by very large errors. One of these is related to the choice of a specific convention for summing up the contributions of solvent charges to the electrostatic potential in the ionic cavity, namely, on the basis of point charges within entire solvent molecules (M scheme) or on the basis of individual point charges (P scheme). The use of an inappropriate convention may lead to a charge-independent offset in the calculated potential, which depends on the details of the summation scheme, on the quadrupole-moment trace of the solvent molecule, and on the approximate form used to represent electrostatic interactions in the system. However, whether the M or P scheme (if any) represents the appropriate convention is still a matter of on-going debate. The goal of the present article is to settle this long-standing controversy by carefully analyzing (both analytically and numerically) the properties of the electrostatic potential in molecular liquids (and inside cavities within them). Restricting the discussion to real liquids of "spherical" solvent molecules (represented by a classical solvent model with a single van der Waals interaction site), it is concluded that (i) for Coulombic (or straight-cutoff truncated) electrostatic interactions, the M scheme is the appropriate way of calculating the electrostatic potential; (ii) for non-Coulombic interactions deriving from a continuously differentiable function, both M and P schemes generally deliver an incorrect result (for which an analytical correction must be applied); and (iii) finite-temperature effects, including intermolecular orientation correlations and a preferential orientational structure in the neighborhood of a liquid-vacuum interface, must be taken into account. Applications of these results to the computation methodology-independent ionic solvation free energies from molecular simulations will be the scope of a forthcoming article.     

Nuclear magnetic shielding constants of liquid water: insights from hybrid quantum mechanics/molecular mechanics models

We present a gauge-origin independent method for the calculation of nuclear magnetic shielding tensors of molecules in a structured and polarizable environment. The method is based on a combination of density functional theory (DFT) or Hartree-Fock wave functions with molecular mechanics. The method is unique in the sense that it includes three important properties that need to be fulfilled in accurate calculations of nuclear magnetic shielding constants: (i) the model includes electron correlation effects, (ii) the model uses gauge-including atomic orbitals to give gauge-origin independent results, and (iii) the effect of the environment is treated self-consistently using a discrete reaction-field methodology. The authors present sample calculations of the isotropic nuclear magnetic shielding constants of liquid water based on a large number of solute-solvent configurations derived from molecular dynamics simulations employing potentials which treat solvent polarization either explicitly or implicitly. For both the (17)O and (1)H isotropic shielding constants the best predicted results compare fairly well with the experimental data, i.e., they reproduce the experimental solvent shifts to within 4 ppm for the (17)O shielding and 1 ppm for the (1)H shielding.     

Effects of long-range electrostatic forces on simulated protein folding kinetics

Molecular dynamics simulations are a useful tool for characterizing protein folding pathways. There are several methods of treating electrostatic forces in these simulations with varying degrees of physical fidelity and computational efficiency. In this article, we compare the reaction field (RF) algorithm, particle-mesh Ewald (PME), and tapered cutoffs with increasing cutoff radii to address the impact of the electrostatics method employed on the folding kinetics. We quantitatively compare different methods by a correlation of quantitative measures of protein folding kinetics. The results of these comparisons show that for protein folding kinetics, the RF algorithm can quantitatively reproduce the kinetics of the more costly PME algorithm. These results not only assist the selection of appropriate algorithms for future simulations, but also give insight on the role that long-range electrostatic forces have in protein folding.     

Performance evaluation of the zero‐multipole summation method in modern molecular dynamics software

The zero‐multiple summation method (ZMM) is a cutoff‐based method for calculating electrostatic interactions in molecular dynamics simulations, utilizing an electrostatic neutralization principle as a physical basis. Since the accuracies of the ZMM have been revealed to be sufficient in previous studies, it is highly desirable to clarify its practical performance. In this paper, the performance of the ZMM is compared with that of the smooth particle mesh Ewald method (SPME), where the both methods are implemented in molecular dynamics software package GROMACS. Extensive performance comparisons against a highly optimized, parameter‐tuned SPME implementation are performed for various‐sized water systems and two protein–water systems. We analyze in detail the dependence of the performance on the potential parameters and the number of CPU cores. Even though the ZMM uses a larger cutoff distance than the SPME does, the performance of the ZMM is comparable to or better than that of the SPME. This is because the ZMM does not require a time‐consuming electrostatic convolution and because the ZMM gains short neighbor‐list distances due to the smooth damping feature of the pairwise potential function near the cutoff length. We found, in particular, that the ZMM with quadrupole or octupole cancellation and no damping factor is an excellent candidate for the fast calculation of electrostatic interactions. © 2018 Wiley Periodicals, Inc.     

pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.     

A density functional theory formulation of the reaction field model of solvent effects

It is possible to reformulate the reaction field (RF ) model of continuum solvent effects, by considering an approximate expression describing the energy changes from one ground state to another, in the frame of density functional theory (DFT ). The energy functional for an arbitrary electronic system coupled to a spin-independent electrostatic external perturbation is used to derive the well-known Born expression giving the electrostatic component of the solvation energy of an atomic ion. The approximate RF –DFT model is illustrated for a series of representative singly positive and negatively charged atomic ions. A Kohn–Sham (KS )-like formalism is then proposed to compute solvation energies within a self-consistent field scheme. The extension of the RF -DFT model to molecular systems is also outlined. © 1995 John Wiley & Sons, Inc.     

Cutoff radius effect of isotropic periodic sum method for transport coefficients of Lennard-Jones liquid

Molecular dynamics simulations of a Lennard-Jones (LJ) liquid were applied to compare the isotropic periodic sum (IPS) method [X. Wu and B. R. Brooks, J. Chem. Phys. 122, 044107 (2005)], which can reduce the calculation cost of long-range interactions, such as the Lennard-Jones and Coulombic ones, with the cutoff method for the transport coefficients which includes the self-diffusion coefficient, bulk viscosity, and thermal conductivity. The self-diffusion coefficient, bulk viscosity, and thermal conductivity were estimated with reasonable accuracy if the cutoff distance of the LJ potential for the IPS method was greater than 3sigma. The IPS method is an effective technique for estimating the transport coefficients of the Lennard-Jones liquid in a homogeneous system.     

Efficient lookup table using a linear function of inverse distance squared

The major bottleneck in molecular dynamics (MD) simulations of biomolecules exist in the calculation of pairwise nonbonded interactions like Lennard‐Jones and long‐range electrostatic interactions. Particle‐mesh Ewald (PME) method is able to evaluate long‐range electrostatic interactions accurately and quickly during MD simulation. However, the evaluation of energy and gradient includes time‐consuming inverse square roots and complementary error functions. To avoid such time‐consuming operations while keeping accuracy, we propose a new lookup table for short‐range interaction in PME by defining energy and gradient as a linear function of inverse distance squared. In our lookup table approach, densities of table points are inversely proportional to squared pair distances, enabling accurate evaluation of energy and gradient at small pair distances. Regardless of the inverse operation here, the new lookup table scheme allows fast pairwise nonbonded calculations owing to efficient usage of cache memory. © 2013 Wiley Periodicals, Inc.     

Effective solvent mediated potentials of Na+ and Cl- ions in aqueous solution: temperature dependence

The effective solvent-mediated potentials for Na(+) and Cl(-) ions in aqueous solution were calculated in a wide range of temperatures from 0 to 100 °C. The potentials have been determined using the inverse Monte Carlo approach, from the ion-ion radial distribution functions computed in 50 ns molecular dynamics simulations of ions and explicit water molecules. We further separated the effective potentials into a short-range part and an electrostatic long-range part represented by a coulombic potential with some dielectric permittivity. We adjusted the value of the dielectric permittivity to provide the fastest possible decay of the short-range potentials at larger distances. The obtained temperature dependence of the dielectric permittivity follows well the experimental data. We show also that the largest part of the temperature dependence of the effective potentials can be attributed to the temperature-dependent dielectric permittivity.     

A long-range electrostatic potential based on the Wolf method charge-neutral condition

Molecular simulations rely heavily on a long range electrostatic Coulomb interaction. The Coulomb potential decays inversely with distance, indicating infinite effective range. In practice, molecular simulations do not directly take into account such an infinite interaction. Therefore, the Ewald, fast multipole, and cutoff methods are frequently used. Although cutoff methods are implemented easily and the calculations are fast, it has been pointed out that they produce serious artifacts. Wolf and coworkers recently discovered one source of the artifacts. They found that when the total charge in a cutoff sphere disappeared, the cutoff error is dramatically suppressed. The Wolf method uses the charge-neutral principle combined with a potential damping that is realized using a complementary error function. To date, many molecular simulation studies have demonstrated the accuracy and reliability of the Wolf method. We propose a novel long-range potential that is constructed only from the charge-neutral condition of the Wolf method without potential damping. We also show that three simulation systems, in which involve liquid sodium-chloride, TIP3P water, and a charged protein in explicit waters with neutralized ions using the new potential, provide accurate statistical and dielectric properties when compared with the particle mesh Ewald method.     

Towards a “Chemical” Hamiltonian

An analysis of the LCAO Hamiltonian is performed in terms of a “mixed” formulation of the second quantization for nonorthogonal orbitals, compressing the different interactions to one- and two-center terms as far as possible by performing appropriate projections. For this purpose an operator of atomic charge is also introduced, the expectation values of which are the Mulliken gross atomic populations on the individual atoms. The LCAO Hamiltonian is decomposed into terms having different physical meaning and significance: (i) sum of effective atomic Hamiltonians; (ii) the electrostatic interactions in the point-charge approximation; (iii) the electrostatic effects connected with the deviation of the actual charge distribution from the pointlike one; (iv) two-center overlap effects; (v) finite basis (“counterpoise”) correction terms related to the individual atoms; and (vi) similar finite basis correction terms with respect to the two-center interactions. Only terms of types (i) to (iv), containing no three- or four-center integrals, are considered as having physical significance. Based on the analysis of the Hamiltonian, an energy partitioning scheme is developed, and explicit expressions are given for one- and two-center (and basis extension) components of the SCF energy. The approach is also applied to the problem of intermolecular interactions, and an explicit formula is given permitting calculation of the “counterpoise” part of the supermolecule energy by properly taking into account that it depends not only on the extension of the basis, but also on the occupation of the additional orbitals in the intervening molecule—a factor completely overlooked in the usual scheme of calculations.     

Enhanced Metallophilicity in Metal–Carbene Systems: Stronger Character of Aurophilic Interactions in Solution

Metallophilicity is an essential concept that builds upon the attraction between closed shell metal ions. We report on the [M₂(bisNHC)₂]²⁺ (M=Au^I, Ag^I; NHC=N-heterocyclic carbene) systems, which display almost identical features in the solid state. However, in solution the Au₂ cation exhibits a significantly higher degree of rigidity owed to the stronger character of the aurophilic interactions. Both Au₂ and Ag₂ cationic constructs are able to accommodate Ag⁺ ions via M–M interactions, despite their inherent Coulombic repulsion. When electrostatic repulsion between host and guest is partially diminished, M–M distances are substantially shortened. Quantum chemical calculations estimate intermetallic bond orders up to 0.2. Although at the limit of (or beyond) the van der Waals radii, metallophilic interactions are responsible for their behavior in solution.     

Toward the correction of effective electrostatic forces in explicit-solvent molecular dynamics simulations: restraints on solvent-generated electrostatic potential and solvent polarization

Despite considerable advances in computing power, atomistic simulations under nonperiodic boundary conditions, with Coulombic electrostatic interactions and in systems large enough to reduce finite-size associated errors in thermodynamic quantities to within the thermal energy, are still not affordable. As a result, periodic boundary conditions, systems of microscopic size and effective electrostatic interaction functions are frequently resorted to. Ensuing artifacts in thermodynamic quantities are nowadays routinely corrected a posteriori, but the underlying configurational sampling still descends from spurious forces. The present study addresses this problem through the introduction of on-the-fly corrections to the physical forces during an atomistic molecular dynamics simulation. Two different approaches are suggested, where the force corrections are derived from special potential energy terms. In the first approach, the solvent-generated electrostatic potential sampled at a given atom site is restrained to a target value involving corrections for electrostatic artifacts. In the second approach, the long-range regime of the solvent polarization around a given atom site is restrained to the Born polarization, i.e., the solvent polarization corresponding to the ideal situation of a macroscopic system under nonperiodic boundary conditions and governed by Coulombic electrostatic interactions. The restraints are applied to the explicit-water simulation of a hydrated sodium ion, and the effect of the restraints on the structural and energetic properties of the solvent is illustrated. Furthermore, by means of the calculation of the charging free energy of a hydrated sodium ion, it is shown how the electrostatic potential restraint translates into the on-the-fly consideration of the corresponding free-energy correction terms. It is discussed how the restraints can be generalized to situations involving several solute particles. Although the present study considers a very simple system only, it is an important step toward the on-the-fly elimination of finite-size and approximate-electrostatic artifacts during atomistic molecular dynamics simulations.     

Electrostatics of proteins in dielectric solvent continua. II. First applications in molecular dynamics simulations

In the preceding paper by Stork and Tavan, [J. Chem. Phys. 126, 165105 (2007)], the authors have reformulated an electrostatic theory which treats proteins surrounded by dielectric solvent continua and approximately solves the associated Poisson equation [B. Egwolf and P. Tavan, J. Chem. Phys. 118, 2039 (2003)]. The resulting solution comprises analytical expressions for the electrostatic reaction field (RF) and potential, which are generated within the protein by the polarization of the surrounding continuum. Here the field and potential are represented in terms of Gaussian RF dipole densities localized at the protein atoms. Quite like in a polarizable force field, also the RF dipole at a given protein atom is induced by the partial charges and RF dipoles at the other atoms. Based on the reformulated theory, the authors have suggested expressions for the RF forces, which obey Newton's third law. Previous continuum approaches, which were also built on solutions of the Poisson equation, used to violate the reactio principle required by this law, and thus were inapplicable to molecular dynamics (MD) simulations. In this paper, the authors suggest a set of techniques by which one can surmount the few remaining hurdles still hampering the application of the theory to MD simulations of soluble proteins and peptides. These techniques comprise the treatment of the RF dipoles within an extended Lagrangian approach and the optimization of the atomic RF polarizabilities. Using the well-studied conformational dynamics of alanine dipeptide as the simplest example, the authors demonstrate the remarkable accuracy and efficiency of the resulting RF-MD approach.     

Treatment of electrostatic interactions in simulations of the interface between transition metal surfaces and organic matter: the discrete classical model and the reaction-field dielectric method

A more efficient molecular dynamics model is devised to study the interface between platinum(111) and liquid isopropanol. The faster reaction-field dielectric method (RFD) is investigated as the alternative to the very accurate discrete classical model (DCM) (M. W. Finnis, Surf. Sci., 1991, 241, 61) to account for the electrostatic interactions at the interface. Differences in static and dynamic properties of the interface are examined. The mass-density and number-density distributions showed only a little sensitivity towards the choice of the method. The orientational structure revealed notable but still small quantitative changes for highly-polarized parts of the isopropanol molecules if the RFD is used instead of the DCM. The general qualitative picture is still preserved. The simulations of the platinum(111)/isopropanol interface with the DCM or the Ewald summation were carried out for longer times than in our previous contribution (K. B. Tarmyshov and F. Müller-Plathe, J. Chem. Phys., 2007, 126, 74702). The analysis disclosed that the influence of chemical or electrostatic interactions at the interface on the relative orientation of the neighbouring hydroxyl groups, which are adsorbed on the surface, is notably smaller than reported before. Overall, the reaction-field dielectric method can be used to study the properties of the interface between transition metals and organic matter, if there appears a chemical bond of 10 kT or higher between the surface and adsorbate molecules or their parts.     

A finite expansion method for the calculation and interpretation of molecular electrostatic potentials

Because it is useful to have the molecular electrostatic potential as an element in a complex scheme to assess the toxicity of large molecules, efficient and reliable methods are needed for the calculation and characterization of these potentials. A multicenter multipole expansion of the molecular electron charge density calculated with a limited Gaussian basis set is shown here to have only a finite number of nonzero terms from which the molecular electrostatic potential can be calculated. The discrete contributions to the electrostatic potentials from the terms of this expansion provide a physically meaningful decomposition of the potential and a means for its characterization. With pyrrole as an example, the electrostatic potential calculated from this finite expansion of the electron density is compared to that obtained from exact calculations from the same wave function. Good agreement is obtained at distances greater than 1.5 A from any atom in the molecule. In contrast, rearrangement of the terms into an expansion corresponding only to Mulliken atomic charges and dipoles yields a decomposition that produces electrostatic potentials which agree less well with the exact potential. This discrepancy is attributable to the neglect of terms due to higher moments.