We examine the ability of two-stage free-energy perturbation methods to yield solid-phase free energies using a system of harmonically coupled particles as a reference. We consider two ways to construct a reference system, one based on derivatives of the intermolecular potential of the target system of interest ͑the conventional choice in lattice dynamics͒, and the other based on analysis of pairwise configurational correlations observed in simulations of the target system. For each case, we

more »

... sider two perturbation techniques that compute the free energy difference between the target and reference systems while avoiding lengthy thermodynamic integration procedures. The methods are overlap sampling as optimized by Bennett, and umbrella sampling optimized in a similar fashion. Such methods require at most two simulations to yield a result, but they can fail if the target and reference do not share a sufficiently large set of relevant configurations. In particular, failure can be expected for large systems, and we examine the question of how large a system can be before this point is reached. Our test case is a face-centered cubic system of r −12 soft spheres, and we find that for systems of up to 108 particles the methods are accurate for all temperatures up to melting; for systems of 256 particles the methods begin to break down at about half the melting temperature. Significantly, we observe that the correction to the harmonic reference is only weakly dependent on system size, suggesting an N-hybrid technique in which perturbation is applied to a small system and the result added to a large-system harmonic reference to obtain a good estimate of the correct large-system free energy. We also examine these approaches, along with thermodynamic integration in temperature, with respect to their computational efficiency. We find that Bennett's method using a derivative-based harmonic reference is the most efficient of all those examined, particularly when employed in the N-hybrid method. Knowledge of the free energy is needed for rigorous evaluation of thermodynamic stability, including that related to phase equilibria. Over the years, many methods have been proposed and developed to compute free energies via molecular simulation. 1 The free energy is rarely computed directly, but rather the relative free energy with respect to a reference system ͑often one with known thermodynamic properties͒ is calculated. One application of particular importance is in the study of crystalline solids. 2,3 The method of choice here usually involves thermodynamic integration, and techniques differ in the formulation of the integration path. Among the earliest approaches is the single-occupancy cell method, 4 in which the solid-state free energy is found by integrating from a low-density state in which noninteracting atoms are each confined in a spherical cell, to the solid state of interest. Probably the most popular method today is due to Frenkel and Ladd, 5 who employed thermodynamic integration ͑TI͒ following a path from the solid of interest to an Einstein crystal, where all the atoms are tethered harmonically to their lattice sites and not interacting with each other. A similar approach can be implemented instead using a reference of coupled harmonic oscillators, with interactions de-termined via a second-order series expansion of the potential of the target system. 3,6 Another recent development is the formulation of lattice-switch methods, which work through a path that reversibly transforms one crystal into another. 7 A somewhat different approach is found in self-referential methods, 8, 9 in which the crystal is grown to some multiple of its original size by reversibly switching on certain degrees of freedom; the extensive nature of the free energy is used to extract an absolute free energy from this free-energy difference. These methodologies have proven their reliability in determining the free energy, but they require significant computational effort and moreover they can be tedious to implement, particularly for molecular crystals. Thus it is a nontrivial task to evaluate solid-phase free energies rigorously, and consequently applications requiring many freeenergy calculations-such as screening of polymorphsrarely make use of them. Instead, crude approximations, such as complete neglect of the entropy, are applied. This leads to uncertainty in identifying the source of error when such methods fail in predictive applications. In this work, we will examine the capability of more direct methods to obtain an absolute free energy. We consider methods that use a system of coupled harmonic oscillators as a reference. The harmonic system includes entropic effects, a͒ Electronic