A fast multipole method combined with a reaction field for long-range electrostatics in molecular dynamics simulations: The effects of truncation on the properties of water

Analytic Gradients for the Electrostatic Embedding QM/MM Model in Periodic Boundary Conditions Using Particle-Mesh Ewald Sums and Electrostatic Potential Fitted Charge Operators

Long-range electrostatic effects are fundamental for describing chemical reactivity in the condensed phase. Here, we present the methodology of an efficient quantum mechanical/molecular mechanical (QM/MM) model in periodic boundary conditions (PBC) compatible with QM/MM boundaries at chemical bonds. The method combines electrostatic potential fitted charge operators and electrostatic potentials derived from the smooth particle-mesh Ewald (PME) sum approach. The total energy and its analytic first derivatives with respect to QM, MM, and lattice vectors allow QM/MM molecular dynamics (MD) in the most common thermodynamic ensembles. We demonstrate the robustness of the method by performing a QM/MM MD equilibration of methanol in water. We simulate the cis/trans isomerization free-energy profiles in water of proline amino acid and a proline-containing oligopeptide, showing a correct description of the reaction barrier. Our PBC-compatible QM/MM model can efficiently be used to study the chemical reactivity in the condensed phase and enzymatic catalysis.

Unsolved problem of long-range interactions: dipolar spin-ice study

Long-range interactions derive various strange phenomena. As illustrated by cutoff simulations of water, increasing cutoff length does not improve the simulation result necessarily; on the contrary, it makes the result worse. In the extreme situation, the structure of water transforms into a layer structure. In this study, to explore the underlying mechanism of this phenomenon, we performed Monte Carlo simulations on dipolar spins arranged on a pyrochlore spin-ice lattice. Like the water case, the present dipolar spin system also showed cutoff-induced dipole ordering and layer formation. The width of the layers depended on the cutoff length; and longer cutoff length led to a broader layer. These features are certainly consistent with the previous water case. This indicates that layer formation is the general behavior of dipolar systems whose interactions are truncated within a finite distance. The result is important for future exploration of the relationship between long-range interactions and resulting structures. In addition, it emphasizes the necessity of rigorous treatment of long-range interactions because increasing the cutoff length prevents convergence and provides an entirely different result from the rigorous Ewald calculation.

Simulations of Amyloid-Forming Peptides in the Crystal State

There still is little treatment available for amyloid diseases, despite their significant impact on individuals and the social and economic implications for society. One reason for this is that the physical nature of amyloid formation is not understood sufficiently well. Therefore, fundamental research at the molecular level remains necessary to support the development of therapeutics. A few structures of short peptides from amyloid-forming proteins have been determined. These can in principle be used as scaffolds for designing aggregation inhibitors. Attempts to this end have often used the tools of computational chemistry, in particular molecular simulation. However, few simulation studies of these peptides in the crystal state have been presented so far. Hence, to validate the capability of common force fields (AMBER19SB, CHARMM36m, and OPLS-AA/M) to yield insight into the dynamics and structural stability of amyloid peptide aggregates, we have performed molecular dynamics simulations of twelve different peptide crystals at two different temperatures. From the simulations, we evaluate the hydrogen bonding patterns, the isotropic B-factors, the change in energy, the Ramachandran plots, and the unit cell parameters and compare the results with the crystal structures. Most crystals are stable in the simulations but for all force fields there is at least one that deviates from the experimental crystal, suggesting more work is needed on these models.

Predicting liposome formulations by the integrated machine learning and molecular modeling approaches

Liposome is one of the most widely used carriers for drug delivery because of the great biocompatibility and biodegradability. Due to the complex formulation components and preparation process, formulation screening mostly relies on trial-and-error process with low efficiency. Here liposome formulation prediction models have been built by machine learning (ML) approaches. The important parameters of liposomes, including size, polydispersity index (PDI), zeta potential and encapsulation, are predicted individually by optimal ML algorithm, while the formulation features are also ranked to provide important guidance for formulation design. The analysis of key parameter reveals that drug molecules with logS [-3, -6], molecular complexity [500, 1000] and XLogP3 (≥2) are priority for preparing liposome with higher encapsulation. In addition, naproxen (NAP) and palmatine HCl (PAL) represented the insoluble and water-soluble molecules are prepared as liposome formulations to validate prediction ability. The consistency between predicted and experimental value verifies the satisfied accuracy of ML models. As the drug properties are critical for liposome particles, the molecular interactions and dynamics of NAP and PAL liposome are further investigated by coarse-grained molecular dynamics simulations. The modeling structure reveals that NAP molecules could distribute into lipid layer, while most PAL molecules aggregate in the inner aqueous phase of liposome. The completely different physical state of NAP and PAL confirms the importance of drug properties for liposome formulations. In summary, the general prediction models are built to predict liposome formulations, and the impacts of key factors are analyzed by combing ML with molecular modeling. The availability and rationality of these intelligent prediction systems have been proved in this study, which could be applied for liposome formulation development in the future.

Molecular simulations of proteins: From simplified physical interactions to complex biological phenomena

Molecular dynamics simulation is the most popular computational technique for investigating the structural and dynamical behaviour of proteins, in search of the molecular basis of their function. Far from being a completely settled field of research, simulations are still evolving to best capture the essential features of the atomic interactions that govern a protein's inner motions. Modern force fields are becoming increasingly accurate in providing a physical description adequate to this purpose, and allow us to model complex biological systems under fairly realistic conditions. Furthermore, the use of accelerated sampling techniques is improving our access to the observation of progressively larger molecular structures, longer time scales, and more hidden functional events. In this review, the basic principles of molecular dynamics simulations and a number of key applications in the area of protein science are summarized, and some of the most important results are discussed. Examples include the study of the structure, dynamics and binding properties of 'difficult' targets, such as intrinsically disordered proteins and membrane receptors, and the investigation of challenging phenomena like hydration-driven processes and protein aggregation. The findings described provide an overall picture of the current state of this research field, and indicate new perspectives on the road ahead to the upcoming future of molecular simulations.

Dielectric continuum model examination of real-space electrostatic treatments

Electrostatic interaction is long ranged; thus, the accurate calculation is not an easy task in molecular dynamics or Monte Carlo simulations. Though the rigorous Ewald method based on the reciprocal space has been established, real-space treatments have recently become an attractive alternative because of the efficient calculation. However, the construction is not yet completed and is now a challenging subject. In an earlier theoretical study, Neumann and Steinhauser employed the Onsager dielectric continuum model to explain how simple real-space cutoff produces artificial dipolar orientation. In the present study, we employ this continuum model to explore the fundamental properties of the recently developed real-space treatments of three shifting schemes. The result of the distance-dependent Kirkwood function G_K(R) showed that the simple bare cutoff produces a well-known hole-shaped artifact, whereas the shift treatments do not. Two-dimensional mapping of electric field well explained how these shift treatments remove the hole-shaped artifact. Still, the shift treatments are not sufficient because they do not produce a flat G_K(R) profile unlike ideal no-cutoff treatment. To test the continuum model results, we also performed Monte Carlo simulations of dipolar particles. The results found that the continuum model could predict the qualitative tendency as to whether each electrostatic treatment produces the hole-shaped artifact of G_K(R) or not. We expect that the present study using the continuum model offers a stringent criterion to judge whether the primitive electrostatic behavior is correctly described or not, which will be useful for future construction of electrostatic treatments.

Molecular Dynamics Simulations Applied to Structural and Dynamical Transitions of the Huntingtin Protein: A Review

Over the recent years, Huntington's disease (HD) has become widely discussed in the scientific literature especially because at the mutant level there are several contradictions regarding the aggregation mechanism. The specific role of the physiological huntingtin protein remains unknown, due to the lack of characterization of its entire crystallographic structure, making the experimental and theoretical research even harder when taking into consideration its involvement in multiple biological functions and its high affinity for different interacting partners. Different types of models, containing fewer (not more than 35 Qs) polyglutamine residues for the WT structure and above 35 Qs for the mutants, were subjected to classical or advanced MD simulations to establish the proteins' structural stability by evaluating their conformational changes. Outside the polyQ tract, there are two other regions of interest (the N17 domain and the polyP rich domain) considered to be essential for the aggregation kinetics at the mutant level. The polymerization process is considered to be dependent on the polyQ length. As the polyQ tract's dimension increases, the structures present more β-sheet conformations. Contrarily, it is also considered that the aggregation stability is not necessarily dependent on the number of Qs, while the initial stage of the aggregation seed might play the decisive role. A general assumption regarding the polyP domain is that it might preserve the polyQ structures soluble by acting as an antagonist for β-sheet formation.

Electrospray molecular dynamics simulations using an octree-based Coulomb interaction method

A new octree-based Coulomb interaction model is developed to model the electrospray of ionic liquids (ILs) in molecular dynamics. Using an octree-based method, Coulomb interactions are categorized as intra- and interleaf Coulomb interactions based on a criterion related to the Bjerrum length of the IL. The octree-based method is found capable of reproducing Coulomb energy in agreement with established and computationally more expensive models, such as the direct Coulomb and the damped shifted force (DSF) method in the absence of an external electric field. In the presence of an external electric field, the octree-based method produces distinctly different results compared to that obtained by the direct Coulomb method. The time required to form Taylor's cone was shorter for the octree method compared to the direct Coulomb approach. While no emission larger than monomers was observed from the direct Coulomb simulation, emission of larger species such as dimers and trimers was observed when the octree-based Coulomb interaction model was used. Furthermore, the octree-based model forms a smaller ion emission cone compared to that from the direct Coulomb method.

A fast and accurate computational method for the linear-combination-based isotropic periodic sum

An isotropic periodic sum (IPS) is a powerful technique to reasonably calculate intermolecular interactions for wide range of molecular systems under periodic boundary conditions. A linear-combination-based IPS (LIPS) has been developed to attain computational accuracy close to an exact lattice sum, such as the Ewald sum. The algorithm of the original LIPS method has a high computational cost because it needs long-range interaction calculations in real space. This becomes a performance bottleneck for long-time molecular simulations. In this work, the combination of an LIPS and fast Fourier transform (FFT) was developed, and evaluated on homogeneous and heterogeneous molecular systems. This combinational approach of LIPS/FFT attained computational efficiency close to that of a smooth particle mesh Ewald while maintaining the same high accuracy as the original LIPS. We concluded that LIPS/FFT has great potential to extend the capability of IPS techniques for the fast and accurate computation of many types of molecular systems.

Anisotropic reaction field correction for long-range electrostatic interactions in molecular dynamics simulations

Reaction-field (RF) methods have been extensively used in molecular dynamics simulations to efficiently compute long-range electrostatic interactions. They assume a continuous dielectric medium outside a certain cutoff, which has shown to be a reasonable approximation in many cases. However, lattice sum or fast multipole methods are nowadays often used instead, which treat long-range interactions explicitly but may introduce different artefacts. In the following work, the major issue of RFs is addressed, i.e., their inability to account for inhomogeneity even in heterogenous environments (e.g., membranes or protein binding pockets). By using a first-order Laplace series expansion of the dielectric permittivity on the cutoff sphere, local anisotropic effects can be described in a simple form. It is shown that the resulting boundary-value problem cannot be solved analytically, but instead a well-behaved approximative anisotropic reaction field (ARF) is introduced, which preserves coordinate invariance and approaches the standard RF solution for homogeneous systems. The comparison of RF to the state-of-the-art particle-particle particle-mesh (P₃M) method shows a difference in the orientation of molecules close to the interface between two different dielectrics (water/chloroform). It was found that the ARF leads to a correction of the orientational distribution toward the P₃M reference for planar and spherical interfaces.

Critical test of isotropic periodic sum techniques with group-based cut-off schemes

Truncation is still chosen for many long-range intermolecular interaction calculations to efficiently compute free-boundary systems, macromolecular systems and net-charge molecular systems, for example. Advanced truncation methods have been developed for long-range intermolecular interactions. Every truncation method can be implemented as one of two basic cut-off schemes, namely either an atom-based or a group-based cut-off scheme. The former computes interactions of “atoms” inside the cut-off radius, whereas the latter computes interactions of “molecules” inside the cut-off radius. In this work, the effect of group-based cut-off is investigated for isotropic periodic sum (IPS) techniques, which are promising cut-off treatments to attain advanced accuracy for many types of molecular system. The effect of group-based cut-off is clearly different from that of atom-based cut-off, and severe artefacts are observed in some cases. However, no severe discrepancy from the Ewald sum is observed with the extended IPS techniques.

Are the Sublimation Thermodynamics of Organic Molecules Predictable?

We compare a range of computational methods for the prediction of sublimation thermodynamics (enthalpy, entropy, and free energy of sublimation). These include a model from theoretical chemistry that utilizes crystal lattice energy minimization (with the DMACRYS program) and quantitative structure property relationship (QSPR) models generated by both machine learning (random forest and support vector machines) and regression (partial least squares) methods. Using these methods we investigate the predictability of the enthalpy, entropy and free energy of sublimation, with consideration of whether such a method may be able to improve solubility prediction schemes. Previous work has suggested that the major source of error in solubility prediction schemes involving a thermodynamic cycle via the solid state is in the modeling of the free energy change away from the solid state. Yet contrary to this conclusion other work has found that the inclusion of terms such as the enthalpy of sublimation in QSPR methods does not improve the predictions of solubility. We suggest the use of theoretical chemistry terms, detailed explicitly in the Methods section, as descriptors for the prediction of the enthalpy and free energy of sublimation. A data set of 158 molecules with experimental sublimation thermodynamics values and some CSD refcodes has been collected from the literature and is provided with their original source references.

A polarizable QM/MM approach to the molecular dynamics of amide groups solvated in water

The infrared (IR) spectra of polypeptides are dominated by the so-called amide bands. Because they originate from the strongly polar and polarizable amide groups (AGs) making up the backbone, their spectral positions sensitively depend on the local electric fields. Aiming at accurate computations of these IR spectra by molecular dynamics (MD) simulations, which derive atomic forces from a hybrid quantum and molecular mechanics (QM/MM) Hamiltonian, here we consider the effects of solvation in bulk liquid water on the amide bands of the AG model compound N-methyl-acetamide (NMA). As QM approach to NMA we choose grid-based density functional theory (DFT). For the surrounding MM water, we develop, largely based on computations, a polarizable molecular mechanics (PMM) model potential called GP6P, which features six Gaussian electrostatic sources (one induced dipole, five static partial charge distributions) and, therefore, avoids spurious distortions of the DFT electron density in hybrid DFT/PMM simulations. Bulk liquid GP6P is shown to have favorable properties at the thermodynamic conditions of the parameterization and beyond. Lennard-Jones (LJ) parameters of the DFT fragment NMA are optimized by comparing radial distribution functions in the surrounding GP6P liquid with reference data obtained from a "first-principles" DFT-MD simulation. Finally, IR spectra of NMA in GP6P water are calculated from extended DFT/PMM-MD trajectories, in which the NMA is treated by three different DFT functionals (BP, BLYP, B3LYP). Method-specific frequency scaling factors are derived from DFT-MD simulations of isolated NMA. The DFT/PMM-MD simulations with GP6P and with the optimized LJ parameters then excellently predict the effects of aqueous solvation and deuteration observed in the IR spectra of NMA. As a result, the methods required to accurately compute such spectra by DFT/PMM-MD also for larger peptides in aqueous solution are now at hand.

Scaling of Multimillion-Atom Biological Molecular Dynamics Simulation on a Petascale Supercomputer

A strategy is described for a fast all-atom molecular dynamics simulation of multimillion-atom biological systems on massively parallel supercomputers. The strategy is developed using benchmark systems of particular interest to bioenergy research, comprising models of cellulose and lignocellulosic biomass in an aqueous solution. The approach involves using the reaction field (RF) method for the computation of long-range electrostatic interactions, which permits efficient scaling on many thousands of cores. Although the range of applicability of the RF method for biomolecular systems remains to be demonstrated, for the benchmark systems the use of the RF produces molecular dipole moments, Kirkwood G factors, other structural properties, and mean-square fluctuations in excellent agreement with those obtained with the commonly used Particle Mesh Ewald method. With RF, three million- and five million-atom biological systems scale well up to ∼30k cores, producing ∼30 ns/day. Atomistic simulations of very large systems for time scales approaching the microsecond would, therefore, appear now to be within reach.

The "Hot-Solvent/Cold-Solute" Problem Revisited

The temperature steers the equilibrium and nonequilibrium conformational dynamics of macromolecules in solution. Therefore, corresponding molecular dynamics simulations require a strategy for temperature control which should guarantee that the experimental statistical ensemble is also sampled in silico. Several algorithms for temperature control have been proposed. All these thermostats interfere with the macromolecule's "natural" dynamics as given by the Newtonian mechanics. Furthermore, using a single thermostat for an inhomogeneous solute-solvent system can lead to stationary temperature gradients. To avoid this "hot solvent/cold solute" problem, two separate thermostats are frequently applied, one to the solute and one to the solvent. However, such a separate temperature control will perturb the dynamics of the macromolecule much more strongly than a global one and, therefore, can introduce large artifacts into its conformational dynamics. Based on the concept that an explicit solvent environment represents an ideal thermostat concerning the magnitude and time correlation of temperature fluctuations of the solute, we propose a temperature control strategy that, on the one hand, provides a homogeneous temperature distribution throughout the system together with the correct statistical ensemble for the solute molecule while, on the other hand, minimally perturbing its dynamics.

Simulated Solute Tempering

For the enhanced conformational sampling in molecular dynamics (MD) simulations, we present "simulated solute tempering" (SST) which is an easy to implement variant of simulated tempering. SST extends conventional simulated tempering (CST) by key concepts of "replica exchange with solute tempering" (REST, Liu et al. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13749). We have applied SST, CST, and REST to molecular dynamics (MD) simulations of an alanine octapeptide in explicit water. The weight parameters required for CST and SST are determined by two different formulas whose performance is compared. For SST only one of them yields a uniform sampling of the temperature space. Compared to CST and REST, SST provides the highest exchange probabilities between neighboring rungs in the temperature ladder. Concomitantly, SST leads to the fastest diffusion of the simulation system through the temperature space, in particular, if the "even-odd" exchange scheme is employed in SST. As a result, SST exhibits the highest sampling speed of the investigated tempering methods.

Linearly scaling and almost Hamiltonian dielectric continuum molecular dynamics simulations through fast multipole expansions

Hamiltonian Dielectric Solvent (HADES) is a recent method [S. Bauer et al., J. Chem. Phys. 140, 104103 (2014)] which enables atomistic Hamiltonian molecular dynamics (MD) simulations of peptides and proteins in dielectric solvent continua. Such simulations become rapidly impractical for large proteins, because the computational effort of HADES scales quadratically with the number N of atoms. If one tries to achieve linear scaling by applying a fast multipole method (FMM) to the computation of the HADES electrostatics, the Hamiltonian character (conservation of total energy, linear, and angular momenta) may get lost. Here, we show that the Hamiltonian character of HADES can be almost completely preserved, if the structure-adapted fast multipole method (SAMM) as recently redesigned by Lorenzen et al. [J. Chem. Theory Comput. 10, 3244-3259 (2014)] is suitably extended and is chosen as the FMM module. By this extension, the HADES/SAMM forces become exact gradients of the HADES/SAMM energy. Their translational and rotational invariance then guarantees (within the limits of numerical accuracy) the exact conservation of the linear and angular momenta. Also, the total energy is essentially conserved-up to residual algorithmic noise, which is caused by the periodically repeated SAMM interaction list updates. These updates entail very small temporal discontinuities of the force description, because the employed SAMM approximations represent deliberately balanced compromises between accuracy and efficiency. The energy-gradient corrected version of SAMM can also be applied, of course, to MD simulations of all-atom solvent-solute systems enclosed by periodic boundary conditions. However, as we demonstrate in passing, this choice does not offer any serious advantages.

Parameterization of the Hamiltonian Dielectric Solvent (HADES) Reaction-Field Method for the Solvation Free Energies of Amino Acid Side-Chain Analogs

Optimization of the Hamiltonian dielectric solvent (HADES) method for biomolecular simulations in a dielectric continuum is presented with the goal of calculating accurate absolute solvation free energies while retaining the model's accuracy in predicting conformational free-energy differences. The solvation free energies of neutral and polar amino acid side-chain analogs calculated by using HADES, which may optionally include nonpolar contributions, were optimized against experimental data to reach a chemical accuracy of about 0.5 kcal mol(-1). The new parameters were evaluated for charged side-chain analogs. The HADES results were compared with explicit-solvent, generalized Born, Poisson-Boltzmann, and QM-based methods. The potentials of mean force (PMFs) between pairs of side-chain analogs obtained by using HADES and explicit-solvent simulations were used to evaluate the effects of the improved parameters optimized for solvation free energies on intermolecular potentials.

Utilizing fast multipole expansions for efficient and accurate quantum-classical molecular dynamics simulations

Recently, a novel approach to hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations has been suggested [Schwörer et al., J. Chem. Phys. 138, 244103 (2013)]. Here, the forces acting on the atoms are calculated by grid-based density functional theory (DFT) for a solute molecule and by a polarizable molecular mechanics (PMM) force field for a large solvent environment composed of several 10(3)-10(5) molecules as negative gradients of a DFT/PMM hybrid Hamiltonian. The electrostatic interactions are efficiently described by a hierarchical fast multipole method (FMM). Adopting recent progress of this FMM technique [Lorenzen et al., J. Chem. Theory Comput. 10, 3244 (2014)], which particularly entails a strictly linear scaling of the computational effort with the system size, and adapting this revised FMM approach to the computation of the interactions between the DFT and PMM fragments of a simulation system, here, we show how one can further enhance the efficiency and accuracy of such DFT/PMM-MD simulations. The resulting gain of total performance, as measured for alanine dipeptide (DFT) embedded in water (PMM) by the product of the gains in efficiency and accuracy, amounts to about one order of magnitude. We also demonstrate that the jointly parallelized implementation of the DFT and PMM-MD parts of the computation enables the efficient use of high-performance computing systems. The associated software is available online.

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.

The zero-multipole summation method for estimating electrostatic interactions in molecular dynamics: analysis of the accuracy and application to liquid systems

In the preceding paper [I. Fukuda, J. Chem. Phys. 139, 174107 (2013)], the zero-multipole (ZM) summation method was proposed for efficiently evaluating the electrostatic Coulombic interactions of a classical point charge system. The summation takes a simple pairwise form, but prevents the electrically non-neutral multipole states that may artificially be generated by a simple cutoff truncation, which often causes large energetic noises and significant artifacts. The purpose of this paper is to judge the ability of the ZM method by investigating the accuracy, parameter dependencies, and stability in applications to liquid systems. To conduct this, first, the energy-functional error was divided into three terms and each term was analyzed by a theoretical error-bound estimation. This estimation gave us a clear basis of the discussions on the numerical investigations. It also gave a new viewpoint between the excess energy error and the damping effect by the damping parameter. Second, with the aid of these analyses, the ZM method was evaluated based on molecular dynamics (MD) simulations of two fundamental liquid systems, a molten sodium-chlorine ion system and a pure water molecule system. In the ion system, the energy accuracy, compared with the Ewald summation, was better for a larger value of multipole moment l currently induced until l ≲ 3 on average. This accuracy improvement with increasing l is due to the enhancement of the excess-energy accuracy. However, this improvement is wholly effective in the total accuracy if the theoretical moment l is smaller than or equal to a system intrinsic moment L. The simulation results thus indicate L ∼ 3 in this system, and we observed less accuracy in l = 4. We demonstrated the origins of parameter dependencies appearing in the crossing behavior and the oscillations of the energy error curves. With raising the moment l we observed, smaller values of the damping parameter provided more accurate results and smoother behaviors with respect to cutoff length were obtained. These features can be explained, on the basis of the theoretical error analyses, such that the excess energy accuracy is improved with increasing l and that the total accuracy improvement within l ⩽ L is facilitated by a small damping parameter. Although the accuracy was fundamentally similar to the ion system, the bulk water system exhibited distinguishable quantitative behaviors. A smaller damping parameter was effective in all the practical cutoff distance, and this fact can be interpreted by the reduction of the excess subset. A lower moment was advantageous in the energy accuracy, where l = 1 was slightly superior to l = 2 in this system. However, the method with l = 2 (viz., the zero-quadrupole sum) gave accurate results for the radial distribution function. We confirmed the stability in the numerical integration for MD simulations employing the ZM scheme. This result is supported by the sufficient smoothness of the energy function. Along with the smoothness, the pairwise feature and the allowance of the atom-based cutoff mode on the energy formula lead to the exact zero total-force, ensuring the total-momentum conservations for typical MD equations of motion.

Electrostatics of proteins in dielectric solvent continua. I. An accurate and efficient reaction field description

We present a reaction field (RF) method which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of an electrostatics calculation with polarizable molecular mechanics (MM) force fields. The method combines an approach originally suggested by Egwolf and Tavan [J. Chem. Phys. 118, 2039 (2003)] with concepts generalizing the Born solution [Z. Phys. 1, 45 (1920)] for a solvated ion. First, we derive an exact representation according to which the sources of the RF potential and energy are inducible atomic anti-polarization densities and atomic shielding charge distributions. Modeling these atomic densities by Gaussians leads to an approximate representation. Here, the strengths of the Gaussian shielding charge distributions are directly given in terms of the static partial charges as defined, e.g., by standard MM force fields for the various atom types, whereas the strengths of the Gaussian anti-polarization densities are calculated by a self-consistency iteration. The atomic volumes are also described by Gaussians. To account for covalently overlapping atoms, their effective volumes are calculated by another self-consistency procedure, which guarantees that the dielectric function ε(r) is close to one everywhere inside the protein. The Gaussian widths σ(i) of the atoms i are parameters of the RF approximation. The remarkable accuracy of the method is demonstrated by comparison with Kirkwood's analytical solution for a spherical protein [J. Chem. Phys. 2, 351 (1934)] and with computationally expensive grid-based numerical solutions for simple model systems in dielectric continua including a di-peptide (Ac-Ala-NHMe) as modeled by a standard MM force field. The latter example shows how weakly the RF conformational free energy landscape depends on the parameters σ(i). A summarizing discussion highlights the achievements of the new theory and of its approximate solution particularly by comparison with so-called generalized Born methods. A follow-up paper describes how the method enables Hamiltonian, efficient, and accurate MM molecular dynamics simulations of proteins in dielectric solvent continua.

Electrostatics of proteins in dielectric solvent continua. II. Hamiltonian reaction field dynamics

In Paper I of this work [S. Bauer, G. Mathias, and P. Tavan, J. Chem. Phys. 140, 104102 (2014)] we have presented a reaction field (RF) method, which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of polarizable molecular mechanics (MM) force fields. Building upon these results, here we suggest a method for linearly scaling Hamiltonian RF/MM molecular dynamics (MD) simulations, which we call "Hamiltonian dielectric solvent" (HADES). First, we derive analytical expressions for the RF forces acting on the solute atoms. These forces properly account for all those conditions, which have to be self-consistently fulfilled by RF quantities introduced in Paper I. Next we provide details on the implementation, i.e., we show how our RF approach is combined with a fast multipole method and how the self-consistency iterations are accelerated by the use of the so-called direct inversion in the iterative subspace. Finally we demonstrate that the method and its implementation enable Hamiltonian, i.e., energy and momentum conserving HADES-MD, and compare in a sample application on Ac-Ala-NHMe the HADES-MD free energy landscape at 300 K with that obtained in Paper I by scanning of configurations and with one obtained from an explicit solvent simulation.

Design of a reaction field using a linear-combination-based isotropic periodic sum method

In our previous study (Takahashi et al., J. Chem. Theory Comput. 2012, 8, 4503), we developed the linear-combination-based isotropic periodic sum (LIPS) method. The LIPS method is based on the extended isotropic periodic sum theory that produces a ubiquitous interaction potential function to estimate homogeneous and heterogeneous systems. The LIPS theory also provides the procedure to design a periodic reaction field. To demonstrate this, in the present work, a novel reaction field of the LIPS method was developed. The novel reaction field was labeled LIPS-SW, because it provides an interaction potential function with a shape that resembles that of the switch function method. To evaluate the ability of the LIPS-SW method to describe in homogeneous and heterogeneous systems, we carried out molecular dynamics (MD) simulations of bulk water and water-vapor interfacial systems using the LIPS-SW method. The results of these simulations show that the LIPS-SW method gives higher accuracy than the conventional interaction potential function of the LIPS method. The accuracy of simulating water-vapor interfacial systems was greatly improved, while that of bulk water systems was maintained using the LIPS-SW method. We conclude that the LIPS-SW method shows great potential for high-accuracy, high-performance computing to allow large scale MD simulations. © 2014 Wiley Periodicals, Inc.

Coupling density functional theory to polarizable force fields for efficient and accurate Hamiltonian molecular dynamics simulations

Hybrid molecular dynamics (MD) simulations, in which the forces acting on the atoms are calculated by grid-based density functional theory (DFT) for a solute molecule and by a polarizable molecular mechanics (PMM) force field for a large solvent environment composed of several 10(3)-10(5) molecules, pose a challenge. A corresponding computational approach should guarantee energy conservation, exclude artificial distortions of the electron density at the interface between the DFT and PMM fragments, and should treat the long-range electrostatic interactions within the hybrid simulation system in a linearly scaling fashion. Here we describe a corresponding Hamiltonian DFT/(P)MM implementation, which accounts for inducible atomic dipoles of a PMM environment in a joint DFT/PMM self-consistency iteration. The long-range parts of the electrostatics are treated by hierarchically nested fast multipole expansions up to a maximum distance dictated by the minimum image convention of toroidal boundary conditions and, beyond that distance, by a reaction field approach such that the computation scales linearly with the number of PMM atoms. Short-range over-polarization artifacts are excluded by using Gaussian inducible dipoles throughout the system and Gaussian partial charges in the PMM region close to the DFT fragment. The Hamiltonian character, the stability, and efficiency of the implementation are investigated by hybrid DFT/PMM-MD simulations treating one molecule of the water dimer and of bulk water by DFT and the respective remainder by PMM.

Polarizable six-point water models from computational and empirical optimization

Tröster et al. (J. Phys. Chem B 2013, 117, 9486-9500) recently suggested a mixed computational and empirical approach to the optimization of polarizable molecular mechanics (PMM) water models. In the empirical part the parameters of Buckingham potentials are optimized by PMM molecular dynamics (MD) simulations. The computational part applies hybrid calculations, which combine the quantum mechanical description of a H2O molecule by density functional theory (DFT) with a PMM model of its liquid phase environment generated by MD. While the static dipole moments and polarizabilities of the PMM water models are fixed at the experimental gas phase values, the DFT/PMM calculations are employed to optimize the remaining electrostatic properties. These properties cover the width of a Gaussian inducible dipole positioned at the oxygen and the locations of massless negative charge points within the molecule (the positive charges are attached to the hydrogens). The authors considered the cases of one and two negative charges rendering the PMM four- and five-point models TL4P and TL5P. Here we extend their approach to three negative charges, thus suggesting the PMM six-point model TL6P. As compared to the predecessors and to other PMM models, which also exhibit partial charges at fixed positions, TL6P turned out to predict all studied properties of liquid water at p0 = 1 bar and T0 = 300 K with a remarkable accuracy. These properties cover, for instance, the diffusion constant, viscosity, isobaric heat capacity, isothermal compressibility, dielectric constant, density, and the isobaric thermal expansion coefficient. This success concurrently provides a microscopic physical explanation of corresponding shortcomings of previous models. It uniquely assigns the failures of previous models to substantial inaccuracies in the description of the higher electrostatic multipole moments of liquid phase water molecules. Resulting favorable properties concerning the transferability to other temperatures and conditions like the melting of ice are also discussed.

The Microscopic Physical Cause for the Density Maximum of Liquid Water

The existence of a density maximum at 277 K is probably the most prominent anomaly among the many very special thermodynamic properties of liquid water. While usually attributed to so-called hydrogen bonding, the microscopic physical cause of this prominent anomaly is still elusive. Here we show that the density anomaly is caused by those short-range electrostatic forces, which are generated by the quadrupole and higher moments of the charge distributions present in liquid-phase water molecules. This conclusion derives from 20 ns replica exchange molecular-dynamics simulations with closely related polarizable four-, five-, and six-point water models. As soon as the model complexity suffices to represent the higher electrostatic moments with sufficient accuracy, the density temperature profile n(T) calculated for T ∈ [250,320] K at the standard pressure 1 bar locks in to the experimental observation. The corresponding six-point model is, therefore, the most simple available cartoon for liquid-phase water molecules.

Polarizable water models from mixed computational and empirical optimization

Here we suggest a mixed computational and empirical approach serving to optimize the parameters of complex and polarizable molecular mechanics (PMM) models for complicated liquids. The computational part of the parameter optimization relies on hybrid calculations combining density functional theory (DFT) for a solute molecule with a PMM treatment of its solvent environment at well-defined thermodynamic conditions. As an application we have developed PMM models for water featuring ν = 3, 4, and 5 points of force action, a Gaussian inducible dipole and a Buckingham potential at the oxygen, the experimental liquid phase geometry, the experimental gas phase polarizability α(exp)(g) = 1.47 ų, and, for ν = 4 and 5, the gas phase value μ(exp)(g) = 1.855 D for the static dipole moment. The widths of the Gaussian dipoles and, for ν = 4 and 5, also the electrostatic geometries of these so-called TLνP models are derived from self-consistent DFT/PMM calculations, and the parameters of the Buckingham potentials (and the static TL3P dipole moment) are estimated from molecular dynamics (MD) simulations. The high quality of the resulting models is demonstrated for the observables targeted during optimization (potential energy per molecule, pressure, radial distribution functions) and a series of predicted properties (quadrupole moments, density at constant pressure, dielectric constant, diffusivity, viscosity, compressibility, heat capacity) at certain standard conditions. Remaining deficiencies and possible ways for their removal are discussed.

Non-Ewald methods: theory and applications to molecular systems

Several non-Ewald methods for calculating electrostatic interactions have recently been developed, such as the Wolf method, the reaction field method, the pre-averaging method, and the zero-dipole summation method, for molecular dynamics simulations of various physical systems, including biomolecular systems. We review the theories of these approaches and their potential applications to molecular simulations, and discuss their relationships.

Optimizing the Accuracy and Efficiency of Fast Hierarchical Multipole Expansions for MD Simulations

Based on p'th order Cartesian Taylor expansions of Coulomb interactions and on hierarchical decompositions of macromolecular simulation systems into hierarchies of nested, structure-adapted, and adaptively formed clusters of increasing size, fast multipole methods are constructed for rapid and accurate calculations of electrostatic interactions. These so-called SAMMp algorithms are formulated through totally symmetric and traceless tensors describing the multipole moments and the coefficients of local Taylor expansions. Simple recursions for the efficient evaluation and shifting of multipole moments are given. The required tensors are explicitly given up to order p = 4. The SAMMp algorithms are shown to guarantee the reaction principle. For systems with periodic boundaries, a reaction field (RF) correction is applied, which introduces at distances beyond the "minimum image convention" boundary a dielectric continuum surrounding each cluster at the top level of coarse graining. The correctness of the present SAMMp implementation is demonstrated by analyzing the scaling of the residuals and by checking the numerical accuracy of the reaction principle for a pair of distant molecular ions in vacuum. Molecular dynamics simulations of pure water and aqueous solutions containing artificial ions, which are enclosed by periodic boundaries, demonstrate the stability and low-noise behavior of SAMMp/RF.

The ELBA force field for coarse-grain modeling of lipid membranes

A new coarse-grain model for molecular dynamics simulation of lipid membranes is presented. Following a simple and conventional approach, lipid molecules are modeled by spherical sites, each representing a group of several atoms. In contrast to common coarse-grain methods, two original (interdependent) features are here adopted. First, the main electrostatics are modeled explicitly by charges and dipoles, which interact realistically through a relative dielectric constant of unity (ε(r) = 1). Second, water molecules are represented individually through a new parametrization of the simple Stockmayer potential for polar fluids; each water molecule is therefore described by a single spherical site embedded with a point dipole. The force field is shown to accurately reproduce the main physical properties of single-species phospholipid bilayers comprising dioleoylphosphatidylcholine (DOPC) and dioleoylphosphatidylethanolamine (DOPE) in the liquid crystal phase, as well as distearoylphosphatidylcholine (DSPC) in the liquid crystal and gel phases. Insights are presented into fundamental properties and phenomena that can be difficult or impossible to study with alternative computational or experimental methods. For example, we investigate the internal pressure distribution, dipole potential, lipid diffusion, and spontaneous self-assembly. Simulations lasting up to 1.5 microseconds were conducted for systems of different sizes (128, 512 and 1058 lipids); this also allowed us to identify size-dependent artifacts that are expected to affect membrane simulations in general. Future extensions and applications are discussed, particularly in relation to the methodology's inherent multiscale capabilities.

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 Hü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.

DFT/MM description of flavin IR spectra in BLUF domains

A class of photoreceptors occurring in various organisms consists of domains that are blue light sensing using flavin (BLUF). The vibrational spectra of the flavin chromophore are spectroscopically well characterized for the dark-adapted resting states and for the light-adapted signaling states of BLUF domains in solution. Here we present a theoretical analysis of such spectra by applying density functional theory (DFT) to the flavin embedded in molecular mechanics (MM) models of its protein and solvent environment. By DFT/MM we calculate flavin spectra for seven different X-ray and NMR structures of BLUF domains occurring in the transcriptional antirepressor AppA and in the blue light receptor B (BlrB) of the purple bacterium Rb. Sphaeroides as well as in the phototaxis photoreceptor Slr1694 of the cyanobacterium Synechocystis. By considering the dynamical stabilities of associated all-atom simulation models and by comparing calculated with observed vibrational spectra, we show that two of the considered structures (both AppA) are obviously erroneous and that specific features of two further crystal structures (BlrB and Slr1694) cannot represent the states of the respective BLUF domains in solution. Thereby, the conformational transitions elicited by solvation are identified. In this context we demonstrate how hydrogen bonds of varying strengths can tune in BLUF domains the C═O stretching frequencies of the flavin chromophore. Furthermore we show that the DFT/MM spectra of the flavin calculated for two different AppA BLUF conformations, which are called Trp(in) and Met(in), fit very well to the spectroscopic data observed for the dark and light states, respectively, if (i) polarized MM force fields, which are calculated by an iterative DFT/MM procedure, are employed for the flavin binding pockets and (ii) the calculated frequencies are properly scaled. Although the associated analysis indicates that the Trp(in) conformation belongs to the dark state, no clear light vs dark distinction emerges for the Met(in) conformation. In this connection, a number of methodological issues relevant for such complex computations are thoroughly discussed showing, in particular, why our current descriptions could not decide the light vs dark question for Met(in).

Highly polar environments catalyze the unfolding of PrP C helix 1

The first alpha-helix (H1) likely plays an important role in the conversion of the cellular prion protein (PrP(C)) into its pathogenic isoform (PrP(Sc)). In this conversion, H1 may either have to unfold or may represent a site of intermolecular contact. A recent molecular dynamics simulation suggested that H1 can unfold if it is detached from the protein core (Hirschberger et al. in Biophys J 90:3908, 2006). It has been hypothesized that the high dielectric constant epsilon (S) of the bulk water environment facilitates the unfolding of H1. To check this hypothesis, we performed a number of replica exchange molecular dynamics simulations of an H1 peptide in solvents of different epsilon (S). We found that the equilibrium helix fraction in water is less than 40%, in agreement with previous experimental findings, and that the helix unfolds much faster in water than in less polar solvents. The kinetically stabilizing effect of the organic solvents is largely unspecific and correlates well with their dielectric constant epsilon (S).