Measure-preserving integrators for molecular dynamics in the isothermal–isobaric ensemble derived from the Liouville operator

Tracking conformational transitions of the gonadotropin hormone receptors in a bilayer of (SDPC) poly-unsaturated lipids from all-atom molecular dynamics simulations

Glycoprotein hormone receptors [thyrotropin (TSHR), luteinizing hormone/chorionic gonadotropin (LHCGR), and follicle stimulating hormone (FSHR) receptors] are rhodopsin-like G protein-coupled receptors. These receptors display common structural features including a prominent extracellular domain with leucine-rich repeats (LRR) stabilized by β-sheets and a long and flexible loop known as the hinge region (HR), and a transmembrane (TM) domain with seven α-helices interconnected by intra- and extracellular loops. Binding of the ligand to the LRR resembles a hand coupling transversally to the α- and β-subunits of the hormone, with the thumb being the HR. The structure of the FSH-FSHR complex suggests an activation mechanism in which Y335 at the HR binds into a pocket between the α- and β-chains of the hormone, leading to an adjustment of the extracellular loops. In this study, we performed molecular dynamics (MD) simulations to identify the conformational changes of the FSHR and LHCGR. We set up a FSHR structure as predicted by AlphaFold (AF-P23945); for the LHCGR structure we took the cryo-electron microscopy structure for the active state (PDB:7FII) as initial coordinates. Specifically, the flexibility of the HR domain and the correlated motions of the LRR and TM domain were analyzed. From the conformational changes of the LRR, TM domain, and HR we explored the conformational landscape by means of MD trajectories in all-atom approximation, including a membrane of polyunsaturated phospholipids. The distances and procedures here defined may be useful to propose reaction coordinates to describe diverse processes, such as the active-to-inactive transition, and to identify intermediaries suited for allosteric regulation and biased binding to cellular transducers in a selective activation strategy.

Fullerenes containing water molecules: a study of reactive molecular dynamics simulations

A different technique was used to investigate fullerenes encapsulating a polar guest species. By reactive molecular dynamics simulations, three types of fullerenes were investigated on a gold surface: an empty C60, a single H2O molecule inside C60 (H2O@C60), and two water molecules inside C60 ((H2O)2@C60). Our findings revealed that despite the free movement of all fullerenes on gold surfaces, confined H2O molecules within the fullerenes result in a distinct pattern of motion in these systems. The (H2O)2@C60 complex had the highest displacement and average velocity, while C60 had the lowest displacement and average velocity. The symmetry of molecules and the polarity of water seem to be crucial in these cases. ReaxFF simulations showed that water molecules in an H2O molecule, H2O@C60, and (H2O)2@C60 have dipole moments of 1.76, 0.42, and 0.47 D, respectively. A combination of the non-polar C60 and polar water demonstrated a significant reduction in the dipole moment of H2O molecules due to encapsulation. The dipole moments of water molecules agreed with those in other studies, which can be useful in the development of biocompatible and high-efficiency nanocars.

The Crucial Role of Solvation Forces in the Steric Stabilization of Nanoplatelets

The steric stability of inorganic colloidal particles in an apolar solvent is usually described in terms of the balance between three contributions: the van der Waals attraction, the free energy of mixing, and the ligand compression. However, in the case of nanoparticles, the discrete nature of the ligand shell and the solvent has to be taken into account. Cadmium selenide nanoplatelets are a special case. They combine a weak van der Waals attraction and a large facet to particle size ratio. We use coarse grained molecular dynamics simulations of nanoplatelets in octane to demonstrate that solvation forces are strong enough to induce the formation of nanoplatelet stacks and by that have a crucial impact on the steric stability. In particular, we demonstrate that for sufficiently large nanoplatelets, solvation forces are proportional to the interacting facet area, and their strength is intrinsically tied to the softness of the ligand shell.

Investigation of fullerene motion on thermally activated gold substrates with different shapes

In the current study, the regime of motion of fullerene molecules on substrates with different shapes at a range of specific temperatures has been investigated. To do so, the potential energy of fullerene molecules was analyzed using the classical molecular dynamics method. C₂₀, C₃₆, C₅₀, C₆₀, C₇₂, C₇₆, C₈₀, and C₉₀ fullerene molecules were selected due to their spherical shapes with different sizes. In addition, to completely analyze the behavior of these molecules, different gold substrates, including flat, concave, the top side of the step (upward step), and the downside of the step (downward step) substrates, were considered. Specifying the regime of the motion at different temperatures is one of the main goals of this study. For this purpose, we have studied the translational and rotational motions of fullerene molecules independently. In the first step of the investigation, Lennard-Jones potential energy of fullerene molecules was calculated. Subsequently, the regime of motion of different fullerenes has been classified, based on their displacement and sliding velocity. Our findings indicated that C₆₀ is appropriate in less than [Formula: see text] of the conditions. However, C₂₀, C₇₆ and C₈₀ molecules were found to be appropriate candidates in most cases in different conditions while they were incompetent only in seven situations. As far as a straight-line movement is considered, the concave geometry demonstrated a better performance compared to the other substrates. In addition, C₇₂ indicated less favorable performance concerning the range of movement and diffusion coefficients. All in all, our investigation helps to understand the performance of different fullerene molecules on gold substrates and find their probable application, especially as a wheel in nano-machine structures.

Molecular Dynamics of Solids at Constant Pressure and Stress Using Anisotropic Stochastic Cell Rescaling

Molecular dynamics simulations of solids are often performed using anisotropic barostats that allow the shape and volume of the periodic cell to change during the simulation. Most existing schemes are based on a second-order differential equation that might lead to undesired oscillatory behaviors and should not be used in the equilibration phase. We recently introduced stochastic cell rescaling, a first-order stochastic barostat that can be used for both the equilibration and production phases. Only the isotropic and semi-isotropic variants have been formulated and implemented so far. In this paper, we develop and implement the equations of motion of the fully anisotropic variant and test them on Lennard-Jones solids, ice, gypsum, and gold. The algorithm has a single parameter that controls the relaxation time of the volume, results in the exponential decay of correlation functions, and can be effectively applied to a wide range of systems.

Dynamics and phase behavior of two-dimensional size-asymmetric binary mixtures of core-softened colloids

The self-assembly of binary colloidal mixtures provides a bottom-up approach to create novel functional materials. To elucidate the effect of composition, temperature, and pressure on the self-assembly behavior of size-asymmetric mixtures, we performed extensive dynamics simulations of a simple model of polymer-grafted colloids. We have used a core-softened interaction potential and extended it to represent attractive interactions between unlike colloids and repulsions between like colloids. Our study focused on size-asymmetric mixtures where the ratio between the sizes of the colloidal cores was fixed at σ_Bσ_A=0.5. We have performed extensive simulations in the isothermal-isobaric and canonical (NVT) ensembles to elucidate the phase behavior and dynamics of mixtures with different stoichiometric ratios. Our simulation results uncovered a rich phase behavior, including the formation of hierarchical structures with many potential applications. For compositions where small colloids are the majority, sublattice melting occurs for a wide range of densities. Under these conditions, large colloids form a well-defined lattice, whereas small colloids can diffuse through the system. As the temperature is decreased, the small colloids localize, akin to a metal-insulator transition, with the small colloids playing a role similar to electrons. Our results are summarized in terms of phase diagrams.

Structure factor lineshape model gives approximate nanoscale size of polar aggregates in the ionic liquid N-methyl-N-propylpyrrolidinium bis(trifluoromethylsulfonyl)imide

A Lorentzian lineshape model is developed and tested for the charge alternation peak in X-ray structure factors calculated from MD simulations for N-methyl-N-propylpyrrolidinium bis(trifluoromethylsulfonyl)imide. Applying the model to published, experimental X-ray scattering data reproduces calculated cation-cation and anion-anion distances within 6% and implies that half of ionic aggregates are larger than 12.7 Å.

Locomotion of the C60-based nanomachines on graphene surfaces

We provide a comprehensive computational characterization of surface motion of two types of nanomachines with four C₆₀ "wheels": a flexible chassis Nanocar and a rigid chassis Nanotruck. We study the nanocars' lateral and rotational diffusion as well as the wheels' rolling motion on two kinds of graphene substrates-flexible single-layer graphene which may form surface ripples and an ideally flat graphene monolayer. We find that the graphene surface ripples facilitate the translational diffusion of Nanocar and Nanotruck, but have little effect on their surface rotation or the rolling of their wheels. The latter two types of motion are strongly affected by the structure of the nanomachines instead. Surface diffusion of both nanomachines occurs preferentially via a sliding mechanism whereas the rolling of the "wheels" contributes little. The axial rotation of all "wheels" is uncorrelated.

Pressure control using stochastic cell rescaling

Molecular dynamics simulations require barostats to be performed at a constant pressure. The usual recipe is to employ the Berendsen barostat first, which displays a first-order volume relaxation efficient in equilibration but results in incorrect volume fluctuations, followed by a second-order or a Monte Carlo barostat for production runs. In this paper, we introduce stochastic cell rescaling, a first-order barostat that samples the correct volume fluctuations by including a suitable noise term. The algorithm is shown to report volume fluctuations compatible with the isobaric ensemble and its anisotropic variant is tested on a membrane simulation. Stochastic cell rescaling can be straightforwardly implemented in the existing codes and can be used effectively in both equilibration and production phases.

The thermodynamics of a liquid-solid interface at extreme conditions: A model close-packed system up to 100 GPa

The first experimental insight into the nature of the liquid-solid interface occurred with the pioneering experiments of Turnbull, which simultaneously demonstrated both that metals could be deeply undercooled (and therefore had relatively large barriers to nucleation) and that the inferred interfacial free energy γ was linearly proportional to the enthalpy of fusion [D. Turnbull, J. Appl. Phys. 21, 1022 (1950)]. By an atomistic simulation of a model face-centered cubic system via adiabatic free energy dynamics, we extend Turnbull's result to the realm of high pressure and demonstrate that the interfacial free energy, evaluated along the melting curve, remains linear with the bulk enthalpy of fusion, even up to 100 GPa. This linear dependence of γ on pressure is shown to be a consequence of the entropy dominating the free energy of the interface in conjunction with the fact that the entropy of fusion does not vary greatly along the melting curve for simple monoatomic metals. Based on this observation, it appears that large undercoolings in liquid metals can be achieved even at very high pressure. Therefore, nucleation rates at high pressure are expected to be non-negligible, resulting in observable solidification kinetics.

Stochastic sampling of the isothermal-isobaric ensemble: Phase diagram of crystalline solids from molecular dynamics simulation

A methodology to sample the isothermal-isobaric ensemble using Langevin dynamics is proposed, which combines novel features of geometric integrators for the equations of motion. By employing the Trotter expansion, the methodology generates a robust, symmetric, and accurate numerical algorithm. In order to show that the proposed method correctly samples the phase-space, simulations in the isotropic NPT ensemble were carried out for two analytical examples. Also this method lets us study a solid-solid phase transition by conducting a fully flexible-cell molecular dynamics simulation. Additionally, we present an efficient method to determine the Gibbs free energy in a wide interval of pressure along an isothermal path, which allows us to determine the transition pressure in a driven by pressure solid-solid phase transition. Our calculations show that the methodology is highly suitable for the study of the phase diagram of crystalline solids.

Structures of FEC-containing electrolytes and the stabilization mechanism at high voltage and elevated temperature

The performance of lithium-ion batteries is strongly dependent on the properties of electrolytes.

Phase diagram of Janus particles: The missing dimension of pressure anisotropy

Brownian dynamics simulations of single-patch Janus particles under sedimentation equilibrium reveal that the phases found at fixed temperature and volume fraction are extremely sensitive to small changes in lateral box dimension. We trace this sensitivity to an uncontrolled parameter, namely, the pressure component parallel to the hexagonally ordered layers formed through sedimentation. We employ a flexible-cell constant-pressure scheme to achieve explicit control over this usually overlooked parameter, enabling the estimation of phase behavior under given pressure anisotropy. Our results show an increase in the stability range of an orientationally ordered lamellar phase with lateral layer compression and suggest a novel mechanism to control solid-solid phase transitions with negligible change in system volume, thus showing prospect for design of novel structures and switchable crystals from anisotropic building blocks.

Evaluation of atomic pressure in the multiple time-step integration algorithm

In molecular dynamics (MD) calculations, reduction in calculation time per MD loop is essential. A multiple time-step (MTS) integration algorithm, the RESPA (Tuckerman and Berne, J. Chem. Phys. 1992, 97, 1990-2001), enables reductions in calculation time by decreasing the frequency of time-consuming long-range interaction calculations. However, the RESPA MTS algorithm involves uncertainties in evaluating the atomic interaction-based pressure (i.e., atomic pressure) of systems with and without holonomic constraints. It is not clear which intermediate forces and constraint forces in the MTS integration procedure should be used to calculate the atomic pressure. In this article, we propose a series of equations to evaluate the atomic pressure in the RESPA MTS integration procedure on the basis of its equivalence to the Velocity-Verlet integration procedure with a single time step (STS). The equations guarantee time-reversibility even for the system with holonomic constrants. Furthermore, we generalize the equations to both (i) arbitrary number of inner time steps and (ii) arbitrary number of force components (RESPA levels). The atomic pressure calculated by our equations with the MTS integration shows excellent agreement with the reference value with the STS, whereas pressures calculated using the conventional ad hoc equations deviated from it. Our equations can be extended straightforwardly to the MTS integration algorithm for the isothermal NVT and isothermal-isobaric NPT ensembles. © 2017 Wiley Periodicals, Inc.

Recent advances in the computational chemistry of soft porous crystals

We highlight recent progress in the field of computational chemistry of nanoporous materials, focusing on methods and studies that address the extraordinary dynamic nature of these systems: the high flexibility of their frameworks, the large-scale structural changes upon external physical or chemical stimulation, and the presence of defects and disorder.

Rotator-to-Lamellar Phase Transition in Janus Colloids Driven by Pressure Anisotropy

We demonstrate through Brownian dynamics simulations a phase transition in plastic crystalline assemblies of Janus spheres through controlled pressure anisotropy. When the pressure in plane with hexagonally ordered layers is increased relative to that normal to the layers, a rapid first-order rotator-to-lamellar transition of Janus sphere orientation occurs at constant temperature. We show that the underlying mechanism closely follows the Maier-Saupe theory, originally developed for isotropic-to-nematic transition in positionally disordered materials but here applied to positionally ordered ones. Since the transition involves almost no translational diffusion or volume change, and occurs rapidly by particle rotation, the results should help guide the design of rapidly switchable colloidal crystals.

Exploring polymorphism of benzene and naphthalene with free energy based enhanced molecular dynamics

Prediction and exploration of possible polymorphism in organic crystal compounds are of great importance for industries ranging from organic electronics to pharmaceuticals to high-energy materials. Here we apply our crystal structure prediction procedure and the enhanced molecular dynamics based sampling approach called the Crystal-Adiabatic Free Energy Dynamics (Crystal-AFED) method to benzene and naphthalene. Crystal-AFED allows the free energy landscape of structures to be explored efficiently at any desired temperature and pressure. For each system, we successfully predict the most stable crystal structures at atmospheric pressure and explore the relative Gibbs free energies of predicted polymorphs at high pressures. Using Crystal-AFED sampling, we find that mixed structures, which typically cannot be discovered by standard crystal structure prediction methods, are prevalent in the solid forms of these compounds at high pressure.

Order-parameter-aided temperature-accelerated sampling for the exploration of crystal polymorphism and solid-liquid phase transitions

The problem of predicting polymorphism in atomic and molecular crystals constitutes a significant challenge both experimentally and theoretically. From the theoretical viewpoint, polymorphism prediction falls into the general class of problems characterized by an underlying rough energy landscape, and consequently, free energy based enhanced sampling approaches can be brought to bear on the problem. In this paper, we build on a scheme previously introduced by two of the authors in which the lengths and angles of the supercell are targeted for enhanced sampling via temperature accelerated adiabatic free energy dynamics [T. Q. Yu and M. E. Tuckerman, Phys. Rev. Lett. 107, 015701 (2011)]. Here, that framework is expanded to include general order parameters that distinguish different crystalline arrangements as target collective variables for enhanced sampling. The resulting free energy surface, being of quite high dimension, is nontrivial to reconstruct, and we discuss one particular strategy for performing the free energy analysis. The method is applied to the study of polymorphism in xenon crystals at high pressure and temperature using the Steinhardt order parameters without and with the supercell included in the set of collective variables. The expected fcc and bcc structures are obtained, and when the supercell parameters are included as collective variables, we also find several new structures, including fcc states with hcp stacking faults. We also apply the new method to the solid-liquid phase transition in copper at 1300 K using the same Steinhardt order parameters. Our method is able to melt and refreeze the system repeatedly, and the free energy profile can be obtained with high efficiency.

Ab initio molecular dynamics study of water at constant pressure using converged basis sets and empirical dispersion corrections

It is generally believed that studies of liquid water using the generalized gradient approximation to density functional theory require dispersion corrections in order to obtain reasonably accurate structural and dynamical properties. Here, we report on an ab initio molecular dynamics study of water in the isothermal-isobaric ensemble using a converged discrete variable representation basis set and an empirical dispersion correction due to Grimme [J. Comp. Chem. 27, 1787 (2006)]. At 300 K and an applied pressure of 1 bar, the density obtained without dispersion corrections is approximately 0.92 g/cm(3) while that obtained with dispersion corrections is 1.07 g/cm(3), indicating that the empirical dispersion correction overestimates the density by almost as much as it is underestimated without the correction for this converged basis. Radial distribution functions exhibit a loss of structure in the second solvation shell. Comparison of our results with other studies using the same empirical correction suggests the cause of the discrepancy: the Grimme dispersion correction is parameterized for use with a particular basis set; this parameterization is sensitive to this choice and, therefore, is not transferable to other basis sets.

Identifying low variance pathways for free energy calculations of molecular transformations in solution phase

Improving the efficiency of free energy calculations is important for many biological and materials design applications, such as protein-ligand binding affinities in drug design, partitioning between immiscible liquids, and determining molecular association in soft materials. We show that for any pair potential, moderately accurate estimation of the radial distribution function for a solute molecule is sufficient to accurately estimate the statistical variance of a sampling along a free energy pathway. This allows inexpensive analytical identification of low statistical error free energy pathways. We employ a variety of methods to estimate the radial distribution function (RDF) and find that the computationally cheap two-body "dilute gas" limit performs as well or better than 3D-RISM theory and other approximations for identifying low variance free energy pathways. With a RDF estimate in hand, we can search for pairwise interaction potentials that produce low variance. We give an example of a search minimizing statistical variance of solvation free energy over the entire parameter space of a generalized "soft core" potential. The free energy pathway arising from this optimization procedure has lower curvature in the variance and reduces the total variance by at least 50% compared to the traditional soft core solvation pathway. We also demonstrate that this optimized pathway allows free energies to be estimated with fewer intermediate states due to its low curvature. This free energy variance optimization technique is generalizable to solvation in any homogeneous fluid and for any type of pairwise potential and can be performed in minutes to hours, depending on the method used to estimate g(r).

A non-polarizable model of water that yields the dielectric constant and the density anomalies of the liquid: TIP4Q

A four-site rigid water model is presented, whose parameters are fitted to reproduce the experimental static dielectric constant at 298 K, the maximum density of liquid water and the equation of state at low pressures. The model has a positive charge on each of the three atomic nuclei and a negative charge located at the bisector of the HOH bending angle. This charge distribution allows increasing the molecular dipole moment relative to four-site models with only three charges and improves the liquid dielectric constant at different temperatures. Several other properties of the liquid and of ice Ih resulting from numerical simulations with the model are in good agreement with experimental values over a wide range of temperatures and pressures. Moreover, the model yields the minimum density of supercooled water at 190 K and the minimum thermal compressibility at 310 K, close to the experimental values. A discussion is presented on the structural changes of liquid water in the supercooled region where the derivative of density with respect to temperature is a maximum.

Temperature-accelerated method for exploring polymorphism in molecular crystals based on free energy

The ability of certain organic molecules to form multiple crystal structures, known as polymorphism, has important ramifications for pharmaceuticals and high energy materials. Here, we introduce an efficient molecular dynamics method for rapidly identifying and thermodynamically ranking polymorphs. The new method employs high temperature and adiabatic decoupling to the simulation cell parameters in order to sample the Gibbs free energy of the polymorphs. Polymorphism in solid benzene is revisited, and a resolution to a long-standing controversy concerning the benzene II structure is proposed.

Constant pressure ab initio molecular dynamics with discrete variable representation basis sets

The use of discrete variable representation (DVR) basis sets within ab initio molecular dynamics calculations allows the latter to be performed with converged energies and, more importantly, converged forces. In this paper, we show how to carry out ab initio molecular dynamics calculations in the isothermal-isobaric ensemble with fully flexible simulation boxes within the DVR basis set framework. In particular, we derive the appropriate DVR based expression for the pressure tensor when the electronic structure is represented using Kohn-Sham density functional theory, and we examine the convergence of this expression as a function of the basis set size. An illustrative example using 64 silicon atoms in a fully flexible box using a combination of the Martyna-Tobias-Klein [Martyna et al., J. Chem. Phys. 101, 4177 (1994)] and Car-Parrinello [Car and Parinello, Phys. Rev. Lett. 55, 2471 (1985)] algorithms is presented to demonstrate the efficacy of the approach.