Analysis of probability of inserting a hard spherical particle with small diameter in hard-sphere fluid.

The probability of inserting, without overlap, a hard spherical particle of diameter σ in a hard-sphere fluid of diameter σ0 and packing fraction η determines its excess chemical potential at infinite dilution, µex(σ, η). In our previous work [R. L. Davidchack and B. B. Laird, J. Chem. Phys. 157, 074701 (2022)], we used Widom's particle insertion method within molecular dynamics simulations to obtain high precision results for µex(σ, η) with σ/σ0 ≤ 4 and η ≤ 0.5. In the current work, we investigate the behavior of this quantity at small σ. In particular, using the inclusion-exclusion principle, we relate the insertion probability to the hard-sphere fluid distribution functions and thus derive the higher-order terms in the Taylor expansion of µex(σ, η) at σ = 0. We also use direct evaluation of the excluded volume for pairs and triplets of hard spheres to obtain simulation results for µex(σ, η) at σ/σ0 ≤ 0.2247 that are of much higher precision than those obtained earlier with Widom's method. These results allow us to improve the quality of the small-σ correction in the empirical expression for µex(σ, η) presented in our previous work.

Chemical potential and surface free energy of a hard spherical particle in hard-sphere fluid over the full range of particle diameters.

The excess chemical potential µex(σ, η) of a test hard spherical particle of diameter σ in a fluid of hard spheres of diameter σ0 and packing fraction η can be computed with high precision using Widom's particle insertion method [B. Widom, J. Chem. Phys. 39, 2808 (1963)] for σ between 0 and just larger than 1 and/or small η. Heyes and Santos [J. Chem. Phys. 145, 214504 (2016)] analytically showed that the only polynomial representation of µex consistent with the limits of σ at zero and infinity has a cubic form. On the other hand, through the solvation free energy relationship between µex and the surface free energy γ of hard-sphere fluids at a hard spherical wall, we can obtain precise measurements of µex for large σ, extending up to infinity (flat wall) [R. L. Davidchack and B. B. Laird, J. Chem. Phys. 149, 174706 (2018)]. Within this approach, the cubic polynomial representation is consistent with the assumptions of morphometric thermodynamics. In this work, we present the measurements of µex that combine the two methods to obtain high-precision results for the full range of σ values from zero to infinity, which show statistically significant deviations from the cubic polynomial form. We propose an empirical functional form for the µex dependence on σ and η, which better fits the measurement data while remaining consistent with the analytical limiting behavior at zero and infinite σ.

Cleaving Method for Molecular Crystals and Its Application to Calculation of the Surface Free Energy of Crystalline ß-d-Mannitol at Room Temperature.

Calculation of the surface free energy (SFE) is an important application of the thermodynamic integration (TI) methodology, which was mainly employed for atomic crystals (such as Lennard-Jones or metals). In this work, we present the calculation of the SFE of a molecular crystal using the cleaving technique which we implemented in the LAMMPS molecular dynamics package. We apply this methodology to a crystal of ß-d-mannitol at room temperature and report the results for two types of force fields belonging to the GROMOS family: all atoms and united atoms. The results show strong dependence on the type of force field used, highlighting the need for the development of better force fields to model the surface properties of molecular crystals. In particular, we observe that the united-atoms force field, despite its higher degree of coarse graining compared to the all-atoms force field, produces SFE results in better agreement with the experimental results from inverse gas chromatography measurements.

Shuttleworth equation: A molecular simulations perspective.

Even though the study of interfacial phenomena can be traced back to Laplace and was given a solid thermodynamic foundation by Gibbs, it appears that some concepts and relations among them are still causing some confusion and debates in the literature, particularly for interfaces involving solids. In particular, the definitions of the concepts of interfacial tension, free energy, and stress and the relationships between them sometimes lack clarity, and some authors even question their validity. So far, the debates about these relationships, in particular the Shuttleworth equation, have taken place within the framework of classical thermodynamics. In this work, we are offering to look at these concepts within the framework of statistical mechanics, which can be readily tested in Molecular Dynamics (MD) simulations. For a simple one component system of particles interacting via the Lennard-Jones potential, we calculate by the cleaving method the excess free energy of a solid-vacuum interface (solid surface) for systems in different states of tangential strain and compare the results to the calculation of surface stress via the difference of normal and tangential forces at the surface. As a result, we demonstrate consistency, within the statistical uncertainty, of the thermodynamic and statistical mechanical definitions of surface free energy and surface stress and how they are expressed via interaction-dependent quantities in MD simulations.

Surface free energy of a hard-sphere fluid at curved walls: Deviations from morphometric thermodynamics.

We report molecular-dynamics (MD) simulation results for the surface free energy of a hard-sphere fluid at cylindrical and spherical hard walls of different radii. The precision of the results is much higher than that in our previous study [B. B. Laird et al., Phys. Rev. E 86, 060602 (2012)], allowing us to estimate the size of deviations from the predictions of Morphometric Thermodynamics (MT). We compare our results to the analytical expressions for the surface energy as a function of wall radius R and fluid density derived from the White Bear II variant of the density functional theory, as well as to the leading terms of the virial expansion. For the cylindrical wall, we observe deviations from MT proportional to R -2 and R -3, which are consistent with the available virial expressions. For the spherical wall, while the precision is not sufficient to detect statistically significant deviations from MT, the MD results indicate the range of densities for which the truncated virial expansions are applicable.

Properties of the hard-sphere fluid at a planar wall using virial series and molecular-dynamics simulation.

We study the hard-sphere fluid in contact with a planar hard wall. By combining the inhomogeneous virial series with simulation results, we achieve a new benchmark of accuracy for the calculation of surface thermodynamics properties such as surface adsorption Γ and the surface free energy (or surface tension), γ. We briefly introduce the problem of choosing a position for the dividing surface and avoid it by proposing the use of alternative functions to Γ and γ that are independent of the adopted frame of reference. Finally, we present analytic expressions for the dependence of system surface thermodynamic properties on packing fraction, ensuring the high accuracy of the parameterized functions for any frame of reference. The proposed parametric expressions for both, Γ and γ, fit the accurate simulation results within the statistical error.

Reduction of SO(2) symmetry for spatially extended dynamical systems.

Spatially extended systems, such as channel or pipe flows, are often equivariant under continuous symmetry transformations, with each state of the flow having an infinite number of equivalent solutions obtained from it by a translation or a rotation. This multitude of equivalent solutions tends to obscure the dynamics of turbulence. Here we describe the "first Fourier mode slice," a very simple, easy to implement reduction of SO(2) symmetry. While the method exhibits rapid variations in phase velocity whenever the magnitude of the first Fourier mode is nearly vanishing, these near singularities can be regularized by a time-scaling transformation. We show that after application of the method, hitherto unseen global structures, for example, Kuramoto-Sivashinsky relative periodic orbits and unstable manifolds of traveling waves, are uncovered.

A molecular dynamics study of Young's modulus change of semi-crystalline polymers during degradation by chain scissions.

This paper presents a molecular dynamics study on the change in Young's modulus of semi-crystalline polymers during degradation by chain scissions, which is relevant to the study of mechanical properties of biodegrading polymers. Using a simple polymer model whose structural and mechanical properties are similar to that of a commonly used biodegrading polymer poly(glycolic acid), we combine molecular dynamics and Monte Carlo to model a system of two polymer crystals separated by an amorphous region between them. The polymer chains in the amorphous region are cut randomly to mimic hydrolysis chain scissions. In a series of virtual tensile tests, the systems with various numbers of chain scissions are subjected to a unidirectional deformation. We find that at temperatures below the glass transition temperature of the model polymer, the Young's modulus of the system reduces quickly with the number of chain scissions, while at temperatures above the glass transition temperature, the Young's modulus reduction lags behind the polymer chain scissions. This observation supports the entropy-spring model of amorphous polymers proposed by Wang et al., which suggests that Young's modulus above the glass transition temperature is dominated by the internal energy of the system, while below the glass transition temperature it is dominated by the entropy of the amorphous phase. The numerical study therefore provides a molecular understanding of the widely observed behaviours of semi-crystalline biodegradable polymers.

Ice Ih-Water Interfacial Free Energy of Simple Water Models with Full Electrostatic Interactions.

We employ the cleaving approach to calculate directly the ice Ih-water interfacial free energy for the simple models of water, TIP4P, TIP4P-Ew, and TIP5P-E, with full electrostatic interactions evaluated via the Ewald sums. The results are in good agreement with experimental values, but lower than previously obtained for TIP4P-Ew and TIP5P-E by indirect methods. We calculate the interfacial free energies for basal, prism, and {112̅0} interfaces and find that the anisotropy of the TIP5P-E model is different from that of the TIP4P models. The effect of including full electrostatic interactions is determined to be smaller than 10% compared to the water models with damped Coulomb interactions, which indicates that the value of the ice-water interfacial free energy is determined predominantly by the short-range packing interaction between water molecules. We also observe a strong linear correlation between the interfacial free energy and the melting temperature of different water models.

Interfacial free energy of a hard-sphere fluid in contact with curved hard surfaces.

Using molecular-dynamics simulation, we have calculated the interfacial free energy γ between a hard-sphere fluid and hard spherical and cylindrical colloidal particles, as functions of the particle radius R and the fluid packing fraction η=ρσ(3)/6, where ρ and σ are the number density and hard-sphere diameter, respectively. These results verify that Hadwiger's theorem from integral geometry, which predicts that γ for a fluid at a surface, with certain restrictions, should be a linear combination of the average mean and Gaussian surface curvatures, is valid within the precision of the calculation for spherical and cylindrical surfaces up to η ≈ 0.42. In addition, earlier results for γ for this system [Bryk et al., Phys. Rev. E 68, 031602 (2003)] using a geometrically based classical density functional theory are in excellent agreement with the current simulation results for packing fractions in the range where Hadwiger's theorem is valid. However, above η ≈ 0.42, γ(R) shows significant deviations from the Hadwiger form indicating limitations to its use for high-density hard-sphere fluids. Using the results of this study together with Hadwiger's theorem allows one, in principle, to determine γ for any sufficiently smooth surface immersed in a hard-sphere fluid.

Hard spheres revisited: accurate calculation of the solid-liquid interfacial free energy.

We revise the earlier [R. L. Davidchack and B. B. Laird, Phys. Rev. Lett. 85, 4751 (2000)] direct calculation of the hard sphere solid-liquid interfacial free energy by the cleaving walls method. The revisions of the method include modified interactions with the cleaving walls and the use of a nonequilibrium work measurements approach, which allows for a more robust control of the accuracy of the obtained results. We find that the new values are lower compared to the original ones, which is consistent with the more recent indirect estimates based on extrapolation from the soft-sphere results [R. L. Davidchack and B. B. Laird, Phys. Rev. Lett. 94, 086102 (2005)], as well as those obtained using the capillary fluctuations method [R. L. Davidchack, J. R. Morris, and B. B. Laird, J. Chem. Phys. 125, 094710 (2006)].

Calculation of the interfacial free energy of a fluid at a static wall by Gibbs-Cahn integration.

The interface between a fluid and a static wall is a useful model for a chemically heterogeneous solid-liquid interface. In this work, we outline the calculation of the wall-fluid interfacial free energy (gamma(wf)) for such systems using molecular simulation combined with adsorption equations based on Cahn's extension of the surface thermodynamics of Gibbs. As an example, we integrate such an adsorption equation to obtain gamma(wf) as a function of pressure for a hard-sphere fluid at a hard wall. The results so obtained are shown to be in excellent agreement in both magnitude and precision with previous calculations of this quantity, but are obtained with significantly lower computational effort.

On the use of stabilizing transformations for detecting unstable periodic orbits in high-dimensional flows.

We explore the possibility of extending the stabilizing transformations approach [J. J. Crofts and R. L. Davidchack, SIAM J. Sci. Comput. (USA) 28, 1275 (2006)]. to the problem of locating large numbers of unstable periodic orbits in high-dimensional flows, in particular those that result from spatial discretization of partial differential equations. The approach has been shown to be highly efficient when detecting large sets of periodic orbits in low-dimensional maps. Extension to low-dimensional flows has been achieved by the use of an appropriate Poincare surface of section [D. Pingel, P. Schmelcher, and F. K. Diakonos, Phys. Rep. 400, 67 (2004)]. For the case of high-dimensional flows, we show that it is more efficient to apply stabilizing transformations directly to the flows without the use of the Poincare surface of section. We use the proposed approach to find many unstable periodic orbits in the model example of a chaotic spatially extended system-the Kuramoto-Sivashinsky equation. The performance of the proposed method is compared against other methods such as Newton-Armijo and Levenberg-Marquardt algorithms. In the latter case, we also argue that the Levenberg-Marquardt algorithm, or any other optimization-based approach, is more efficient and simpler in implementation when applied directly to the detection of periodic orbits in high-dimensional flows without the use of the Poincare surface of section or other additional constraints.

Determination of the solid-liquid interfacial free energy along a coexistence line by Gibbs-Cahn integration.

We calculate the solid-liquid interfacial free energy gamma(sl) for the Lennard-Jones (LJ) system at several points along the pressure-temperature coexistence curve using molecular-dynamics simulation and Gibbs-Cahn integration. This method uses the excess interfacial energy (e) and stress (tau) along the coexistence curve to determine a differential equation for gamma(sl) as a function of temperature. Given the values of gamma(sl) for the (100), (110), and (111) LJ interfaces at the triple-point temperature (T( *)=kT/varepsilon=0.618), previously obtained using the cleaving method by Davidchack and Laird [J. Chem. Phys. 118, 7657 (2003)], this differential equation can be integrated to obtain gamma(sl) for these interfaces at higher coexistence temperatures. Our values for gamma(sl) calculated in this way at T( *)=1.0 and 1.5 are in good agreement with those determined previously by cleaving, but were obtained with significantly less computational effort than required by either the cleaving method or the capillary fluctuation method of Hoyt, Asta, and Karma [Phys. Rev. Lett. 86, 5530 (2001)]. In addition, the orientational anisotropy in the excess interface energy, stress and entropy, calculated using the conventional Gibbs dividing surface, are seen to be significantly larger than the relatively small anisotropies in gamma(sl) itself.

Langevin thermostat for rigid body dynamics.

We present a new method for isothermal rigid body simulations using the quaternion representation and Langevin dynamics. It can be combined with the traditional Langevin or gradient (Brownian) dynamics for the translational degrees of freedom to correctly sample the canonical distribution in a simulation of rigid molecules. We propose simple, quasisymplectic second-order numerical integrators and test their performance on the TIP4P model of water. We also investigate the optimal choice of thermostat parameters.

Direct calculation of solid-liquid interfacial free energy for molecular systems: TIP4P ice-water interface.

By extending the cleaving method to molecular systems, we perform direct calculations of the ice Ih-water interfacial free energy for the TIP4P model. The values for the basal, prism, and {112[over]0} faces are 23.3+/-0.8 mJ m{-2}, 23.6+/-1.0 mJ m{-2}, and 24.7+/-0.8 mJ m{-2}, respectively. The closeness of these values implies a minimal role of thermodynamic factors in the anisotropic growth of ice crystals. These results are about 20% lower than the best experimental estimates. However, the Turnbull coefficient is about 50% higher than for real water, indicating a possible limitation of the TIP4P model in describing freezing.

The anisotropic hard-sphere crystal-melt interfacial free energy from fluctuations.

We have calculated the interfacial free energy for the hard-sphere system, as a function of crystal interface orientation, using a method that examines the fluctuations in the height of the interface during molecular dynamics simulations. The approach is particularly sensitive for the anisotropy of the interfacial free energy. We find an average interfacial free energy of gamma=0.56+/-0.02k(B)Tsigma(-2). This value is lower than earlier results based upon direct calculations of the free energy [R. L. Davidchack and B. B. Laird, Phys. Rev. Lett. 85, 4751 (2000)]. However, both the average value and the anisotropy agree with the recent values obtained by extrapolation from direct calculations for a series of the inverse-power potentials [R. L. Davidchack and B. B. Laird, Phys. Rev. Lett. 94, 086102 (2005)].

Crystal structure and interaction dependence of the crystal-melt interfacial free energy.

We examine via molecular simulation the dependence of the crystal-melt interfacial free energy gamma on molecular interaction and crystal structure (fcc vs bcc) for systems interacting with inverse-power repulsive potentials, u(r)=epsilon(sigma/r)(n), 6< or =n< or =100. Both the magnitude and anisotropy of gamma are found to increase as the range of the potential increases. Also we find that gamma(bcc)

Direct calculation of the crystal-melt interfacial free energy via molecular dynamics computer simulation.

We review our recent work on the direct calculation of the interfacial free energy, gamma, of the crystal-melt interface via molecular dynamics computer simulation for a number of model systems. The value of gamma as a function of crystal orientation is determined using a thermodynamic integration technique employing moving cleaving walls [Phys. Rev. Lett. 2000, 85, 4751]. The calculation is sufficiently precise to resolve the small anisotropy in gamma, which is crucial in determining the kinetics and morphology of dendritic growth. We report values of gamma for the hard-sphere and Lennard-Jones systems, as well as recent results on the series of inverse-power potentials. For the inverse sixth-, seventh-, and eighth-power systems, we determine gamma for both fcc and bcc crystal structures. For these systems, the bcc-melt gamma is lower than that for fcc crystals by about 25%, consistent with recent experiments and computer simulations on fcc-forming systems that show preferential formation of bcc nuclei in the initial stages of crystallization. In addition, we show that our results give a molecular interpretation of Turnbull's rule, which is the empirical relationship between gamma and the enthalpy of fusion.