Calculation of free energy landscapes: A histogram reweighted metadynamics approach

We present an efficient method for the calculation of free energy landscapes. Our approach involves a history‐dependent bias potential, which is evaluated on a grid. The corresponding free energy landscape is constructed via a histogram reweighting procedure a posteriori. Because of the presence of the bias potential, it can be also used to accelerate rare events. In addition, the calculated free energy landscape is not restricted to the actual choice of collective variables and can in principle be extended to auxiliary variables of interest without further numerical effort. The applicability is shown for several examples. We present numerical results for the alanine dipeptide and the Met‐Enkephalin in explicit solution to illustrate our approach. Furthermore, we derive an empirical formula that allows the prediction of the computational cost for the ordinary metadynamics variant in comparison with our approach, which is validated by a dimensionless representation. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Calculation of solid-liquid interfacial free energy: a classical nucleation theory based approach

We present a simple approach to calculate the solid-liquid interfacial free energy. This new method is based on the classical nucleation theory. Using the molecular dynamics simulation, we employ spherical crystal nuclei embedded in the supercooled liquids to create an ideal model of a homogeneous nucleation. The interfacial free energy is extracted by fitting the relation between the critical nucleus size and the reciprocal of the critical undercooling temperature. The orientationally averaged interfacial free energy is found to be 0.302+/-0.002 (in standard LJ unit). The temperature dependence of the interfacial free energy is also obtained in this work. We find that the interfacial free energy increases slightly with increasing temperature. The positive temperature coefficient of the interfacial free energy is in qualitative agreement with Spaepen's analysis [Solid State Phys. 47, FS181 (1994)] and Turnbull's empirical estimation [J. Appl. Phys. 21, 1022 (1950)].     

Hit and miss of classical nucleation theory as revealed by a molecular simulation study of crystal nucleation in supercooled sulfur hexafluoride

Classical nucleation theory pictures the homogeneous nucleation of a crystal as the formation of a spherical crystalline embryo, possessing the properties of the macroscopic crystal, inside a parent supercooled liquid. In this work we study crystal nucleation in moderately supercooled sulfur hexafluoride by umbrella sampling simulations. The nucleation free energy evolves from 5.2kBT at T=170 K to 39.1kBT at T=195 K. The corresponding critical nucleus size ranges from 40 molecules at T=170 K to 266 molecules at T=195 K. Both nucleation free energy and critical nucleus size are shown to evolve with temperature according to the equations derived from the classical nucleation theory. Inspecting the obtained nuclei we show, however, that they present quite anisotropic shapes in opposition to the spherical assumption of the theory. Moreover, even though the critical nuclei possess the structure of the stable bcc plastic phase, the only mechanically stable crystal phase for SF6 in the temperature range investigated, they are shown to be less ordered than the corresponding macroscopic crystal. Their crystalline order is nevertheless shown to increase regularly with their size. This is confirmed by a study of a nucleus growth from a critical size to a size of the order of 10(4) molecules. Similarly to the fact that it does not affect the temperature dependence of the nucleation free energy and of the critical nucleus size, the ordering of the nucleus with size does not affect the growth rate of the nucleus.     

New method for determining size of critical nucleus of fibril formation of polypeptide chains

A new method for determining the size of critical nucleus of fibril formation of polypeptide chains is proposed. Based on the hypothesis that the fibril grows by addition of a nascent peptide to the preformed template, the nucleus size N(c) is defined as the number of forming template peptides above which the time to add a new monomer becomes independent of the template size. Using lattice models one can show that our method and the standard method which is based on calculation of the free energy, provide the same result for N(c).     

Crystal nucleation and the solid-liquid interfacial free energy

We present the results of molecular dynamics simulation of crystal nucleation in a supercooled Lennard-Jones liquid. Temperature and baric dependences of the nucleation rate, the Zeldovich factor, nucleus size diffusion coefficient, the radius, and the pressure in a critical crystal nucleus are defined in computer simulation. The data obtained have been used in the framework of classical nucleation theory to calculate the effective surface energy of crystal nuclei γ(e). It is shown that the value of γ(e) at T = const exceeds the value of the interfacial free energy at a flat crystal-liquid interface γ(∞) and γ(e) < γ(∞) at p = const.     

Atomistic simulation of the homogeneous nucleation and of the growth of N2 crystallites

We report on a computer simulation study of the early stages of the crystallization of molecular nitrogen. First, we study how homogeneous nucleation takes place in supercooled liquid N(2) for a moderate degree of supercooling. Using the umbrella sampling technique, we determine the free energy barrier of formation for a critical nucleus of N(2). We show that, in accord with Ostwald's rule of stages, the structure of the critical nucleus is predominantly that of a metastable polymorph (alpha-N(2) for the state point investigated). We then monitor the evolution of several critical nuclei through a series of unbiased molecular dynamics trajectories. The growth of N(2) crystallites is accompanied by a structural evolution toward the stable polymorph beta-N(2). The microscopic mechanism underlying this evolution qualitatively differs from that reported previously. We do not observe any dissolution or reorganization of the alpha-like core of the nucleus. On the contrary, we show that alpha-like and beta-like blocks coexist in postcritical nuclei. We relate the structural evolution to a greater adsorption rate of beta-like molecules on the surface and show that this transition actually starts well within the precritical regime. We also carefully investigate the effect of the system size on the height of the free energy barrier of nucleation and on the structure and size of the critical nucleus.     

Multidimensional view of amyloid fibril nucleation in atomistic detail

Starting from a disordered aggregate, we have simulated the formation of ordered amyloid-like beta structures in a system formed by 18 polyvaline chains in explicit solvent, employing molecular dynamics accelerated by bias-exchange metadynamics. We exploited 8 different collective variables to compute the free energy of hundreds of putative aggregate structures, with variable content of parallel and antiparallel β-sheets and different packing among the sheets. This allowed characterizing in detail a possible nucleation pathway for the formation of amyloid fibrils: first the system forms a relatively large ordered nucleus of antiparallel β-sheets, and then a few parallel sheets start appearing. The relevant nucleation process culminates at this point: when a sufficient number of parallel sheets is formed, the free energy starts to decrease toward a new minimum in which this structure is predominant. The complex nucleation pathway we found cannot be described within classical nucleation theory, namely employing a unique simple reaction coordinate like the total content of β-sheets.     

Harmonic Fourier beads method for studying rare events on rugged energy surfaces

We present a robust, distributable method for computing minimum free energy paths of large molecular systems with rugged energy landscapes. The method, which we call harmonic Fourier beads (HFB), exploits the Fourier representation of a path in an appropriate coordinate space and proceeds iteratively by evolving a discrete set of harmonically restrained path points-beads-to generate positions for the next path. The HFB method does not require explicit knowledge of the free energy to locate the path. To compute the free energy profile along the final path we employ an umbrella sampling method in two generalized dimensions. The proposed HFB method is anticipated to aid the study of rare events in biomolecular systems. Its utility is demonstrated with an application to conformational isomerization of the alanine dipeptide in gas phase.     

Formation mechanism of critical nucleus during nucleation process of liquid metal sodium

To deeply understand the formation mechanism of a critical nucleus during the nucleation process of liquid metal sodium, a system consisting of 10 000 Na atoms has been simulated by using molecular dynamics method. The evolutions of nuclei are traced directly, adopting the cluster-type index method. It is found that the energies of clusters and their geometrical constraints interplay to form the favorable microstructures during the nucleation process. The nucleus can be formed through many different pathways, and the critical size of the nucleus would be different for each pathway. It is also found that the critical nucleus is nonspherical and may include some metastable structures. Furthermore, the size of the cluster and its internal structure both play a crucial role in determining whether it is a critical nucleus, and this is in agreement with the simulations by computing the free energy of the Lennard-Jones system [D. Moroni, P. R. ten Wolde, and P. G. Bolhuis, Phys. Rev. Lett. 94, 235703 (2005)].     

GENESIS 1.1: A hybrid‐parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms

GENeralized‐Ensemble SImulation System (GENESIS) is a software package for molecular dynamics (MD) simulation of biological systems. It is designed to extend limitations in system size and accessible time scale by adopting highly parallelized schemes and enhanced conformational sampling algorithms. In this new version, GENESIS 1.1, new functions and advanced algorithms have been added. The all‐atom and coarse‐grained potential energy functions used in AMBER and GROMACS packages now become available in addition to CHARMM energy functions. The performance of MD simulations has been greatly improved by further optimization, multiple time‐step integration, and hybrid (CPU + GPU) computing. The string method and replica‐exchange umbrella sampling with flexible collective variable choice are used for finding the minimum free‐energy pathway and obtaining free‐energy profiles for conformational changes of a macromolecule. These new features increase the usefulness and power of GENESIS for modeling and simulation in biological research. © 2017 Wiley Periodicals, Inc.     

String method in collective variables: minimum free energy paths and isocommittor surfaces

A computational technique is proposed which combines the string method with a sampling technique to determine minimum free energy paths. The technique only requires to compute the mean force and another conditional expectation locally along the string, and therefore can be applied even if the number of collective variables kept in the free energy calculation is large. This is in contrast with other free energy sampling techniques which aim at mapping the full free energy landscape and whose cost increases exponentially with the number of collective variables kept in the free energy. Provided that the number of collective variables is large enough, the new technique captures the mechanism of transition in that it allows to determine the committor function for the reaction and, in particular, the transition state region. The new technique is illustrated on the example of alanine dipeptide, in which we compute the minimum free energy path for the isomerization transition using either two or four dihedral angles as collective variables. It is shown that the mechanism of transition can be captured using the four dihedral angles, but it cannot be captured using only two of them.     

Metadynamics in essential coordinates: free energy simulation of conformational changes

We propose an approach that combines an extraction of collective motions of a molecular system with a sampling of its free energy surface. A recently introduced method of metadynamics allows exploration of the free energy surface of a molecular system by means of coarse-grained dynamics with flooding of free energy minima. This free energy surface is defined as a function of a set of collective variables (e.g., interatomic distances, angles, torsions, and others). In this study, essential coordinates determined by essential dynamics (principle component analysis) were used as collective variables in metadynamics. First, dynamics of the model system (explicitly solvated alanine dipeptide, Ace-Ala-Nme) was simulated by a classical molecular dynamics simulation. The trajectory (1 ns) was then analyzed by essential dynamics to obtain essential coordinates. The free energy surface as a function of the first and second essential coordinates was then explored by metadynamics. The resulting free energy surface is in agreement with other studies of this system. We propose that a combination of these two methods (metadynamics and essential dynamics) has great potential in studies of conformational changes in peptides and proteins.     

A limit of stability in supercooled liquid clusters

We examine the metastable liquid phase of a supercooled gold nanocluster by studying the free energy landscape using the largest solidlike embryo as an order parameter. Just below freezing, the free energy exhibits a local minimum at small embryo sizes and a maximum at a larger critical embryo size. At T=660 K the free energy becomes a monotonically decreasing function of the order parameter as the liquid phase becomes unstable, indicating that we have reached a limit of stability. In contrast to the mean-field theory predictions for a spinodal, the size of the critical embryo remains finite as the limit of stability is approached. We also calculate the rate of nucleation, independently from our free energy calculations, and observe a rapid increase in its temperature dependence when the free energy barrier is on the order of kT. We suggest that this supports the idea that freezing becomes a barrierless process at low temperatures.     

Free energy landscapes for homogeneous nucleation of ice for a monatomic water model

We simulate the homogeneous nucleation of ice from supercooled liquid water at 220 K in the isobaric-isothermal ensemble using the MW monatomic water potential. Monte Carlo simulations using umbrella sampling are performed in order to determine the nucleation free energy barrier. We find the Gibbs energy profile to be relatively consistent with that predicted by classical nucleation theory; the free energy barrier to nucleation was determined to be ~18 k(B)T and the critical nucleus comprised ~85 ice particles. Growth from the supercooled liquid gives clusters that are predominantly cubic, whilst starting with a pre-formed subcritical nucleus of cubic or hexagonal ice results in the growth of predominantly that phase of ice only.     

An efficient method for computing excess free energy of liquid

We present a new theoretical method for efficient calculation of free energy of liquid. This interaction entropy method allows one to compute entropy and free energy of liquid from standard single step MD (molecular dynamics) simulation directly in liquid state without the need to perform MD simulations at many intermediate states as required in thermodynamic integration or free energy perturbation methods. In this new approach, one only needs to evaluate the interaction energy of a single (fixed) liquid molecule with the rest of liquid molecules as a function of time from a standard MD simulation of liquid and the fluctuation of distribution of this interaction energy is then used to calculate the interaction entropy of the liquid. Explicit theoretical derivation of this interaction entropy approach is provided and numerical calculations for the benchmark liquid water system were carried out using three different water models. Numerical analysis of the result was performed and comparison of the computational result with experimental data and other theoretical results were provided. Excellent agreement of calculated free energies with the experimental data using TIP4P model is obtained for liquid water.     

Spreading and retraction dynamics of smectic liquid crystal domains at interfaces

The spreading of a liquid drop over liquid subphase can be driven by change in interfacial tension mediated through a surfactant, volatile solvent or photoinduced reaction. In contrast to the spreading dynamics of a liquid drop, a liquid crystal drop with anisotropic structure can lead to interesting behaviour due to its viscoelasticity and anchoring at the interfaces. Recently, we have reported studies on unusual spreading and retraction dynamics of a smectic domain doped with a fluorescent dye in the collapsed state of a Langmuir monolayer. Under epifluorescence microscope, during excitation, a stack of layers of the dye-doped smectic domain gets sheared causing the domain to spread asymmetrically. Further, due to line tension, the domain transforms into a circular shape. We also find the domain size to be about twice that of the initial size. Interestingly, in the absence of excitation, the domain retracts to a smaller area. During retraction of the domain, successive generation of edge dislocation loops arising from a nucleus results in an increase in the domain thickness. The dynamics of spreading and retraction of the domain can be understood by invoking changes in the spreading coefficient due to photoinduced modification of the interfacial tension.     

A new procedure for analyzing the nucleation kinetics of freezing in computer simulation

A new method for deriving the size of the critical nucleus and the Zeldovich factor directly from kinetic data is presented. Moreover, in principle, the form of G(n), the free energy of formation of nuclei consisting of n molecules, can be inferred. The method involves measuring times of first appearance of nuclei of size n in the transient regime and applying the Becker-Doring theory. Times of first appearance exhibit the same characteristics as the conventional times associated with N(n,t), the number of nuclei of at least size n per unit volume that have materialized at time t. That is, they are well represented by three nucleation parameters, the reduced moment, the time lag, and the steady state nucleation rate. But unlike the conventional steady state rate which is independent of n, the steady state times of first appearance vary with n. In order to characterize the three nucleation parameters with precision, however, thousands of independent stochastic events with known n are required. Such sets of data are readily generated in molecular dynamic simulations but, so far, not in laboratory experiments. Results are illustrated by an analysis of simulations of the spontaneous freezing of large clusters of SeF6.     

Revisiting the Frenkel-Ladd method to compute the free energy of solids: the Einstein molecule approach

In this paper a new method to evaluate the free energy of solids is proposed. The method can be regarded as a variant of the method proposed by Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. The main equations of the method can be derived in a simple way. The method can be easily implemented within a Monte Carlo program. We have applied the method to determine the free energy of hard spheres in the solid phase for several system sizes. The obtained free energies agree within the numerical uncertainty with those obtained by Polson et al. [J. Chem. Phys. 112, 5339 (2000)]. The fluid-solid equilibria has been determined for several system sizes and compared to the values published previously by Wilding and Bruce [Phys. Rev. Lett. 85, 5138 (2000)] using the phase switch methodology. It is shown that both the free energies and the coexistence pressures present a strong size dependence and that the results obtained from free energy calculations agree with those obtained using the phase switch method, which constitutes a cross-check of both methodologies. From the results of this work we estimate the coexistence pressure of the fluid-solid transition of hard spheres in the thermodynamic limit to be p*=11.54(4), which is slightly lower than the classical value of Hoover and Ree (p*=11.70) [J. Chem. Phys. 49, 3609 (1968)]. Taking into account the strong size dependence of the free energy of the solid phase, we propose to introduce finite size corrections, which allow us to estimate approximately the free energy of the solid phase in the thermodynamic limit from the known value of the free energy of the solid phase with N molecules. We have also determined the free energy of a Lennard-Jones solid by using both the methodology of this work and the finite size correction. It is shown how a relatively good estimate of the free energy of the system in the thermodynamic limit is obtained even from the free energy of a relatively small system.     

Collective oscillations and the linear and two-dimensional infrared spectra of inhomogeneous beta-sheets

We numerically calculate the collective amide I oscillations and the associated linear and two-dimensional infrared (2DIR) spectra for model antiparallel beta-sheets and study the effect of inhomogeneity. To visualize the collective vibrational exciton states, a new method is introduced, which proves very useful in classifying the optically dominant states with respect to their symmetry properties and phase relations, even in the absence of exact symmetries. We find that energy (diagonal) and interaction (off-diagonal) disorder may have profoundly different effects on the main peaks in the linear spectrum. We also show that in the 2DIR spectra energy disorder leads to diagonal stretching of the diagonal peaks, while the cross-peaks are typically stretched more horizontally. This offers an explanation for the recently observed overall Z-shape in experimental spectra. Finally, we find that the anharmonic splitting between associated positive and negative features in the 2DIR spectra scales inversely proportionally with the exciton delocalization size imposed by the disorder, thus offering a spectroscopic ruler for this size.