Mold: a LAMMPS package to compute interfacial free energies and nucleation rates

The determination of the free energy cost of forming an interface between two phases, i


Statement of need
Liquid-to-solid phase transitions are present in many physical and industrial processes such as freezing (Espinosa, Vega, et al., 2016), nucleation (Mer, 1952), casting (Asta et al., 2009) or welding (Asta et al., 2009).A thorough knowledge of the underlying thermodynamic processes and mechanisms is essential for understanding and controlling such phenomena.In particular, the interfacial free energy, defined as the reversible work required to form a unit area of the interface between two coexisting phases, is crucial for determining the conditions that lead to a phase transition.The calculation of the interfacial free energy of liquid-vapour or liquid-liquid phases can be easily attained by analyzing the anisotropy of the stress tensor (Kirkwood & Buff, 1949;Walton et al., 1983).However, the same method leads to flawed results when used in systems involving solid phases (Di Pasquale & Davidchack, 2020).For that reason, here we present two different computational methods that can be directly incorporated in the LAMMPS code to facilitate their functionality.

Scientific background
The Mold Integration method was developed by Espinosa et al. (Espinosa et al., 2014) as a simple and effective approach to calculate the interfacial free energy of crystalline phases in coexistence with a liquid phase.This technique uses potential wells, that are arranged in the simulation box as a mold that replicates a given crystallographic plane of the solid of interest.The attractive wells are modeled through a square-like well potential that is characterized by the well radius   , the well depth   and the well slope .The Mold Integration method consists of the following steps: 1. Preparation of the configuration, by introducing the mold into the liquid phase.2. Simulation of the liquid-mold system at constant well depth and different well radii to determine the limit of well radius (optimal radius) which supplies the exact free energy to form a crystal slab within the liquid.3. Calculation of the reversible work required to form a crystal slab at different well radii above the optimal well radius.4. Linear extrapolation of the interfacial free energy to the optimal well radius.In the top part of Figure 1 we depict a sketch of the Mold Integration method with a generic solid plane and a liquid phase, showing the induction of the solid structure inside the liquid by the interactions of the liquid particles with the potential wells forming the mold (empty black circles).In stage 2, the depth of the wells is fixed to a value such that each well is filled with a fluid particle.Changes in the well radius leads to a modulation of the required free energy to create a crystal slab.A larger radius implies a longer induction time to form the interface so that the optimal radius is defined as the largest radius for which there is no induction time (Espinosa et al., 2014).The thermodynamic integration (stage 3) is performed then for values above the optimal radius that allow reversibility (in the formation of the crystal slab) by varying the well depth for each radius.Extrapolating the interfacial energy to the optimal radius (stage 4) provides the estimation of the interfacial free energy.
The nucleation rate is also a fundamental observable required to characterize nucleation phenomena.In this respect, the idea of using a mold to induce a crystal nucleus in the fluid leads to the Lattice Mold technique (Espinosa, Sampedro, et al., 2016).In this method the nucleation rate is calculated by introducing a mold made of attractive wells whose positions are given by the crystallographic position of the solid at the same temperature for which the nucleation rate aims to be computed.In contrast with the Mold Integration method, the mold is not a crystal plane but a crystal cluster, typically with spherical shape, since this shape minimizes the interfacial area of the emerging cluster within the metastable fluid.In this case, the optimal radius of the well is reached in the limit beyond which more than a single fluid particle occupies a well.Hence, the thermodynamic integration is carried out for smaller radii than the optimal one.The mold must allow the system to nucleate (at maximum well depth) after a sufficiently long induction time, therefore ensuring reversibility across the thermodynamic integration step.This method can be summarized as follows: 1. Preparation of the configuration by introducing the crystal mold nucleus into the liquid phase.2. Simulation with different well radii and a constant value of the well depth (which ensures that wells are filled) in order to determine the optimal well radius.3. Evaluation of the free energy required to fill the mold with liquid particles and generate a pre-critical cluster at several well radii below the optimal one.4. Calculation of the average nucleation time for each well radius over different independent trajectories.5. Extrapolation of the nucleation rate to the optimal well radius.
In the bottom part of Figure 1 we depict a schematic picture of the Lattice Mold technique, showing the different simulation steps to compute the nucleation rate of the system (see Sanchez-Burgos et al. (2022) for further details).The nucleation free energy to form a sub-critical cluster is calculated from the well-occupancy vs well depth curve for each well radius.The average induction time is then calculated by running independent trajectories (at the highest well depth used in the previous step) and measure the time it takes for the pre-critical crystal nucleus to overcome the free energy barrier and percolate the whole system.
To date, the Mold Integration method has been widely used to calculate solid-liquid interfacial free energy of hard spheres (Espinosa et al., 2014), pseudo-hard spheres (Sanchez-Burgos, Sanz, et al., 2021), Lennard-Jones (Espinosa et al., 2014), SPC/E water with JC NaCl (Sanchez-Burgos & Espinosa, 2023), NaCl (Tosi-Fumi) (Espinosa et al., 2015), ice-water interfaces for TIP4P, TPI4P/2005 TIP4P/Ice and mW models (Espinosa, Vega, et al., 2016), and CO 2 hydrate water interfacial energy (Algaba et al., 2022;Romero-Guzmán et al., 2023;Zerón et al., 2022).The Lattice Mold method has been used to calculate the nucleation rate of pseudo-hard spheres and NaCl from its melt (Espinosa, Sampedro, et al., 2016), and mW and TIP4P/Ice water models (Sanchez-Burgos et al., 2022).All these simulations have been carried out with GROMACS v4.6.7 (Abraham et al., 2015) and using a tabulated potential for the square-like attractive interaction between the wells and the fluid particles.Additionally, the typical run is more efficient in terms of computational resources when running on LAMMPS with the analytic form of the potential (Sanchez-Burgos et al., 2022), and the speed increases by a factor of ∼ 1.8.In fact, when using tabulated potentials in GROMACS, long-range corrections in energy and pressure cannot be applied across the calculation in GROMACS v4.6.7, and hence, larger cut-offs in the interaction potential need to be employed (e.g., from 8.5Å to 14Å in water models or electrolyte solutions).This is even more notable when studying complex systems such as atomistic models of water (Sanchez-Burgos et al., 2022) or hydrates (Zerón et al., 2022).Moreover, a tabulated potential needs to be generated for each run with different well parameters, which makes the process of simulation more complex and tedious.This package implements the square-like potential (please see Fig. 3 in (Espinosa et al., 2014) for further discussion on the square-like well potential) in a explicit form on LAMMPS to speed up the simulation performance and simplify the implementation of both methods.
Further to the applications of this method, the same well potential has been used for simulating 'patchy particles' (Russo et al., 2009) that consist of hard spheres with several potential wells on their surface.The present LAMMPS package enables performing efficient simulation of patchy particles using an explicit potential and thus facilitating reproducibility.These works include the study of liquid-liquid phase separation of proteins and other biomolecules (Espinosa et al., 2019(Espinosa et al., , 2020;;Joseph et al., 2021;Sanchez-Burgos, Espinosa, et al., 2021;Sanchez-Burgos, Joseph, et al., 2021;Sanchez-Burgos et al., 2023), as well as crystallization of such particles with different topologies (Garaizar et al., 2022;Romano et al., 2011;Russo et al., 2009;Sciortino & Zaccarelli, 2011).

Functionality
MOLD is a LAMMPS package for the calculation of interfacial free energy and nucleation rates through both the Mold Integration and the Lattice Mold techniques, respectively.The package includes the implementation of the pair potential (with the characteristic syntax of LAMMPS) to simulate the square-like potential that is used in both methods as well as in the simulation of patchy particles.
The code presented here can be found on Github, and it provides a detailed manual with the instructions to compile and use the package as well as with step-by-step examples to implement both the Mold Integration and the Lattice Mold techniques using the explicit square well potential.The package uses the pair_style square/well with three parameters that control the radius, the depth and the slope of the interaction potential.Additionally, we provide two detailed worked examples to reproduce already published results.
The first example explores the Mold Integration technique for the (100) plane of the fcc crystal for the Lennard-Jones potential.This example reproduces the results of the works of (Davidchack & Laird, 2003) using the cleaving technique and (Espinosa et al., 2014) using the Mold Integration method, and uses the Broughton-Gilmer modified Lennard-Jones potential (Broughton & Gilmer, 1983) that is included in this pacakge (please see the CLEAVING package here for further details).The second example shows step-by-step instructions to perform Lattice Mold calculations of mW Ice Ih nucleation rate at T=220K, reproducing the results of (Sanchez-Burgos et al., 2022) using the same method which provides compatible results with other methods (Haji-Akbari et al., 2014;Haji-Akbari, 2018;Li et al., 2011).Both examples include all the LAMMPS scripts, the configuration files and simple bash and Python scripts to run the simulations and analyze the results.
Overall, this simulation package aims to simplify the simulation process of the presented methods and allow the community to perform simulations of liquid-to-solid phase transitions in a more efficient way.
Horizon 2020 research and innovation programme (grant agreement 803326), and from the Ramon y Cajal fellowship.This work has been performed using resources provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service (http://www.hpc.cam.ac.uk) funded by EPSRC Tier-2 capital grant EP/P020259/1 and by the Sulis High-Performance Computational Center under the grant EP/T022108/1.This project made use of time on HPC granted via the UK High-End Computing Consortium for Biomolecular Simulation, HECBioSim (http://hecbiosim.ac.uk), supported by EPSRC (grant no.EP/X035603/1).F. J. B. also acknowledges Ministerio de Ciencia e Innovación (Grant No. PID2021-125081NB-I00), Junta de Andalucía (P20-00363), and Universidad de Huelva (P.O.FEDER UHU-1255522 and FEDER-UHU-202034), all four cofinanced by EU FEDER funds.

Figure 1 :
Figure 1: Top: Sketch of the thermodynamic pathway to compute the interfacial free energy using the Mold Integration method.Bottom: Schematic representation of the different steps in the Lattice Mold technique to compute the nucleation rate.In both schemes the empty black circles represent the potential wells that form the mold and the blue spheres represent the fluid particles.An extrapolation to a given well radius based on calculations performed at different radii is required in both methods.