Efficient, Precise and Accurate Methods of Calculating Solid-Phase Free Energies by Molecular Simulation

Free energy is an important element in evaluating the stability of a solid crystal. A stable form of crystalline solid occupies an atomic/ molecular packing that gives the lowest free energy. The existing reliable free-energy calculation methods usually involve lengthy thermodynamic integration procedures. To improve the calculation efficiency, we present a novel approach that enables us to calculate the absolute free energy of a crystalline solid by applying targeted free-energy perturbation staged using overlap sampling. The underlying idea of this method is based upon the knowledge of phase-space relation between the two perturbing systems and introduction of appropriate scaling parameter (to the atomic displacement and/ or molecular orientation) to improve the overlap of configuration space of the two systems and eventually the free-energy results.    

We have examined this idea to compute the change in free energy with temperature of the soft-sphere crystalline solid. In this method, the free energy difference between nearby temperatures is calculated via overlap-sampling free-energy perturbation with the Bennett’s optimization. Coupled to this is a harmonically targeted perturbation that displaces the atoms in a manner consistent with the temperature change, such that for a harmonic system the free-energy difference would be recovered with no error. A series of such perturbations can be assembled to bridge larger gaps in temperature. An absolute free energy is then computed by implementing the series to near-zero temperature, where the harmonic model becomes very accurate. This method is shown to provide very precise and accurate results. An extension of this method is the absolute free-energy calculation of the realistic linear-molecular nitrogen model, specifically the orientationally-ordered a- and orientationally-disordered β-phase structures. Through this method, the absolute free energies for both the orientationally-ordered and disordered structures are successfully calculated and the coexistence curve for the two phases is traced from zero to 0.12 GPa.