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.

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.

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.

Global Free-Energy Landscapes as a Smoothly Joined Collection of Local Maps

Enhanced sampling techniques have become an essential tool in computational chemistry and physics, where they are applied to sample activated processes that occur on a time scale that is inaccessible to conventional simulations. Despite their popularity, it is well known that they have constraints that hinder their application to complex problems. The core issue lies in the need to describe the system using a small number of collective variables (CVs). Any slow degree of freedom that is not properly described by the chosen CVs will hinder sampling efficiency. However, the exploration of configuration space is also hampered by including variables that are not relevant for the activated process under study. This paper presents the Adaptive Topography of Landscape for Accelerated Sampling (ATLAS), a new biasing method capable of working with many CVs. The root idea of ATLAS is to apply a divide-and-conquer strategy, where the high-dimensional CVs space is divided into basins, each of which is described by an automatically determined, low-dimensional set of variables. A well-tempered metadynamics-like bias is constructed as a function of these local variables. Indicator functions associated with the basins switch on and off the local biases so that the sampling is performed on a collection of low-dimensional CV spaces that are smoothly combined to generate an effectively high-dimensional bias. The unbiased Boltzmann distribution is recovered through reweighing, making the evaluation of conformational and thermodynamic properties straightforward. The decomposition of the free-energy landscape in local basins can be updated iteratively as the simulation discovers new (meta)stable states.

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.

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.

Using Dimensionality Reduction to Analyze Protein Trajectories

In recent years the analysis of molecular dynamics trajectories using dimensionality reduction algorithms has become commonplace. These algorithms seek to find a low-dimensional representation of a trajectory that is, according to a well-defined criterion, optimal. A number of different strategies for generating projections of trajectories have been proposed but little has been done to systematically compare how these various approaches fare when it comes to analysing trajectories for biomolecules in explicit solvent. In the following paper, we have thus analyzed a molecular dynamics trajectory of the C-terminal fragment of the immunoglobulin binding domain B1 of protein G of Streptococcus modeled in explicit solvent using a range of different dimensionality reduction algorithms. We have then tried to systematically compare the projections generated using each of these algorithms by using a clustering algorithm to find the positions and extents of the basins in the high-dimensional energy landscape. We find that no algorithm outshines all the other in terms of the quality of the projection it generates. Instead, all the algorithms do a reasonable job when it comes to building a projection that separates some of the configurations that lie in different basins. Having said that, however, all the algorithms struggle to project the basins because they all have a large intrinsic dimensionality.

Insights into mechanochemical reactions at the molecular level: simulated indentations of aspirin and meloxicam crystals

Although solvent-free mechanochemical synthesis continues to gain ever greater importance, the molecular scale processes that occur during such reactions remain largely uncharacterised. Here, we apply computational modelling to indentations between particles of crystals of aspirin and meloxicam under a variety of conditions to mimic the early stages of their mechanochemical cocrystallisation reaction. The study also extends to the effects of the presence of small amounts of solvent. It is found that, despite the solid crystalline nature of the reactants and the presence of little or no solvent, mixing occurs readily at the molecular level even during relatively low-energy collisions. When indented crystals are subsequently drawn apart, a connective neck formed by a mixture of the reactant molecules is observed, suggesting plastic-like behaviour of the reacting materials. Overall the work reveals some striking new insights including (i) relatively facile mixing of crystals under solvent-free conditions, (ii) no appreciable local temperature increases, (iii) localised amorphisation at the contact region and neck of the reacting crystals, and (iv) small amounts of solvent have relatively little effect during this early stage of the reaction, suggesting that their accelerating effect on the reaction may be exerted at later stages.

Solvation Effects on Dissociative Electron Attachment to Thymine

Ionizing radiation can excite the cellular medium to produce secondary electrons that can subsequently cause damage to DNA. The damage is believed to occur via dissociative electron attachment (DEA). In DEA, the electron is captured by a molecule in a resonant antibonding state and a transient negative ion is formed. If this ion survives against electron autodetachment, then bonds within the molecule may dissociate as energy is transferred from the electronic degrees of freedom into vibrational modes of the molecule. We present a model for studying the effect that transferring kinetic energy into the vibrational modes of a molecule in this way has on a DNA nucleobase. We show that when the base is in an aqueous environment, dissociation is affected by interactions with the surrounding water molecules. In particular, hydrogen bonding between the nucleobase and the solvent can suppress the dissociative channel.

Analyzing and Biasing Simulations with PLUMED

This chapter discusses how the PLUMED plugin for molecular dynamics can be used to analyze and bias molecular dynamics trajectories. The chapter begins by introducing the notion of a collective variable and by then explaining how the free energy can be computed as a function of one or more collective variables. A number of practical issues mostly around periodic boundary conditions that arise when these types of calculations are performed using PLUMED are then discussed. Later parts of the chapter discuss how PLUMED can be used to perform enhanced sampling simulations that introduce simulation biases or multiple replicas of the system and Monte Carlo exchanges between these replicas. This section is then followed by a discussion on how free-energy surfaces and associated error bars can be extracted from such simulations by using weighted histogram and block averaging techniques.

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.

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.

Extracting the interfacial free energy and anisotropy from a smooth fluctuating dividing surface

Interfaces between different materials and phases play a crucial role in many physical and chemical phenomena. When performing simulations of matter at the atomic scale, however, it is often not trivial to characterize these interfaces, particularly when they are rough or diffuse. Here we discuss a generalization of a construction, due to Willard and Chandler, that allows one to obtain a smooth dividing surface that follows the irregular, ever changing shape of these fluctuating interfaces. We show how this construction can be used to study the surface that separates a solid material from its melt and how analyses of the Fourier modes for the capillary fluctuations of this instantaneous dividing surface can be performed. This particular analysis is useful as one can compute the specific free energy excess of the interface, and its dependence on orientation relative to the bulk phases, from the average amplitude of the Fourier modes. We therefore discuss the efficiency of this approach, both in terms of system size and statistical sampling.

Interactions between low energy electrons and DNA: a perspective from first-principles simulations

DNA damage caused by irradiation has been studied for many decades. Such studies allow us to better assess the dangers posed by radiation, and to increase the efficiency of the radiotherapies that are used to combat cancer. A full description of the irradiation process involves multiple size and time scales. It starts with the interaction of radiation-either photons or swift ions-and the biological medium, which causes electronic excitation and ionisation. The two main products of ionising radiation are thus electrons and radicals. Both of these species can cause damage to biological molecules, in particular DNA. In the long run, this molecular level damage can prevent cells from replicating and can hence lead to cell death. For a long time it was assumed that the main actors in the damage process were the radicals. However, experiments in a seminal paper by the group of Leon Sanche in 2000 showed that low-energy electrons (LEE), such as those generated when ionising biological targets, can also cause bond breaks in biomolecules, and strand breaks in plasmid DNA in particular (Boudaiffa et al 2000 Science 287 1658-60). These results prompted a significant amount of experimental and theoretical work aimed at elucidating the role played by LEE in DNA damage. In this Topical Review we provide a general overview of the problem. We discuss experimental findings and theoretical results hand in hand with the aim of describing the physics and chemistry that occurs during the process of radiation damage, from the initial stages of electronic excitation, through the inelastic propagation of electrons in the medium, the interaction of electrons with DNA, and the chemical end-point effects on DNA. A very important aspect of this discussion is the consideration of a realistic, physiological environment. The role played by the aqueous solution and the amino acids from the histones in chromatin must be considered. Moreover, thermal fluctuations must be incorporated when studying these phenomena. Hence, a special place in this Topical Review is occupied by our recent first-principles molecular dynamics simulations that address the issue of how the environment favours or prevents LEEs from causing damage to DNA. We finish by summarising the conclusions achieved so far, and by suggesting a number of possible directions for further study.

Predicting solvent effects on the structure of porous organic molecules

A computational approach for the prediction of the open, metastable, conformations of porous organic molecules in the presence of solvent is developed.

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.

Understanding the Interaction between Low-Energy Electrons and DNA Nucleotides in Aqueous Solution

Reactions that can damage DNA have been simulated using a combination of molecular dynamics and density functional theory. In particular, the damage caused by the attachment of a low energy electron to the nucleobase. Simulations of anionic single nucleotides of DNA in an aqueous environment that was modeled explicitly have been performed. This has allowed us to examine the role played by the water molecules that surround the DNA in radiation damage mechanisms. Our simulations show that hydrogen bonding and protonation of the nucleotide by the water can have a significant effect on the barriers to strand breaking reactions. Furthermore, these effects are not the same for all four of the bases.

Surface energy and surface proton order of the ice Ih basal and prism surfaces

Density-functional theory (DFT) is used to examine the basal and prism surfaces of ice Ih. Similar surface energies are obtained for the two surfaces; however, in each case a strong dependence of the surface energy on surface proton order is identified. This dependence, which can be as much as 50% of the absolute surface energy, is significantly larger than the bulk dependence (<1%) on proton order, suggesting that the thermodynamic ground state of the ice surface will remain proton ordered well above the bulk order-disorder temperature of about 72 K. On the basal surface this suggestion is supported by Monte Carlo simulations with an empirical potential and solution of a 2D Ising model with nearest neighbor interactions taken from DFT. Order parameters that define the surface energy of each surface in terms of nearest neighbor interactions between dangling OH bonds (those which point out of the surface into vacuum) have been identified and are discussed. Overall, these results suggest that proton order-disorder effects have a profound impact on the stability of ice surfaces and will most likely have an effect on ice surface reactivity as well as ice crystal growth and morphology.

Surface energy and surface proton order of ice Ih

Ice Ih is comprised of orientationally disordered water molecules giving rise to positional disorder of the hydrogen atoms in the hydrogen bonded network of the lattice. Here we arrive at a first principles determination of the surface energy of ice Ih and suggest that the surface of ice is significantly more proton ordered than the bulk. We predict that the proton order-disorder transition, which occurs in the bulk at approximately 72 K, will not occur at the surface at any temperature below surface melting. An order parameter which defines the surface energy of ice Ih surfaces is also identified.

A Blind Structure Prediction of Ice XIV

Evidence points to several hydrogen ordered ice phases of ice XII that exhibit topologically distinct hydrogen networks from one another and which are formed under kinetically controlled conditions.