Measuring friction from simulations of folded graphene sheets.

We ran molecular dynamics simulations of folded graphene sheets and present a procedure to measure the sliding friction in these systems based on the rate of decay of a damped harmonic oscillator. This procedure allowed us to study the effect that the size, geometry, and the temperature of the graphene sheet had on the ability to propagate the initial fold and the rate at which it settles to a final "fully folded" equilibrium state. We offer simple rationalizations for the relationships between the initial geometries of our simulations and the friction values that emerge.

A fully quantum-mechanical treatment for kaolinite.

Neural network potentials for kaolinite minerals have been fitted to data extracted from density functional theory calculations that were performed using the revPBE + D3 and revPBE + vdW functionals. These potentials have then been used to calculate the static and dynamic properties of the mineral. We show that revPBE + vdW is better at reproducing the static properties. However, revPBE + D3 does a better job of reproducing the experimental IR spectrum. We also consider what happens to these properties when a fully quantum treatment of the nuclei is employed. We find that nuclear quantum effects (NQEs) do not make a substantial difference to the static properties. However, when NQEs are included, the dynamic properties of the material change substantially.

Interplay of multiple clusters and initial interface positioning for forward flux sampling simulations of crystal nucleation.

Forward flux sampling (FFS) is a path sampling technique widely used in computer simulations of crystal nucleation from the melt. In such studies, the order parameter underpinning the progress of the FFS algorithm is often the size of the largest crystalline nucleus. In this work, we investigate the effects of two computational aspects of FFS simulations, using the prototypical Lennard-Jones liquid as our computational test bed. First, we quantify the impact of the positioning of the liquid basin and first interface in the space of the order parameter. In particular, we demonstrate that these choices are key to ensuring the consistency of the FFS results. Second, we focus on the frequently encountered scenario where the population of crystalline nuclei is such that there are multiple clusters of size comparable to the largest one. We demonstrate the contribution of clusters other than the largest cluster to the initial flux; however, we also show that they can be safely ignored for the purposes of converging a full FFS calculation. We also investigate the impact of different clusters merging, a process that appears to be facilitated by substantial spatial correlations-at least at the supercooling considered here. Importantly, all of our results have been obtained as a function of system size, thus contributing to the ongoing discussion on the impact of finite size effects on simulations of crystal nucleation. Overall, this work either provides or justifies several practical guidelines for performing FFS simulations that can also be applied to more complex and/or computationally expensive models.

Can molecular dynamics be used to simulate biomolecular recognition?

There are many problems in biochemistry that are difficult to study experimentally. Simulation methods are appealing due to direct availability of atomic coordinates as a function of time. However, direct molecular simulations are challenged by the size of systems and the time scales needed to describe relevant motions. In theory, enhanced sampling algorithms can help to overcome some of the limitations of molecular simulations. Here, we discuss a problem in biochemistry that offers a significant challenge for enhanced sampling methods and that could, therefore, serve as a benchmark for comparing approaches that use machine learning to find suitable collective variables. In particular, we study the transitions LacI undergoes upon moving between being non-specifically and specifically bound to DNA. Many degrees of freedom change during this transition and that the transition does not occur reversibly in simulations if only a subset of these degrees of freedom are biased. We also explain why this problem is so important to biologists and the transformative impact that a simulation of it would have on the understanding of DNA regulation.

Classical nucleation theory predicts the shape of the nucleus in homogeneous solidification.

Macroscopic models of nucleation provide powerful tools for understanding activated phase transition processes. These models do not provide atomistic insights and can thus sometimes lack material-specific descriptions. Here, we provide a comprehensive framework for constructing a continuum picture from an atomistic simulation of homogeneous nucleation. We use this framework to determine the equilibrium shape of the solid nucleus that forms inside bulk liquid for a Lennard-Jones potential. From this shape, we then extract the anisotropy of the solid-liquid interfacial free energy, by performing a reverse Wulff construction in the space of spherical harmonic expansions. We find that the shape of the nucleus is nearly spherical and that its anisotropy can be perfectly described using classical models.

Building maps in collective variable space.

Enhanced sampling techniques such as umbrella sampling and metadynamics are now routinely used to provide information on how the thermodynamic potential, or free energy, depends on a small number of collective variables (CVs). The free energy surfaces that one extracts by using these techniques provide a simplified or coarse-grained representation of the configurational ensemble. In this work, we discuss how auxiliary variables can be mapped in CV space. We show that maps of auxiliary variables allow one to analyze both the physics of the molecular system under investigation and the quality of the reduced representation of the system that is encoded in a set of CVs. We apply this approach to analyze the degeneracy of CVs and to compute entropy and enthalpy surfaces in CV space both for conformational transitions in alanine dipeptide and for phase transitions in carbon dioxide molecular crystals under pressure.

The Gibbs free energy of homogeneous nucleation: From atomistic nuclei to the planar limit.

In this paper we discuss how the information contained in atomistic simulations of homogeneous nucleation should be used when fitting the parameters in macroscopic nucleation models. We show how the number of solid and liquid atoms in such simulations can be determined unambiguously by using a Gibbs dividing surface and how the free energy as a function of the number of solid atoms in the nucleus can thus be extracted. We then show that the parameters (the chemical potential, the interfacial free energy, and a Tolman correction) of a model based on classical nucleation theory can be fitted using the information contained in these free-energy profiles but that the parameters in such models are highly correlated. This correlation is unfortunate as it ensures that small errors in the computed free energy surface can give rise to large errors in the extrapolated properties of the fitted model. To resolve this problem we thus propose a method for fitting macroscopic nucleation models that uses simulations of planar interfaces and simulations of three-dimensional nuclei in tandem. We show that when the chemical potentials and the interface energy are pinned to their planar-interface values, more precise estimates for the Tolman length are obtained. Extrapolating the free energy profile obtained from small simulation boxes to larger nuclei is thus more reliable.

Porous Organic Cages for Sulfur Hexafluoride Separation.

A series of porous organic cages is examined for the selective adsorption of sulfur hexafluoride (SF6) over nitrogen. Despite lacking any metal sites, a porous cage, CC3, shows the highest SF6/N2 selectivity reported for any material at ambient temperature and pressure, which translates to real separations in a gas breakthrough column. The SF6 uptake of these materials is considerably higher than would be expected from the static pore structures. The location of SF6 within these materials is elucidated by X-ray crystallography, and it is shown that cooperative diffusion and structural rearrangements in these molecular crystals can rationalize their superior SF6/N2 selectivity.

Ice formation on kaolinite: Insights from molecular dynamics simulations.

The formation of ice affects many aspects of our everyday life as well as important technologies such as cryotherapy and cryopreservation. Foreign substances almost always aid water freezing through heterogeneous ice nucleation, but the molecular details of this process remain largely unknown. In fact, insight into the microscopic mechanism of ice formation on different substrates is difficult to obtain even if state-of-the-art experimental techniques are used. At the same time, atomistic simulations of heterogeneous ice nucleation frequently face extraordinary challenges due to the complexity of the water-substrate interaction and the long time scales that characterize nucleation events. Here, we have investigated several aspects of molecular dynamics simulations of heterogeneous ice nucleation considering as a prototypical ice nucleating material the clay mineral kaolinite, which is of relevance in atmospheric science. We show via seeded molecular dynamics simulations that ice nucleation on the hydroxylated (001) face of kaolinite proceeds exclusively via the formation of the hexagonal ice polytype. The critical nucleus size is two times smaller than that obtained for homogeneous nucleation at the same supercooling. Previous findings suggested that the flexibility of the kaolinite surface can alter the time scale for ice nucleation within molecular dynamics simulations. However, we here demonstrate that equally flexible (or non flexible) kaolinite surfaces can lead to very different outcomes in terms of ice formation, according to whether or not the surface relaxation of the clay is taken into account. We show that very small structural changes upon relaxation dramatically alter the ability of kaolinite to provide a template for the formation of a hexagonal overlayer of water molecules at the water-kaolinite interface, and that this relaxation therefore determines the nucleation ability of this mineral.

Locating binding poses in protein-ligand systems using reconnaissance metadynamics.

A molecular dynamics-based protocol is proposed for finding and scoring protein-ligand binding poses. This protocol uses the recently developed reconnaissance metadynamics method, which employs a self-learning algorithm to construct a bias that pushes the system away from the kinetic traps where it would otherwise remain. The exploration of phase space with this algorithm is shown to be roughly six to eight times faster than unbiased molecular dynamics and is only limited by the time taken to diffuse about the surface of the protein. We apply this method to the well-studied trypsin-benzamidine system and show that we are able to refind all the poses obtained from a reference EADock blind docking calculation. These poses can be scored based on the length of time the system remains trapped in the pose. Alternatively, one can perform dimensionality reduction on the output trajectory and obtain a map of phase space that can be used in more expensive free-energy calculations.

Using sketch-map coordinates to analyze and bias molecular dynamics simulations.

When examining complex problems, such as the folding of proteins, coarse grained descriptions of the system drive our investigation and help us to rationalize the results. Oftentimes collective variables (CVs), derived through some chemical intuition about the process of interest, serve this purpose. Because finding these CVs is the most difficult part of any investigation, we recently developed a dimensionality reduction algorithm, sketch-map, that can be used to build a low-dimensional map of a phase space of high-dimensionality. In this paper we discuss how these machine-generated CVs can be used to accelerate the exploration of phase space and to reconstruct free-energy landscapes. To do so, we develop a formalism in which high-dimensional configurations are no longer represented by low-dimensional position vectors. Instead, for each configuration we calculate a probability distribution, which has a domain that encompasses the entirety of the low-dimensional space. To construct a biasing potential, we exploit an analogy with metadynamics and use the trajectory to adaptively construct a repulsive, history-dependent bias from the distributions that correspond to the previously visited configurations. This potential forces the system to explore more of phase space by making it desirable to adopt configurations whose distributions do not overlap with the bias. We apply this algorithm to a small model protein and succeed in reproducing the free-energy surface that we obtain from a parallel tempering calculation.

From the Cover: Simplifying the representation of complex free-energy landscapes using sketch-map.

A new scheme, sketch-map, for obtaining a low-dimensional representation of the region of phase space explored during an enhanced dynamics simulation is proposed. We show evidence, from an examination of the distribution of pairwise distances between frames, that some features of the free-energy surface are inherently high-dimensional. This makes dimensionality reduction problematic because the data does not satisfy the assumptions made in conventional manifold learning algorithms We therefore propose that when dimensionality reduction is performed on trajectory data one should think of the resultant embedding as a quickly sketched set of directions rather than a road map. In other words, the embedding tells one about the connectivity between states but does not provide the vectors that correspond to the slow degrees of freedom. This realization informs the development of sketch-map, which endeavors to reproduce the proximity information from the high-dimensionality description in a space of lower dimensionality even when a faithful embedding is not possible.

A self-learning algorithm for biased molecular dynamics.

A new self-learning algorithm for accelerated dynamics, reconnaissance metadynamics, is proposed that is able to work with a very large number of collective coordinates. Acceleration of the dynamics is achieved by constructing a bias potential in terms of a patchwork of one-dimensional, locally valid collective coordinates. These collective coordinates are obtained from trajectory analyses so that they adapt to any new features encountered during the simulation. We show how this methodology can be used to enhance sampling in real chemical systems citing examples both from the physics of clusters and from the biological sciences.

Exploring the free energy surfaces of clusters using reconnaissance metadynamics.

A new approach is proposed for exploring the low-energy structures of small to medium-sized aggregates of atoms and molecules. This approach uses the recently proposed reconnaissance metadynamics method [G. A. Tribello, M. Ceriotti, and M. Parrinello. Proc. Natl. Acad. Sci. U.S.A. 107(41), 17509 (2010)] in tandem with collective variables that describe the average structure of the coordination sphere around the atoms/molecules. We demonstrate this method on both Lennard-Jones and water clusters and show how it is able to quickly find the global minimum in the potential energy surface, while exploring the finite temperature free energy surface.

Asprich peptides are occluded in calcite and permanently disorder biomineral crystals.

Macromolecules are a minority but important component of the minerals formed by living organisms, or biominerals. The role these macromolecules play at the early stages of biomineral formation, as well as their long-term and long-range effects on the mature biomineral, is poorly understood. A 42-amino acid peptide, asp2, was derived from the Asprich family of proteins. In this study we present X-ray absorption near-edge structure spectroscopy and X-ray photoelectron emission microscopy data from the asp2 peptide, the calcite (CaCO(3)) crystals, and the peptide + crystal composites. The results clearly show that asp2 is occluded in fully formed biomineral crystals and slightly but permanently disorders the crystal structure at short- and long-range distances.

Isomorphism between ice and silica.

Both ice and silica crystallize into solid-state structures composed of tetrahedral building units that are joined together to form an infinite four-connected net. Mathematical considerations suggest that there is a vast number of such nets and thus potential crystal structures. It is therefore perhaps surprising to discover that, despite the differences in the nature of interatomic interactions in these materials, a fair number of commonly observed ice and silica phases are based on common nets. Here we use computer simulation to investigate the origin of this symmetry between the structures formed for ice and silica and to attempt to understand why it is not complete. We start from a comparison of the dense phases and then move to the relationship between the different open (zeolitic and clathratic) structures formed for both materials. We show that there is a remarkably strong correlation between the energetics of isomorphic silica and water ice structures and that this correlation arises because of the strong link between the total energy of a material and its local geometric features. Finally, we discuss a number of as yet unsynthesized low-energy structures which include a phase of ice based on quartz, a silica based on the structure of ice VI, and an ice clathrate that is isomorphic to the silicate structure nonasil.

Combining Machine Learning and Enhanced Sampling Techniques for Efficient and Accurate Calculation of Absolute Binding Free Energies.

Calculating absolute binding free energies is challenging and important. In this paper, we test some recently developed metadynamics-based methods and develop a new combination with a Hamiltonian replica-exchange approach. The methods were tested on 18 chemically diverse ligands with a wide range of different binding affinities to a complex target; namely, human soluble epoxide hydrolase. The results suggest that metadynamics with a funnel-shaped restraint can be used to calculate, in a computationally affordable and relatively accurate way, the absolute binding free energy for small fragments. When used in combination with an optimal pathlike variable obtained using machine learning or with the Hamiltonian replica-exchange algorithm SWISH, this method can achieve reasonably accurate results for increasingly complex ligands, with a good balance of computational cost and speed. An additional benefit of using the combination of metadynamics and SWISH is that it also provides useful information about the role of water in the binding mechanism.

Binding of calcium and carbonate to polyacrylates.

Polyacrylate molecules can be used to slow the growth of calcium carbonate. However, little is known about the mechanism by which the molecules impede the growth rate. A recent computational study (Bulo et al. Macromolecules 2007, 40, 3437) used metadynamics to investigate the binding of calcium to polyacrylate chains and has thrown some light on the coiling and precipitation of these polymers. We extend these simulations to examine the binding of calcium and carbonate to polyacrylate chains. We show that calcium complexed with both carbonate and polyacrylate is a very stable species. The free energies of calcium-carbonate-polyacrylate complexes, with different polymer configurations, are calculated, and differences in the free energy of the binding of carbonate are shown to be due to differences in the amount of steric hindrance about the calcium, which prevents the approach of the carbonate ion.

A molecular dynamics study of the early stages of calcium carbonate growth.

The precipitation of calcium carbonate in water has been examined using a combination of molecular dynamics and umbrella sampling. During 20 ns molecular dynamics trajectories at elevated calcium carbonate concentrations, amorphous particles are observed to form and appear to be composed of misaligned domains of vaterite and aragonite. The addition of further calcium ions to these clusters is found to be energetically favorable and virtually barrierless. By contrast, there is a large barrier to the addition of calcium to small calcite crystals. Thus, even though calcite nanocrystals are stable in solution, at high supersaturations, particles of amorphous material form because this material grows much faster than ordered calcite nanocrystals.

A theoretical examination of known and hypothetical clathrate hydrate materials.

The recent synthesis of a new hydrogen binary hydrate with the sH structure has highlighted the potential storage capabilities of water clathrates [T. A. Strobel, C. A. Koh, and E. D. Sloan, J. Phys. Chem. B 112, 1885 (2008) and A. R. C. Duarte, A. Shariati, L. J. Rovetto, and C. J. Peters, J. Phys. Chem. B 112, 1888 (2008)]. In this work, the absorption of hydrogen and the promoters used in the experimental work are considered using a simplified model for the host-guest interaction, which allows one to understand the stabilizing effects of multiple help molecules. Two further hypothetical clathrates, which are isostructural with known zeolite structures, are also investigated. It is shown that the energy gained by absorbing adamantane into these two frameworks is far greater than that gained upon absorption of adamantane into the sH structure. Hence, a clathrate with the same topology as the DDR (Sigma 1) zeolite may be synthesizable with adamantane and hydrogen as guest molecules as, in the conditions explored here, this phase appears to be more stable than the sH structure.