The multiscale coarse-graining method. X. Improved algorithms for constructing coarse-grained potentials for molecular systems

The multiscale coarse-graining (MS-CG) method uses simulation data for an atomistic model of a system to construct a coarse-grained (CG) potential for a coarse-grained model of the system. The CG potential is a variational approximation for the true potential of mean force of the degrees of freedom retained in the CG model. The variational calculation uses information about the atomistic positions and forces in the simulation data. In principle, the resulting MS-CG potential will be an accurate representation of the true CG potential if the basis set for the variational calculation is complete enough and the canonical distribution of atomistic states is well sampled by the data set. In practice, atomistic configurations that have very high potential energy are not sampled. As a result there usually is a region of CG configuration space that is not sampled and about which the data set contains no information regarding the gradient of the true potential. The MS-CG potential obtained from a variational calculation will not necessarily be accurate in this unsampled region. A priori considerations make it clear that the true CG potential of mean force must be very large and positive in that region. To obtain an MS-CG potential whose behavior in the sampled region is determined by the atomistic data set, and whose behavior in the unsampled region is large and positive, it is necessary to intervene in the variational calculation in some way. In this paper, we discuss and compare two such methods of intervention, which have been used in previous MS-CG calculations for dealing with nonbonded interactions. For the test systems studied, the two methods give similar results and yield MS-CG potentials that are limited in accuracy only by the incompleteness of the basis set and the statistical error of associated with the set of atomistic configurations used. The use of such methods is important for obtaining accurate CG potentials.     

The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models

Coarse-grained (CG) models provide a computationally efficient method for rapidly investigating the long time- and length-scale processes that play a critical role in many important biological and soft matter processes. Recently, Izvekov and Voth introduced a new multiscale coarse-graining (MS-CG) method [J. Phys. Chem. B 109, 2469 (2005); J. Chem. Phys. 123, 134105 (2005)] for determining the effective interactions between CG sites using information from simulations of atomically detailed models. The present work develops a formal statistical mechanical framework for the MS-CG method and demonstrates that the variational principle underlying the method may, in principle, be employed to determine the many-body potential of mean force (PMF) that governs the equilibrium distribution of positions of the CG sites for the MS-CG models. A CG model that employs such a PMF as a "potential energy function" will generate an equilibrium probability distribution of CG sites that is consistent with the atomically detailed model from which the PMF is derived. Consequently, the MS-CG method provides a formal multiscale bridge rigorously connecting the equilibrium ensembles generated with atomistic and CG models. The variational principle also suggests a class of practical algorithms for calculating approximations to this many-body PMF that are optimal. These algorithms use computer simulation data from the atomically detailed model. Finally, important generalizations of the MS-CG method are introduced for treating systems with rigid intramolecular constraints and for developing CG models whose equilibrium momentum distribution is consistent with that of an atomically detailed model.     

The multiscale coarse-graining method. VII. Free energy decomposition of coarse-grained effective potentials

The potential of mean force (PMF) with respect to coarse-grained (CG) coordinates is often calculated in order to study the molecular interactions in atomistic molecular dynamics (MD) simulations. The multiscale coarse-graining (MS-CG) approach enables the computation of the many-body PMF of an atomistic system in terms of the CG coordinates, which can be used to parameterize CG models based on all-atom configurations. We demonstrate here that the MS-CG method can also be used to analyze the CG interactions from atomistic MD trajectories via PMF calculations. In addition, MS-CG calculations at different temperatures are performed to decompose the PMF values into energetic and entropic contributions as a function of the CG coordinates, which provides more thermodynamic information regarding the atomistic system. Two numerical examples, liquid methanol and a dimyristoylphosphatidylcholine lipid bilayer, are presented. The results show that MS-CG can be used as an analysis tool, comparable to various free energy computation methods. The differences between the MS-CG approach and other PMF calculation methods, as well as the characteristics and advantages of MS-CG, are also discussed.     

The multiscale coarse-graining method. IX. A general method for construction of three body coarse-grained force fields

The multiscale coarse-graining (MS-CG) method is a method for constructing a coarse-grained (CG) model of a system using data obtained from molecular dynamics simulations of the corresponding atomically detailed model. The formal statistical mechanical derivation of the method shows that the potential energy function extracted from an MS-CG calculation is a variational approximation for the true potential of mean force of the CG sites, one that becomes exact in the limit that a complete basis set is used in the variational calculation if enough data are obtained from the atomistic simulations. Most applications of the MS-CG method have employed a representation for the nonbonded part of the CG potential that is a sum of all possible pair interactions. This approach, despite being quite successful for some CG models, is inadequate for some others. Here we propose a systematic method for including three body terms as well as two body terms in the nonbonded part of the CG potential energy. The current method is more general than a previous version presented in a recent paper of this series [L. Larini, L. Lu, and G. A. Voth, J. Chem. Phys. 132, 164107 (2010)], in the sense that it does not make any restrictive choices for the functional form of the three body potential. We use hierarchical multiresolution functions that are similar to wavelets to develop very flexible basis function expansions with both two and three body basis functions. The variational problem is solved by a numerical technique that is capable of automatically selecting an appropriate subset of basis functions from a large initial set. We apply the method to two very different coarse-grained models: a solvent free model of a two component solution made of identical Lennard-Jones particles and a one site model of SPC/E water where a site is placed at the center of mass of each water molecule. These calculations show that the inclusion of three body terms in the nonbonded CG potential can lead to significant improvement in the accuracy of CG potentials and hence of CG simulations.     

A multiscale coarse-graining method for biomolecular systems

A new approach is presented for obtaining coarse-grained (CG) force fields from fully atomistic molecular dynamics (MD) trajectories. The method is demonstrated by applying it to derive a CG model for the dimyristoylphosphatidylcholine (DMPC) lipid bilayer. The coarse-graining of the interparticle force field is accomplished by an application of a force-matching procedure to the force data obtained from an explicit atomistic MD simulation of the biomolecular system of interest. Hence, the method is termed a "multiscale" CG (MS-CG) approach in which explicit atomistic-level forces are propagated upward in scale to the coarse-grained level. The CG sites in the lipid bilayer application were associated with the centers-of-mass of atomic groups because of the simplicity in the evaluation of the forces acting on them from the atomistic data. The resulting CG lipid bilayer model is shown to accurately reproduce the structural properties of the phospholipid bilayer.     

The multiscale coarse-graining method. VIII. Multiresolution hierarchical basis functions and basis function selection in the construction of coarse-grained force fields

The multiscale coarse-graining (MS-CG) method is a method for determining the effective potential energy function for a coarse-grained (CG) model of a molecular system using data obtained from molecular dynamics simulation of the corresponding atomically detailed model. The coarse-grained potential obtained using the MS-CG method is a variational approximation for the exact many-body potential of mean force for the coarse-grained sites. Here we propose a new numerical algorithm with noise suppression capabilities and enhanced numerical stability for the solution of the MS-CG variational problem. The new method, which is a variant of the elastic net method [Friedman et al., Ann. Appl. Stat. 1, 302 (2007)], allows us to construct a large basis set, and for each value of a so-called "penalty parameter" the method automatically chooses a subset of the basis that is most important for representing the MS-CG potential. The size of the subset increases as the penalty parameter is decreased. The appropriate value to choose for the penalty parameter is the one that gives a basis set that is large enough to fit the data in the simulation data set without fitting the noise. This procedure provides regularization to mitigate potential numerical problems in the associated linear least squares calculation, and it provides a way to avoid fitting statistical error. We also develop new basis functions that are similar to multiresolution Haar functions and that have the differentiability properties that are appropriate for representing CG potentials. We demonstrate the feasibility of the combined use of the elastic net method and the multiresolution basis functions by performing a variational calculation of the CG potential for a relatively simple system. We develop a method to choose the appropriate value of the penalty parameter to give the optimal basis set. The combined effect of the new basis functions and the regularization provided by the elastic net method opens the possibility of using very large basis sets for complicated CG systems with many interaction potentials without encountering numerical problems in the variational calculation.     

Peptide folding using multiscale coarse-grained models

The multiscale coarse-graining (MS-CG) method has been previously used to describe the equilibrium properties of peptides. The present study reveals that MS-CG models of alpha-helical polyalanine and the beta-hairpin V 5PGV 5 possess the capacity to efficiently refold in simulations initiated from unfolded configurations. The MS-CG peptides exhibit free energy landscapes that are funneled toward folded configurations and two-state folding behavior, consistent with the known characteristics of small, rapidly folding peptides. Moreover, the models demonstrate enhanced sampling capabilities when compared to systems with full atomic detail. The significance of these observations with respect to the theoretical basis of the MS-CG approach is discussed. The MS-CG peptides were used to reconstruct atomically detailed configurations in order to evaluate the extent to which MS-CG ensembles embody all-atom peptide free energy landscapes. Ensembles obtained from these reconstructed configurations display good agreement with the all-atom simulation data used to generate the MS-CG models and also corroborate the presence of features observed in the MS-CG peptide free energy landscapes. These findings suggest that MS-CG models may be of significant utility in the study of peptide folding.     

Non-Hamiltonian molecular dynamics implementation of the Gibbs ensemble method. II. Molecular liquid-vapor results for carbon dioxide

The Gibbs ensemble molecular dynamics algorithm introduced in the preceding paper (paper I) [C. Bratschi and H. Huber, J. Chem. Phys. v126, 164104 (2007)] is applied to two recently published CO2 ab initio pair potentials, the Bock-Bich-Vogel and symmetry-adapted perturbation theory site-site potentials. The critical properties of these potentials are calculated for the first time. Critical values and points in the single and two-phase zones are compared with Monte Carlo results to demonstrate the accuracy of the molecular dynamics algorithm, and are compared with experiment to test the accuracy of the potentials. Pressure calculations in the liquid, gas, and supercritical states are carried out and are used to explain potential-related effects and systematic discrepancies. The best ab initio potential yields results in good agreement with experiment.     

Statistical temperature molecular dynamics: application to coarse-grained beta-barrel-forming protein models

Recently the authors proposed a novel sampling algorithm, "statistical temperature molecular dynamics" (STMD) [J. Kim et al., Phys. Rev. Lett. 97, 050601 (2006)], which combines ingredients of multicanonical molecular dynamics and Wang-Landau sampling. Exploiting the relation between the statistical temperature and the density of states, STMD generates a flat energy distribution and efficient sampling with a dynamic update of the statistical temperature, transforming an initial constant estimate to the true statistical temperature T(U), with U being the potential energy. Here, the performance of STMD is examined in the Lennard-Jones fluid with diverse simulation conditions, and in the coarse-grained, off-lattice BLN 46-mer and 69-mer protein models, exhibiting rugged potential energy landscapes with a high degree of frustration. STMD simulations combined with inherent structure (IS) analysis allow an accurate determination of protein thermodynamics down to very low temperatures, overcoming quasiergodicity, and illuminate the transitions occurring in folding in terms of the energy landscape. It is found that a thermodynamic signature of folding is significantly suppressed by accurate sampling, due to an incoherent contribution from low-lying non-native IS in multifunneled landscapes. It is also shown that preferred accessibility to such IS during the collapse transition is intimately related to misfolding or poor foldability.     

Comparison of atomic-level and coarse-grained models for liquid hydrocarbons from molecular dynamics configurational entropy estimates

Molecular liquids can be modeled at different levels of spatial resolution. In atomic-level (AL) models, all (heavy) atoms can be explicitly simulated. In coarse-grained (CG) models, particles (beads) that represent groups of covalently bound atoms are used as elementary units. Ideally, a CG model should reproduce the thermodynamic and structural properties of the corresponding AL model after mapping to the lower-resolution scale. In the present work, two such models are investigated: (i) the classical GROMOS atomic-level model; (ii) a CG model recently proposed by Marrink et al., which maps approximately four non-hydrogen atoms to one bead [J. Phys. Chem. B 2004, 108, 750]. The study is restricted to n-alkanes whose aliphatic fragments are abundantly found in lipids of biological interest. Additionally, cis-9-octadecene is included, as a template chain of the lipid dioleoylphosphatidylcholine (DOPC). The two representations of molecules in the liquid phase are compared in terms of average molecular structures, extent of configurational space sampled, and single-molecule entropies. An approximate method is used to estimate the rotational contributions to the absolute configurational entropy. Good correspondence between the AL and CG representations is found. The loss in configurational entropy due to the reduction in degrees of freedom upon coarse-graining of the model is estimated.     

Dynamics in coarse-grained models for oligomer-grafted silica nanoparticles

Coarse-grained models of poly(ethylene oxide) oligomer-grafted nanoparticles are established by matching their structural distribution functions to atomistic simulation data. Coarse-grained force fields for bulk oligomer chains show excellent transferability with respect to chain lengths and temperature, but structure and dynamics of grafted nanoparticle systems exhibit a strong dependence on the core-core interactions. This leads to poor transferability of the core potential to conditions different from the state point at which the potential was optimized. Remarkably, coarse graining of grafted nanoparticles can either accelerate or slowdown the core motions, depending on the length of the grafted chains. This stands in sharp contrast to linear polymer systems, for which coarse graining always accelerates the dynamics. Diffusivity data suggest that the grafting topology is one cause of slower motions of the cores for short-chain oligomer-grafted nanoparticles; an estimation based on transition-state theory shows the coarse-grained core-core potential also has a slowing-down effect on the nanoparticle organic hybrid materials motions; both effects diminish as grafted chains become longer.     

Polymer chain generation for coarse-grained models using radical-like polymerization

A versatile method is proposed to generate configurations of coarse-grained models for polymer melts. This method, largely inspired by chemical "radical polymerization," is divided in three stages: (i) nucleation of radicals (reacting molecules caching monomers), (ii) growth of chains within a solvent of monomers and (iii) termination: annihilation of radicals and removal of residual monomers. The main interest of this method is that relaxation is performed while chains are generated. Pure mono and polydisperse polymer melts are generated and compared to the configurations generated by the push off method from Auhl et al. [J. Chem. Phys. 119, 12718 (2003)]. A detailed study of the static properties (radius of gyration, mean square internal distance, entanglement length) confirms that the radical-like polymerization technique is suitable to generate equilibrated melts. Moreover, the method is flexible and can be adapted to generate nanostructured polymers, namely, diblock and triblock copolymers.     

Treatment of symmetry in MO calculations. II. Numerical projection operators for molecular orbitals

This article describes the numerical application of projection operators to restore the symmetry of molecular orbitals in self-consistent field (SCF) calculations when the symmetry is lost because of degeneracy or near degeneracy. The application of projection operators is particularly useful in cases of near degeneracies of three or more molecular orbitals, where it is difficult to find an effective algorithm for restoring the symmetry of molecular orbitals by orthonormal transformations.     

Non-Hamiltonian molecular dynamics implementation of the Gibbs ensemble method. I. Algorithm

In this paper (paper I) and the following paper (paper II) [C. Bratschi, H. Huber, and D. J. Searles, J. Chem. Phys. 126, 164105 (2007)], a new molecular dynamics algorithm implementing the Gibbs ensemble will be presented and then on applied to the liquid-vapor coexistence curve for two ab initio CO2 potentials. In paper I, the Gibbs ensemble molecular dynamics algorithm using non-Hamiltonian molecular dynamics techniques is introduced. It is shown that states of the correct probability density function are sampled and the correct exchange probability is generated. The extended system Nose-Hoover formalism is used to generate a constant temperature ensemble with equal pressures in the subsystems, combined with single particle transfers between the subsystems, over several time steps, to get equal chemical potentials.     

Molecular tight-binding method. II. One- and many-electron method for excess electron bands of molecular crystals

An ab initio method for calculating the energies of excess-electron bands in a molecular crystal is developed. These energies represent the electron affinities of a crystal. The present method employs a tight-binding approximation to describe the molecules in a crystal and a set of local functions to describe the excess electron. One- and many-electron formulations of the method are given, the latter takes into account the effect of polarization of all the molecules in the crystal by the excess electron. A scheme for approximate determination of the electronic-correlation corrections to the one-electron bands is developed and applied to calculation of the excess-electron bands in the HCP helium crystal.     

On developing coarse-grained models for biomolecular simulation: a review

So-called coarse-grained models are a popular type of model for accessing long time scales in simulations of biomolecular processes. Such models are coarse-grained with respect to atomic models. But any modelling of processes or substances involves coarse-graining, i.e. the elimination of non-essential degrees of freedom and interactions from a more fine-grained level of modelling. The basic ingredients of developing coarse-grained models based on the properties of fine-grained models are reviewed, together with the conditions that must be satisfied in order to preserve the correct physical mechanisms in the coarse-graining process. This overview should help the reader to determine how realistic a coarse-grained model of a biomolecular system is, i.e. whether it reflects the underlying physical mechanisms or merely provides a set of pretty pictures of the process or substances of interest.     

Comparison of thermodynamic properties of coarse-grained and atomic-level simulation models.

Thermodynamic data are often used to calibrate or test amomic-level (AL) force fields for molecular dynamics (MD) simulations. In contrast, the majority of coarse-grained (CG) force fields do not rely extensively on thermodynamic quantities. Recently, a CG force field for lipids, hydrocarbons, ions, and water, in which approximately four non-hydrogen atoms are mapped onto one interaction site, has been proposed and applied to study various aspects of lipid systems. To date, no extensive investigation of its capability to describe salvation thermodynamics has been undertaken. In the present study, a detailed picture of vaporization, solvation, and phase-partitioning thermodynamics for liquid hydrocarbons and water was obtained at CG and AL resolutions, in order to compare the two types or models and evaluate their ability to describe thermodynamic properties in the temperature range between 263 and 343 K. Both CG and AL models capture the experimental dependence of the thermodynamic properties on the temperature, albeit a systematically weaker dependence is found for the CG model. Moreover, deviations are found for solvation thermodynamics and for the corresponding enthalpy-entropy compensation for the CG model. Particularly water/oil repulsion seems to be overestimated. However, the results suggest that the thermodynamic properties considered should be reproducible by a CG model provided it is reparametrized on the basis of these liquid-phase properties.     

IBIsCO: a molecular dynamics simulation package for coarse-grained simulation

IBIsCO is a parallel molecular dynamics simulation package developed specially for coarse-grained simulations with numerical potentials derived by the iterative Boltzmann inversion (IBI) method (Reith et al., J Comput Chem 2003, 24, 1624). In addition to common features of molecular dynamics programs, the techniques of dissipative particle dynamics (Groot and Warren, J Chem Phys 1997, 107, 4423) and Lowe-Andersen dynamics (Lowe, Europhys Lett 1999, 47, 145) are implemented, which can be used both as thermostats and as sources of friction to compensate the loss of degrees of freedom by coarse-graining. The reverse nonequilibrium molecular dynamics simulation method (Müller-Plathe, Phys Rev E 1999, 59, 4894) for the calculation of viscosities is also implemented. Details of the algorithms, functionalities, implementation, user interfaces, and file formats are described. The code is parallelized using PE_MPI on PowerPC architecture. The execution time scales satisfactorily with the number of processors.     

Symmetry in the Heitler-London method. II. The determination of the allowed molecular multiplets

In the framework of the Heitler-London method, a method for determining the allowed molecular multiplets is proposed. The method is based on the connection of the total molecular spin with the permutation symmetry of the coordinate wave function and on the isomorphism of the molecular point group with a certain subgroup of the electron permutation group. The method does not depend on the approximation in which the molecular ions are considered.