Efficient calculation of temperature dependence of ... - Semantic Scholar
THE JOURNAL OF CHEMICAL PHYSICS 133, 134104 共2010兲
Efficient calculation of temperature dependence of solid-phase free energies by overlap sampling coupled with harmonically targeted perturbation Tai Boon Tan, Andrew J. Schultz, and David A. Kofkea兲 Department of Chemical and Biological Engineering, University at Buffalo, The State University of New York, Buffalo, New York 14260-4200, USA
共Received 5 July 2010; accepted 9 August 2010; published online 4 October 2010兲 We examine a method for computing the change in free energy with temperature of a crystalline solid. In the method, the free-energy difference between nearby temperatures is calculated via overlap-sampling free-energy perturbation with 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. We test this harmonically targeted temperature perturbation 共HTTP兲 method through the application to the inverse-power soft potential, u共r兲 = 共 / r兲n, over a range of temperatures up to the melting condition. Three exponent values 共n = 12, 9, and 6兲 for the potential are studied with different crystal structures, specifically face-centered cubic 共fcc兲, body-centered cubic 共bcc兲, and hexagonal close packing. Absolute free energies 共classical only兲 for each system are obtained by implementing the series to near-zero temperature, where the harmonic model becomes very accurate. The HTTP method is shown to provide very precise results, with errors in the free energy smaller than two parts in 105. An analysis of the thermodynamic stability of the various structures in the infinite-system limit confirms previous findings. In particular, for n = 12 and 9, the fcc structure is stable for all temperatures up to melting, and for n = 6, the bcc crystal becomes stable relative to fcc for temperatures above kT / = 0.802⫾ 0.001. The effects of vacancies and other defects are not considered in the analysis. © 2010 American Institute of Physics. 关doi:10.1063/1.3483899兴 I. BACKGROUND AND INTRODUCTION
Free energy is one of the most fundamental quantities in thermodynamics. Knowledge of free energy is important in understanding the thermodynamic behavior and stability of a system, e.g., the study of phase equilibria and the prediction of crystal stability. Thermodynamic systems are driven to a state where the appropriate free energy is at its minimum with respect to any available degrees of freedom. This forms a foundation to determine the stability of a system at a particular thermodynamic state of interest, for example, in studies of protein folding.1 Very often, free-energy calculations performed in the context of molecular simulation involve great computational expense. Due to the significance of this thermodynamic property, considerable efforts have been devoted to developing new approaches to its calculation and to perfecting the existing methods. Typically, free-energy methods compute a relative rather than the absolute free energy. This is because the absolute free energy is directly related to the phase-space volume that is accessible to the system, and in general, there is no straightforward way to measure this quantity directly and reliably in a simulation. However, there are still exceptional cases where the absolute free energy of a system can be gauged directly, for example, the ideal gas and the harmonia兲
Electronic mail: [email protected].
0021-9606/2010/133共13兲/134104/11/$30.00
cally coupled system, for which there are analytical expression for the free energy. Still, most of the free-energy methods are designed to measure a free-energy difference instead of the absolute quantity itself. These calculation methods can be broadly categorized into two general approaches,2 which are density-of-state methods and work-based methods. Among the latter is thermodynamic integration3 共TI兲, which is a popular method because it is reliable and simple to apply. The efficiency of the method depends on the nature of the integration, and with a smooth, nondivergent integrand, good results can be obtained. Any inaccuracy can be made small by using a fine integration step, but this defeats some of the appeal of the method. Instead, we focus here on another work-based technique, specifically that based on instantaneous work processes and known as free-energy perturbation 共FEP兲.4 The FEP method has been reviewed and studied5–7 in great detail in recent years. The basic FEP working equation to calculate the free-energy difference between two systems 共labeled here system 0 and system 1兲 is given as exp关− 共A1 − A0兲兴 = 具exp关− 共U1 − U0兲兴典0 ,
共1兲
where  is the reciprocal of kT, with k being the Boltzmann constant and T the absolute temperature. Ui is the potential energy of system i. The angle bracket represents the canonical ensemble average, sampled according to the Boltzmann weight of the system specified by the subscript. The FEP
133, 134104-1
© 2010 American Institute of Physics
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-2
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
method is especially convenient to apply, even more so than TI. However, the method is highly prone to bias and inaccuracy, resulting from an inability to fully explore the 1 system when sampling the 0 system 共or vice versa兲.5,6 In this regard, it should be noted that for an unbiased result, it is not sufficient that the important configurations of the 0 and 1 system merely overlap; instead, it is necessary that the configurations of the perturbation or target system 共system 1兲 be a wholly contained subset of the sampling or reference system 共system 0兲, where “contained” is defined with respect to the amount of sampling contributing to the average.7 It is also important that the subset not be a vanishingly small fraction of the sampled configurations. The bias inherent in simple FEP can be severe and difficult to detect, so much so that it should in almost no case be applied in the single-stage form suggested by Eq. 共1兲. A very effective way to improve the accuracy of FEP is to use a two-stage technique, such as overlap or umbrella sampling. In a two-stage perturbation technique, one introduces an intermediate stage and calculates the free-energy difference between each system and the intermediate stage. With such approaches, mere overlap of the important configurations is sufficient to permit an accurate value to be obtained. The intermediate stage can be defined in such a way that the subset relation holds for both perturbation stages, so the overall calculation yields a precise and accurate result. In this work, we apply overlap sampling as optimized by Bennett6,8 to calculate the relative free energy between systems of different temperatures. Another way to improve the overlap of the sample and perturbation systems is to apply targeted perturbation.9 The idea is that an invertible mapping of configurations from one system to the other can be applied such that the important configurations of one can be made to coincide better with those of the other. This is difficult to achieve in general, but we propose here an approach that is highly effective for temperature perturbations of crystals, exploiting the harmonic character of these systems. The method presented here can be applied between any two temperatures to obtain the corresponding free-energy difference. If the solid remains stable to zero temperature 共while treating it entirely classically兲, the perturbation results can yield an absolute free energy for a solid crystalline structure for any temperature that connects reversibly to the lowtemperature system, where the 共analytically solvable兲 harmonic approximation becomes increasingly valid. Absolute free energies are usually needed when comparing phases of different symmetry, for which it is difficult to form a reversible path connecting them. In this work, we demonstrate how this information can be used to gauge the relative stability of different crystal structures. To demonstrate these schemes, we use the classical softsphere potential 共also known as the inverse-power potential兲10,11 as our test case. The interaction energy between two spherical-symmetric atoms given by the pairwiseadditive inverse-power potential is
u共r兲 =
冉冊 r
n
,
where and are the model energy and size parameters, respectively, while r is the atomic separation distance. n is the potential model parameter that dictates the “softness” of the potential. Increasing the value of exponent n makes the potential harder and more short-ranged, and as n → ⬁, the atomic interaction approaches the hard-sphere limit. An appealing feature of this potential is its thermodynamic scaling properties. Dimensionless thermodynamic properties of the inverse-power potential depend on only a single dimensionless state parameter, conventionally defined as ⌫ ⬅ 3共kT / 兲−3/n, where is the atomic density. The n values that are studied for the solid-phase free-energy calculation in this work are 12, 9, and 6. In addition, in the context of studying the relative stability between different crystal packings, we consider for each value of n free energies of face-centered cubic 共fcc兲, body-centered cubic 共bcc兲, and hexagonal close packed 共hcp兲 arrangements. The solid-fluid10–14 and the solid-solid12,13 phase transitions of inverse-power soft-sphere systems have been studied extensively in the past. The phase diagram of this model is well understood and characterized.13,14 For large n values, for example, n = 12 and 9, it is found that the bcc structure is unstable with respect to shear deformation and the fluid phase of such a system freezes into the fcc crystalline structure. However, the bcc phase becomes mechanically stable as n decreases. The structure is also reported to be thermodynamically stable for n less than 7.13 These previous works10–13 provide guidelines in our polymorphic transition study, especially for the n = 6 model. In addition to demonstrating the methods, the present work advances our understanding of these systems by providing free-energy values that are much more precise than previously obtained. The presentation of this paper is arranged such that we first review the formalism and the methods that we employed in this study 共Sec. II兲. Section III present discussion on our free-energy calculations results, and we will then end with concluding remarks for our findings in this work in Sec. IV.
II. FORMALISM AND METHOD
The key idea in this work is the application of targeted perturbation in a way that exploits the approximately harmonic character of the solid phase. This idea is further synthesized with the Bennett-optimized overlap-sampling FEP method5,6,8 to yield a very efficient technique for evaluating the difference in free energy of a crystalline solid at different temperatures. In cases where the solid is stable at very low temperatures, this method can further be used to evaluate absolute free-energy values. We present or review these various elements of the method in the following. It should be noted that in this development and application, we employ an entirely classical framework for the underlying mechanics. We also ignore momentum contributions to the partition function, as these cancel when computing free-energy differences.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-3
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
A. Harmonically targeted temperature perturbation 9
Targeted perturbation was introduced by Jarzynski as a means to improve the performance of FEP calculations by increasing the phase-space overlap of the two systems of interest via an invertible mapping of their coordinates. In applying the idea here, we are guided by the expectation that as the temperature is increased, each atom explores positions displaced further from its lattice site. Specifically, under the assumption that the potential well is quadratic, the magnitude of the displacement is proportional to the square-root of temperature. This suggests that when simulating a system at T0 and perturbing into the system at T1, we will find a configuration more typical of T1 if we accompany the temperature perturbation by a scaling of atomic displacements according to the equation ⌬x1 = ⌬x0
冉 冊 T1 T0
1/2
,
共2兲
where ⌬x is the atom’s Cartesian displacement from its nominal lattice site. Jarzynski points out that in a targeted FEP calculation, the Jacobian of the transformation must be included in the ensemble average. For the scaling described here, the Jacobian J = 兩x1 / x0兩 is independent of atom positions, J=
冉 冊 T1 T0
formulation given by Bennett. The details of overlapsampling method have been well studied and outlined in previous works,5,6,8 so we will provide only a brief review of the main equations here, particularly for this application. For two systems at temperature T0 and T1, respectively, we define an overlap system with weight function as given below,
共3兲
so it may be applied apart from the ensemble-averaging process; here, 3共N − 1兲 appears as the number of degrees of freedom 共not including the center-of-mass兲. Thus, for simple single-stage FEP, we have 共A兲1 = 共A兲0 + ⌬共AC兲 − 23 共N − 1兲ln共T1/T0兲.
Z␥ =
B. Overlap sampling
While this harmonically targeted temperature perturbation 共HTTP兲 method is formulated in the context of a singlestage FEP, it can just as well be applied in a two-stage implementation. In this study, we employ overlap sampling in the
冕
dx0␥共x0兲 =
so Z␥ 1 = Z0 Z0
冕
Z␥ 1 = Z1 Z1
冕
and
dx0
dx1
冕
1 dx1 ␥共x1兲, J
␥共x0兲 e0共x0兲 = e0共x0兲
共6兲
冓冔
1 ␥共x1兲 e1共x1兲 = J e1共x1兲
␥ e0
0
冓 冔 ␥ Je1
. 1
Combining we have
冉 冊
T0 Z0 = Z1 T1
共4兲
The quantity ⌬共AC兲 is the free-energy difference 共1 minus 0兲 between systems T0 and T1, which is measured in the presence of the transformation 关i.e., via Eq. 共1兲 but with the perturbation to U1 accompanied by the transformation in Eq. 共2兲兴, and the term added to it is the application of the Jacobian. This term is, not coincidentally, equal to the difference in free energies of a harmonic system at the two temperatures T1 and T0. In a recent work,15 we examined methods for computing the difference in free energy between a crystal and a harmonic reference system, and there we identified AC ⬅ A − Aharm as the focus of the calculation. We found that this term was much less dependent on system size 共N兲 than either A or Aharm, and we suggested that it could be computed for small systems 共where free-energy methods are easier to apply兲 and the result added to the free energy of a large-N harmonic system to yield a good estimate of the thermodynamic free energy of the target system. It is convenient then that in the proposed method, AC emerges as a direct result of the FEP calculation.
共5兲
where ei = exp共−U共xi兲 / kTi兲 is the Boltzmann factor for the system i. The 3N coordinates x1 are directly given by x0 via the invertible scaling relation in Eq. 共2兲, so ␥ is fully specified by either coordinate set. The parameter ␣ is the optimized Bennett’s parameter, where the optimum value of ␣ is given as ␣ ⬇ exp关−⌬共AC兲兴. The partition function for the overlap system can be expressed in terms of x0 or x1,
3共N−1兲/2
,
e0共x0兲e1共x1兲 , e1共x1兲 + ␣e0共x0兲
␥共x0兲 =
3共N−1兲/2
具␥/e1典1 , 具␥/e0典0
共7兲
so
−1 ⬅ exp关⌬共AC兲兴 =
具␥/e1典1 , 具␥/e0典0
共8兲
again showing, now for the two-stage implementation, that the anharmonic contribution to the free-energy difference is given directly from the targeted-perturbation averages. For notational convenience in the subsequent presentation, we have introduced the symbol defined as in this equation. In the overlap-sampling method, we would normally conduct two simulations to determine the free-energy difference for each adjacent pair of temperatures, one simulation conducted at each temperature and perturbing into their common intermediate region. In the HTTP technique examined in the present work, we are interested in doing this over a series of temperatures, so we can perturb the system simulated at one temperature T into the overlap regions on either side—one for the next lower temperature Tl and one for the next higher temperature, Th. We note that in doing so, the measured values of ␥OS,Tl / eT and ␥OS,Th / eT from the two perturbation averages taken in each direction from the simulation at temperature T are strongly correlated. In order to account for this issue in the error analysis, we measure the covariance and include its contribution to the uncertainty of
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-4
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
the free-energy difference. Other than this, the free energy for two temperatures separated by several overlap-sampling intervals is obtained by summing the differences for the individual overlap-sampling pairs. In each perturbation simulation, multiple initial guesses are made for the ␣ parameters 共11 ␣-values to be exact兲, each of which will yield a different estimate of 关defined in Eq. 共8兲兴. These ␣ values span a range from ␣C · exp共−1兲 to ␣C · exp共1兲, with ␣C being the center value. Bennett8 showed that the variance in i is minimized when i = ␣. The subscript i denotes the perturbation of adjacent systems Ti and Ti+1. We estimate the i value by interpolating ln共i兲 as a function of ln共␣兲 to locate i where ln共i兲 is equal to ln共␣兲. The so-determined i is then used as ␣C for any subsequent simulations having more simulation steps. The temperature interval for the perturbations should be chosen to balance two considerations. The first is that the interval is small enough that the perturbation between the two systems is not so difficult, or in other words, the important phase space of the two systems is not too different, as indicated by i being close to 1.0. The second consideration is that the temperature interval is not so small that spanning a given temperature range requires too many simulations, each of which must be equilibrated for some number of steps.
C. Absolute free energy
The foregoing development describes a procedure to determine the relative free-energy difference between two systems at different temperatures. If the crystalline form of the solid is stable down to zero temperature, we can use the method to evaluate absolute free energies as well. We assume that as we approach T = 0, a harmonic representation of the system becomes increasingly accurate, such that a temperature T sufficiently close to zero can be found to make arbitrarily small the difference in free energy ␦A between the 共classical兲 system of interest and a 共classical兲 harmonic approximation. In practice, ␦A is set by some accuracy goal for the calculation, typically established such that any inaccuracy is less than the standard error in the result. There are two ways to proceed from this point. First, we could imagine performing a temperature perturbation in the manner outlined above, from T to the lowest temperature T0 reached in the perturbation series. Given that the system at this temperature is by definition practically indistinguishable from a harmonic system, one could instead perturb a harmonic system at T to the system of interest at T0. However, if this is done using targeted perturbation as described above, the temperature of the harmonic reference becomes irrelevant. The method is in fact equivalent to a perturbation from a harmonic system to the system of interest, both at T0. This method was described in our recent work,15 and it should be expected to perform satisfactorily for this purpose. An alternative approach is to simply extrapolate the temperature perturbation data to T = 0. For this purpose, we choose a functional form for ⌬AC as an M th-order polyno-
mial in T. For the sake of generality, rather than just temperature, we present the fitting form in terms of the state parameter ⌼ = ⌫−n/3 = 共kT / 兲共3兲−n/3; thus M
ACfit共⌼兲
= 兺 c i⌼ i ,
共9兲
i=1
which enforces the requirement that AC = 0 at kT = 0 共Ref. 11兲 关we also verify that as T → 0, AC → 0, e.g., at kT = 0.001, AC is calculated to be 0.000 193共4兲 共previous work15兲 and 0.000 186 5共3兲 共this work, through fitting兲兴, and we minimize the weighted-least square objective function, fit ⍀ = 兺 共⌬AC,i − ⌬AC,i兲2/␦2i ,
共10兲
i
where ␦i is the uncertainty for ⌬AC,i. With AC expressed in terms of ⌼, the coefficients ci are independent of density and temperature. Results of the fit are presented in the Appendix. The fit may be used simply to establish an absolute free energy for the system at T0 and from this the values at all other temperatures connected to it via the perturbation series, or it may itself provide a surrogate for the free energies over the entire temperature range. Given values for AC at all temperatures of interest, the absolute free energy is obtained by adding the free energy of the harmonic system at each corresponding temperature. This free energy is obtained through the diagonalization of the 3N ⫻ 3N Hessian matrix, as detailed in many places.15,16 The result for the harmonic free energy can be expressed as
Aharm Ulat 1 = + N N 2N
3共N−1兲
兺i
冉 冊
ln
3 ln N 1 i + ln , − 2 N 2kT N 共11兲
where i are the eigenvalues of the Hessian matrix and Ulat is the lattice energy. The last two terms on the right-hand side of the equation are contributions from center-of-mass motion and vanish in the thermodynamic limit N → ⬁. D. Simulation details
We demonstrate the temperature perturbation method by computing free energies for three crystalline structures—fcc, bcc, and hcp—for inverse-power potentials of varying softness, namely, n = 12, 9, and 6. For the ABCABCABC. . .-planes stacking fcc structure, there are four basis atoms in each unit cell, with each basis atom located at scaled positions: 共0, 0, 0兲, 共0, 1/2, 1/2兲, 共1/2, 1/2, 0兲, and 共1/2, 0, 1/2兲, respectively. The bcc structure has two basis atoms at scaled positions: 共0, 0, 0兲 and 共1/2, 1/2, 1/2兲. Both of these structures retain a cubic packing with unit cell length, a given as a = 冑3 k / 3, where k and 3 are the basis number and atomic density respectively. The hcp has two basis atoms at 共0, 0, 0兲 and 共2/3, 1/3, 1/2兲 with the unit cell dimensions defined by vectors: 共a, 0, 0兲, 共−a cos共2 / 3兲, a cos共2 / 3兲, 0兲 and 共0, 0, c兲, where a = 冑2 / 3 is the unit cell length, while c = a冑8 / 3 is the height of the unit cell. Multiple system sizes are examined, specifically N = 500, 864, 1372, 2048, 2916, and 4000 for the fcc crystal 共with four basis atoms兲 and N = 686, 1024, 1458, 2662, 3456, and 4394 for both bcc and
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-5
Solid-phase free energies
hcp crystals, which have two basis atoms in each unit cell. The number particles are selected such that the number of unit cells in each direction are equal 共km3 with m as the number of unit cells in one direction兲. To gauge the effect of the shape of the simulation cell, we also simulated the n = 12 fcc system in a rhombohedral simulation box composed of rhombohedral unit cells, each containing one basis atom. The unit cell length on each edge is a = 冑2 / 3 with angle parameters ␣ =  = ␥ = 2 / 3. The system sizes studied with this simulation box geometry are N = 512, 1000, 1728, 2744, and 4096 particles, respectively. The potential parameters ⑀ and are each taken to be 1 in this study. Each simulation is run in the canonical 共NVT兲 ensemble at densities 3 = 1.1964, 1.3732, and 2.206 17 for n = 12, 9, and 6, respectively. Each density corresponds to melting condition10,11,17 at temperature kT / = 1.0. Each HTTP series starts at a low temperature 共kT / = 0.01兲 and perturbs into a higher temperature with temperature interval of 0.01, proceeding up to temperatures above the respective melting condition, ending at kT / = 1.2. In order to prevent the crystal structures from melting at temperatures above melting, we introduce a constraint potential on each atom, where the atom is always contained in the “enclosure” of its firstneighbor atoms’ shell, and there is no plane sliding allowed. This constraint potential is applied to simulations at all temperatures. We conduct 108 Monte Carlo 共MC兲 trial steps for each HTTP simulation. 50N MC trial moves are proposed as equilibration steps prior to any data collection. In each MC trial move, a pair of atoms selected at random is moved in directions opposite to one another so that the center-of-mass of the system is conserved. To improve the sampling time for systems with softness exponent n = 12 and 9, we apply a potential cutoff of 2.2 and 3.1, respectively. For n = 12, we test the truncation effect on the free-energy calculation for 864-particle system size at melting condition using different truncation radii of 2.2, 2.8, and 3.2, while 3.1, 3.3, and 3.6 are examined as the atomic truncation distance for the n = 9 soft potential. We find that the differences between the free energies calculated using these different atomic truncation radii are negligible, approximately 2 ⫻ 10−5 kT. However, for a longer-range soft potential 共n = 6兲, we apply halfa-simulation-box-length truncation 共=0.48ma兲 to sample the system. Also, in the context of minimizing the finite-size effect, the lattice sum used to calculate the lattice energy of each system and the Hessian for the harmonic free energy are computing using a cutoff of 58 and 15, respectively. The ␣C values are first estimated through a short simulation with 1 ⫻ 106 MC steps before computing the freeenergy difference with a longer simulation run. 108 MC trial steps are proposed in each simulation with the interval for data collection being the number of particles, N. The data are grouped into 100 blocks for error analysis. III. RESULTS AND DISCUSSION A. Results for correction term, AC calculation
In Fig. 1共a兲, we present the i values for the fcc crystal in a cubic box as a function of temperature kT / for different
J. Chem. Phys. 133, 134104 共2010兲
FIG. 1. 共a兲 Plot of i vs kT / for cubic fcc n = 12 soft potential at 3 = 1.1964. The symbols and lines represent the raw data and the fitted equation, respectively. The error bar in 共a兲 is smaller than the symbols. 共b兲 Plot of relative deviation of the simulation results and the fitted data, 共i − i,fitted兲 / i,fitted as a function of temperature, kT / , for N = 500 共black兲 and N = 4000 共violet兲 particles.
system sizes. The value of i is a measure of the free-energy difference between adjacent systems with temperatures Ti and Ti+1, and this does not necessarily relate to the degree of overlap of their important configuration spaces. However, the HTTP scheme developed here is designed to provide perfect overlap when applied to a harmonic system, so it is reasonable to expect that we will have a good degree of configuration-space overlap in the real systems to which the method is applied. In such a circumstance, we might look to i by itself as a rough measure of the difficulty of the HTTP calculation because a lot of the complications will be related to the asymmetry of the sizes of the two configuration spaces, and i can quantify this. A i value close to 1.0 indicates comparably sized configuration spaces, which coupled with an expectation of good overlap for the HTTP method should indicate an easier calculation. Figure 1共a兲 shows that i decreases with system size, indicating that, as expected, the perturbation between temperatures is harder for larger systems. The relative uncertainty, i / i, for each perturbation at any specific temperature 共with temperature inter-
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-6
Tan, Schultz, and Kofke
FIG. 2. AC / N vs kT / for cubic fcc exponent-12 soft potential at 3 = 1.1964. The inset is the blow-up of the curve at melting condition, kT / = 1.0.
val of 0.01兲 is approximately proportional to the system size. For example, for the perturbation between systems at kT / = 0.01 and 0.02, the relative uncertainties for i are about 0.3%, 0.8%, and 3% for N = 500 共smallest system兲, 1372, and 4000 共largest system size兲, respectively. Together with i from raw data, to demonstrate the quality of the fit described by Eq. 共9兲, we also plot i,fitted = exp共 −⌬AC兲 in Fig. 1共a兲, where AC comes from the fit. The simulation data are well described by the polynomial equation as can be seen for 500 and 4000 particles in Fig. 1共b兲. The relative deviation between the simulation data and the fitted line, 共i − i,fitted兲 / i,fitted, fluctuates about zero along the x-axis. The fluctuations are commensurate with the relative uncertainties. The root mean square of the deviation divided by the uncertainty ranges from 0.87 to 1.03. The results for the free-energy correction term, AC / N for different system sizes calculated using Eq. 共9兲, are presented in Fig. 2. We note that the correction term is nearly independent of system size, and in fact the results for different sizes cannot be distinguished on the plot. This behavior is consistent with what we found in our recent work.15 Nevertheless, an expanded scale reveals that AC / N decreases slightly with system size and appears to converge for the large system as shown in the inset of Fig. 2. To examine the system-size effect in more detail, we plot in Fig. 3 the correction term, AC / N, at the melting point as a function of 1 / N. We find the slopes to be 0.10共1兲 and 0.13共2兲 for the cubic and rhombohedral simulation box, respectively 共the last digit in parentheses represents the 67% uncertainty limit兲, which are about ten times smaller than the slopes when plotting the full free energy including the harmonic contribution.15,17 In Fig. 3, we observe that not only the correction term is barely dependent on the size of the system, the uncertainty for this term—as measured using the HTTP method—also does not vary much with N. We find that the values extrapolated to infinite-system size 共1 / N = 0兲 for both simulation box shapes agree well with each other, within their error bars.
J. Chem. Phys. 133, 134104 共2010兲
FIG. 3. n = 12. AC / N vs 1 / N at melting condition kT / = 1.0. The last digit in the parentheses is the 67% confidence limit.
B. Free-energy calculation results
In Fig. 4, we compare our results calculated using the HTTP method with previous work done by Polson et al.17 The plot presents the excess free energy per particle at the melting condition as a function of reciprocal of system size, where the excess free energy includes the harmonic contribution and is given with respect to an ideal gas as Aex / N = 共A − AIG兲 / N; the ideal gas free energy is AIG = ln − 1 + ln共2N兲 / 2N. The plot is constructed adding the term ln共N兲 / N to remove the leading-N dependence.17 We present results using our two box geometries, as well as the results reported by Polson et al., which used a rhombohedral box. Each result is then extrapolated to 1 / N = 0 to obtain its excess free energy for an infinite system. As we can see, the figure shows that the HTTP results for the two box shapes differ for finite-sized systems, but upon extrapolation, they are in excellent mutual agreement. Comparing our results with those of Polson et al., there is a slight disagreement of 0.002 observed. While this is a very small difference, it is significantly larger than 共more
FIG. 4. Excess free energy vs 1 / N at melting condition kT / = 1.0. The last digit in the parentheses and the error bars in the plot represent the 67% confidence limits; the error bars are not easily visible for the present results.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-7
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
TABLE I. Values and estimates of excess free energy 共Aex / N兲N→⬁ of hcp and fcc crystals for soft potential n = 12 at conditions 3 = 1.1964 and kT / = 1.0. The values reported in “This work” are results from hybrid method. Method This work Polson et al.a BenDBb HybridInfb
hcp
fcc
Difference, hcp− fcc
9.231 11共2兲 ¯ 9.235共2兲 9.231共2兲
9.228 31共2兲 9.226 08共7兲 9.231共2兲 9.226共2兲
0.002 80共3兲 0.0028共8兲 0.004共3兲 0.005共3兲
a
Reference 17. Reference 15.
b
than 30⫻兲 the confidence limits of either set of results. We have no definitive explanation for this discrepancy, though our suspicion is that it is due to the finite-step error inherent in the quadrature used for the Frenkel–Ladd method employed by Polson et al.—the difficulty may come, in particular, in handling the integrable singularity that arises in application to soft particles. Barring a mistake in our derivations or calculations and given the good phase-space overlap observed in our overlap-sampling calculations, there is no place in the HTTP method that could lead to an uncharacterized error that might cause this difference. Our preferred approach for calculating the infinitesystem free energy is the “hybrid method” that we presented in our previous study.15 In this approach, the free energy is calculated by performing separate extrapolations of the harmonic and anharmonic contributions to the free energy. The harmonic contribution has a greater system-size dependence, but we can calculate values for much larger systems, thereby improving the extrapolation. In contrast, it is much more difficult to calculate the anharmonic correction AC / N for large systems, but this contribution is much less system-size dependent 共at least for the soft-sphere systems examined so far; see Fig. 3兲. By separating the extrapolations, we can exploit the positive features of each to offset the respective difficulties. The extrapolated excess free energies 共1 / N = 0兲 from this hybrid method 共using HTTP for AC / N兲 are 9.228 32共2兲 and 9.228 30共2兲 for the cubic and rhombohedral simulation boxes, respectively. These results are consistent with the results presented in Fig. 4. We also calculate the free energy for n = 12 hcp solid at corresponding melting condition for the fcc crystal 共3 = 1.1964; kT / = 1.0兲. In similar fashion, we regress the data versus 1 / N as in Fig. 3 to yield the line 0.154 69共2兲 − 0.02共3兲 / N; adding in the extrapolated 共1 / N = 0兲 harmonic contribution yields the value given in Table I. Table I presents the infinite-system results from this work along with the findings15,17 that were previously reported. At melting condition, the fcc crystal is found to be more stable than the hcp crystal by an amount of 0.002 80共3兲. It is interesting to observe that this finding is in perfect agreement with the value reported by Polson et al.,17 even though their individual free energies differ from ours by an amount roughly equal to the hcp-fcc free-energy difference. We note that the results found using the HTTP method are more precise and conclusive than what we reported previously,15 especially in this application of determining the relative stability between different crystalline phases.
FIG. 5. The free-energy difference, 共⌬A / N兲N→⬁ between the competing crystal structures 共fcc, bcc, and hcp for different softness 共n兲 values. The Helmholtz free-energy difference, 共⌬A / N兲N→⬁, e.g., hcp− fcc, denotes 共A / N兲hcp,N→⬁ − 共A / N兲fcc,N→⬁. The subscript N → ⬁ indicates the extrapolated value to infinite-system size. The inset shows the bcc-fcc difference for n = 6 soft potential as a function of temperature. All results are collected at the same density for corresponding softness n values, which are 3 = 1.1964, 1.3732, and 2.20617 for n = 12, 9 and 6 respectively.
C. Relative stability of crystals
In this section, we use the HTTP method to examine the stability of different crystal structures 共fcc, bcc, and hcp兲 for different softness values, n. Lattice dynamics calculations10,11 show that the fcc and hcp crystals are stable 共or metastable兲 for all the softness values that we study here 共n = 12, 9, and 6兲, while the bcc structure is mechanically stable only for n = 6. For each crystalline phase, a constraint potential is imposed during simulation to the prevent the loss of its structure. To test the effect of the constraint potential on the free-energy calculation, we compare AC / N values calculated with and without the constraint at melting condition 共kT / = 1.0兲 using the 864-particle fcc n = 12 system. The effect is found to be negligible, with the calculated values being 0.157 10共3兲 共with兲 and 0.157 06共3兲 共without兲 for the two cases. The results for the crystal-stability calculation using the HTTP method are presented in Fig. 5. Also, see the Appendix for more calculations. Each curve reported in the plot is the infinite-system free-energy difference per particle between two crystals, where the absolute free energy of each crystal structure is determined as described above and the values at kT / = 0.0 indicate the difference in lattice energy between the two structures. We observe that the fcc crystal structure is more stable than the hcp structure for all values of the softness and temperatures, with fcc stability increasing with temperature. Not surprisingly, the free-energy difference between the fcc and hcp is small in comparison to the bcc-fcc difference. The free energy for the n = 6 system favors fcc at low temperatures, and the system undergoes a solid-solid phase transition to bcc at kT / = 0.802, in agreement with previous work, which found a transition temperature of approximately 0.822.13 Also, as can be seen from the inset of Fig. 5, the fcc-bcc difference, ⌬A = ⌬U − T⌬S, slowly decreases with temperature, and this indicates that the bcc
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-8
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
TABLE II. Values and estimates of excess free energy 共Aex / N兲N→⬁ for different crystalline structures 共fcc, bcc, and hcp兲. The excess free energy reported for n = 9 and 6 is at the melting condition 共3 = 1.3732 and kT / = 1.0兲 and the bcc-fcc phase transition 共3 = 2.206 17 and kT / = 0.802兲, respectively. Corrections AC are the results of this work 关Eq. 共9兲兴. n=9 Method Lattice energy Lattice energy+ harmonic Lattice energy+ harmonic+ correction共AC兲
n=6
hcp
fcc
Difference, hcp− fcc
fcc
bcc
Difference, fcc− bcc
6.401 648 46 11.631 837 06 11.689 03共2兲
6.401 293 69 11.627 759 84 11.683 72共2兲
0.000 354 77 0.004 077 22 0.005 31共3兲
22.612 353 77 28.332 440 12 28.314 05共2兲
22.716 955 46 28.303 564 27 28.314 07共4兲
⫺0.104 602 0.028 875 9 ⫺0.000 02共5兲
structure has more thermal fluctuation or higher entropy than the fcc structure over the temperature range. This effect becomes sufficiently important at high temperatures, where the transition occurs. With this HTTP method, we are capable of calculating the free energy for each crystalline phase with high precision and capturing the solid-solid transition with confidence. The largest uncertainty in the Helmholtz freeenergy difference between the crystalline structures shown in Fig. 5 is of the order of 10−5, which is very small on the scale of the plot. To better illustrate the entropic contribution to the relative stability between different crystalline structures, we tabulate in Table II the excess free energy for n = 9 and 6. For n = 9, we determine the relative stability between the hcp and fcc crystal at the fcc melting condition 共3 = 1.3732, kT / = 1.0兲, while the excess free energies of fcc and bcc crystals are calculated for the n = 6 system at the fcc-bcc phase transition 共3 = 2.206 17, kT / = 0.802兲. Table II decomposes each free-energy difference by considering contributions from the lattice energy alone, Ulat, and the lattice energy together with the harmonic contribution 关Eq. 共11兲兴. All the results shown in the table are infinite-limit system size and in excess of the ideal gas contribution, ln共兲 − 1. Corresponding results for n = 12 were presented in previous work.15 Table II shows for the n = 9 system that the lattice energy by itself favors the fcc structure by a very small amount 共0.005%兲. The addition of the entropic contribution to the free energy via lattice dynamics 关Eq. 共11兲兴 amplifies the difference significantly, increasing it by more than a factor of 10. The addition of the correction term further stabilizes the fcc structure but only marginally. The significance of these contributions becomes more prominent when we examine the n = 6 system at the fcc-bcc phase transition. The latticeenergy calculation, a crude approximation with complete neglect of the entropy, shows that the fcc crystal is more stable than the bcc structure by approximately 0.5%. However, the addition of the harmonic contribution switches the stability to bcc by 0.1%. The introduction of the correction term shifts stability back toward fcc, yielding the zero difference. The oscillation in polymorph stability with the introduction of various contributions demonstrates the potential importance of the entropic term in this application. D. Computational efficiency
In this part of the study, we investigate the computational efficiency of HTTP method as a function of temperature interval and number of simulation MC steps for the n = 12 system. The effect of simulation length is gauged by
collecting results using simulations of three lengths 共referred to here as short, medium, and long兲. Temperature intervals examined here are 0.002, 0.01, 0.05, and 0.2. For consistency, we pick the number of MC steps for each temperature interval such that the simulation runs for different temperature intervals employ equal amounts of computational effort by balancing the number of MC steps in each simulation and the total simulations required. Each series starts from the value of the temperature interval, e.g., 0.002, 0.01, 0.05, and 0.2, and proceeds in equal perturbation steps up to temperature, kT / = 1.2 共where again melting is at kT / = 1.0兲. Thus, the series for temperature intervals of 0.002 关0.2, 2, 20兴, 0.01 关1, 10, 100兴, 0.05 关5, 50, 500兴, and 0.2 关20, 200, 2000兴 comprises 600, 120, 24, and 6 simulations, respectively, across the temperature range. The numbers in the square brackets here indicate the total numbers of simulation steps 共in millions兲 considered for the short, medium, and long simulation runs. The calculations were run on an Intel Core2 Quad Q9550 Yorkfield 2.83 GHz processor. The computational efficiency results reported in Fig. 6 are for the 256- and 4000-particle system sizes at kT /
FIG. 6. Error measurement from Eq. 共12兲 vs computation time. The results are presented such that each curve represents results using different temperature intervals, e.g., 256-Int0.002 indicates the error measurement for 256particle system size with temperature interval of 0.002. The filled circle and triangle are for 256- and 4000-particle systems, respectively. Lines join points that apply the same method but use increasing 共from left to right兲 amounts of sampling at each temperature. TI and BenDB indicate the thermodynamic integration and harmonic perturbation methods 共Ref. 15兲, respectively. Each point within the TI curve differs by the number of MC steps at each temperature 共from left to right兲: 62.5, 125, 250, and 500⫻ 106. The solid and dashed lines are aids to the eye. The solid line with “1/2” indicates slope of magnitude 0.5.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-9
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
= 1.0. The accuracy of the HTTP method is calculated by comparing the correction term, AC, determined by the method under consideration, with what we take as the true value determined using temperature interval of 0.01 with 1 ⫻ 109 MC trial steps for each simulation. The error measure, ⌺, inclusive of accuracy and precision is given by
⌺=
冑 共AC − AC,true兲2 + 共SE兲2 , N
共12兲
where SE is the standard error of the mean in the final result. In Fig. 6, the three data points that are connected by the solid 共N = 256兲 or dashed lines 共N = 4000兲 are the error measurements determined from the short 共leftmost point兲, medium 共center兲, and long 共rightmost point兲 simulation runs, respectively, for different temperature intervals. As we can observe from the plot, the short, medium, and long runs 共of the same temperature interval兲 for both system sizes correspondingly amount to almost equal total computation hours of approximately 0.8, 8, and 80 h, respectively. However, a slight variation of about 7 h is observed in the total computation time between the 256- and 4000-particle system sizes for the long simulation run. We believe that the computation expense of the latter system has been greatly improved by the use of a potential cutoff of 2.2, for which the atomic interactions beyond the truncation radius are neglected. Figure 6 shows that 0.002, 0.01, and 0.05 are equally effective temperature intervals to determine the correction term, AC, using the HTTP method. The error measurement, ⌺, of these intervals decreases with the increasing number of MC steps, with a slope of 1/2. This suggests that the simulation results become more precise and converge to the true value 共more accurate兲, and the error for the long simulation comes primarily from the imprecision of the results. In contrast, for the temperature interval of 0.2, the computation efficiency does not improve with increasing calculation efforts. The error is mainly dominated by the inaccuracy, which records differences of 0.0007 and 0.004 from the true value, for N = 256 and 4000, respectively, for the long simulation run. The inaccuracies of these results outweigh their standard error by an order of magnitude, and the large error of this interval is well illustrated in Fig. 6. This large error is mainly caused by the poor perturbation between the two systems at different temperatures, where its temperature interval is apparently too large to have a good overlap in phase space between two systems. The problem worsens with increasing system size. Another reason is the extrapolation of the data to zero from temperature at kT / = 0.2, which seems to be a less appropriate procedure to follow in this case. One way to overcome the latter problem is to introduce an extra point at kT / = 0.05 in order to have a better regression of the data to zero. Judging from the plot, the temperature intervals of 0.01 and 0.05 are reasonable choices for this study because both intervals have a balance between 共1兲 a good amount of overlap in phase space between the perturbing systems and 共2兲 the number of simulations required. Also in the plot, we present the results calculated using
TI in temperature and the harmonic perturbation method 共BenDB兲, both of which were demonstrated in a previous work.15 The error measures at the melting condition from both methods are approximately two orders of magnitude larger than the one calculated using the HTTP method. In the BenDB method, for the 108-particle system size, the error is mainly contributed by its uncertainty, which is 0.0017. The result for the 256-particle system suffers from both its inaccuracy and imprecision; this is where we believe that the limit of the BenDB method is reached. TI shows an error of the order of 10−3. It improves slightly with increased computation, reflecting that it is prone to error introduced by both the quadrature method as well as stochastic errors in the integrated data. It should be noted that the BenDB method provides the free energy at only one point, in contrast to the TI and HTTP methods, which provide the free energy over a range of temperatures as a by-product of the calculation examined in Fig. 6. On the other hand, the TI and HTTP methods are unable provide absolute free energies without taking this path, while the BenDB method yields an absolute free energy without requiring the possibility of connecting to zero temperature. To demonstrate the importance of harmonic targeting as part of the overall HTTP method, we also employed the Bennett-optimized overlap-sampling FEP method to calculate AC / N between systems at different temperatures without performing scaling of the atomic positions. For both 500and 4000-particle system sizes with a temperature step of 0.01, we find that relative uncertainty of i 共or equivalently, the absolute error in the free energy per particle兲 is about an order of magnitude larger than that found using the HTTP method applied to the same pair of temperatures. This is true for perturbation steps over almost all of the temperature range and both system sizes. Further, the nontargeted overlap-sampling perturbation fails completely—with relative uncertainty of i being more than 70%—in the lowtemperature region 共e.g., perturbing from temperatures of 0.01–0.02兲, where the ratio of the phase spaces of the two perturbing systems is large. In addition to the poor precision, we also observe that the low-temperature perturbation results 共without harmonic targeting兲 are prone to inaccuracy. The comparison between targeted and nontargeted temperature perturbation is likely to be only worse for temperature steps larger than the 0.01 increment examined here. IV. CONCLUSIONS
We have demonstrated the effectiveness of harmonically targeted temperature perturbation, coupled with Bennettoptimized overlap sampling, to calculate the free-energy difference between systems at different temperatures, and further to compute the 共classical兲 absolute free energy by connecting to a harmonic system at zero temperature. The demonstration here is for a simple soft-sphere model, considering three values of the softness parameter. The method gives directly the free energy in excess of a harmonic system at the same temperature, and we show again that this quantity 共AC兲 is much less dependent on system size than the harmonic contribution is to the free energy and thus the total
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-10
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
TABLE III. Values for the lattice energy and vibrational contribution to the free energy.
n = 12
fcc hcp fcc hcp fcc bcc hcp
n=9 n=6
共Ulat / N兲N→⬁
共Avib,0 / N兲N→⬁
1.516 485 025 1.516 536 721 2.208 391 122 2.208 528 128 3.613 475 382 3.630 711 328 3.613 718 378
3.888 452 45 3.891 589 92 3.482 175 24 3.485 900 58 2.698 198 34 2.564 731 10 2.702 699 96
free energy itself. This HTTP method is capable of providing results of very high precision. The good precision obtained here allows us to determine the relative stability between different crystalline structures with high confidence. The results obtained here neglect the effects of thermodynamically stable defects, vacancies in particular, which have contributions greater than the errors in our calculated free energies.15 The effects of vacancy defects are however negligible compared to the harmonic and anharmonic contributions 共which are at least an order of magnitude larger兲. Also, a previous study18 has demonstrated that the relative stability of these crystal structures is not affected by the presence of defects. For soft potential with n = 12 and 9, we show that the fcc crystal is more stable than the hcp structure at all temperatures up to melting. The hcp-fcc free-energy difference increases substantially with temperature, indicating the significance of the entropic effects in determining the relative stability of crystal structures at nonzero temperature condition. This conclusion is further supported by the emergence of the bcc crystal structure at kT / = 0.802 for the n = 6 soft potential, where the fcc structure gives way to bcc crystal. The focus of this work is on the soft-sphere model, first to demonstrate the method, and second so that we can establish highly precise values for its free energy. Extension to more complex multiatomics is an important next step in demonstrating the utility of this approach.
ACKNOWLEDGMENTS
Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for support of this research. Computational resources were provided by the University at Buffalo Center for Computational Research.
APPENDIX: FREE ENERGY FORMULAS
For completeness, we present the expressions that we have determined for the free energy of the defect-free softsphere system of different exponent n values: 12, 9, and 6 with various crystal packings 共fcc, bcc, and hcp兲. In the thermodynamic limit, the configurational free energy 共apart from momentum contributions兲 can be expressed as
A Ulat Avib AC = + + . N N N N
共A1兲
The static-lattice-energy contribution is given in terms of density and temperature according to
冉 冊
Ulat 1 Ulat,0 = N ⌼ N
共A2兲
, N→⬁
where again ⌼ = 共kT / 兲共3兲−n/3 is the single independent state parameter for the model. The term in parentheses is the infinite-system lattice sum energy at unit density and is tabulated for the various cases in Table III. The vibrational contribution in Eq. 共A1兲 is the limiting value of the sum in Eq. 共11兲, and it depends on temperature and density according to
冉
Avib Avib,0 = N N
冊
3 − ln ⌼ + ln共3兲. 2 N→⬁
共A3兲
The temperature- and density-independent term Avib,0 / N also is presented in Table III. Finally, the anharmonic contribution Ac is given as a function of temperature and density via the group ⌼ in Eq. 共9兲. The state-independent constants ci are presented in Table IV. The infinite-system configurational free energy can be computed using these formulas for the corresponding crystal at any temperature and density up to the melting point 共values of ⌼ below about 0.5, 0.4, and 0.2 for n = 12, 9, and 6, respectively兲, with the coefficients given in Tables III and IV.
TABLE IV. Fitted coefficients for 共AC / N兲N→⬁ as defined in Eq. 共9兲.
n = 12 n=9 n=6
fcc hcp fcc hcp fcc bcc hcp
c1
c2
c3
c4
c5
c6
c7
0.464 148 0.462 065 0.249 526 0.250 099 ⫺0.024 270 0.283 913 ⫺0.023 687
⫺0.423 821 ⫺0.413 562 ⫺0.245 724 ⫺0.282 838 ⫺0.385 254 ⫺2.421 458 ⫺0.351 963
0.466 322 0.429 175 ⫺0.500 979 0.007 132 ⫺1.989 540 14.558 777 ⫺2.934 478
⫺0.526 094 ⫺0.461 716 4.258 323 0.922 672 24.366 184 ⫺77.105 689 37.070 278
0.148 125 0.111 491 ⫺16.542 027 ⫺5.218 832 ⫺208.856 577 229.711 416 ⫺280.028 273
0.401 893 0.368 595 30.811 960 10.890 988 860.698 810 ⫺391.161 702 953.287 096
⫺0.384 095 ⫺0.352 244 ⫺23.208 566 ⫺7.3695 08 ⫺1469.431 853 335.089 915 ⫺1140.820 226
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-11 1
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
P. A. Bash, U. C. Singh, F. K. Brown, R. Langridge, and P. A. Kollman, Science 235, 574 共1987兲; L. X. Dang, K. M. Merz, and P. A. Kollman, J. Am. Chem. Soc. 111, 8505 共1989兲. 2 D. A. Kofke, Fluid Phase Equilib. 228–229, 41 共2005兲. 3 D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. 共Academic, San Diego, 2002兲. 4 R. W. Zwanzig, J. Chem. Phys. 22, 1420 共1954兲. 5 N. D. Lu, D. A. Kofke, and T. B. Woolf, J. Comput. Chem. 25, 28 共2004兲. 6 N. D. Lu, J. K. Singh, and D. A. Kofke, J. Chem. Phys. 118, 2977 共2003兲; N. D. Lu, D. Wu, T. B. Woolf, and D. A. Kofke, Phys. Rev. E 69, 057702 共2004兲. 7 D. Wu and D. A. Kofke, J. Chem. Phys. 123, 054103 共2005兲; 123, 084109 共2005兲. 8 C. H. Bennett, J. Comput. Phys. 22, 245 共1976兲. 9 C. Jarzynski, Phys. Rev. E 65, 046122 共2002兲.
10
W. G. Hoover, S. G. Gray, and K. W. Johnson, J. Chem. Phys. 55, 1128 共1971兲. 11 W. G. Hoover, M. Ross, K. W. Johnson, D. Henderson, J. A. Barker, and B. C. Brown, J. Chem. Phys. 52, 4931 共1970兲. 12 R. Agrawal and D. A. Kofke, Mol. Phys. 85, 23 共1995兲. 13 S. Prestipino, F. Saija, and P. V. Giaquinta, J. Chem. Phys. 123, 144110 共2005兲. 14 D. C. Wang and A. P. Gast, J. Phys.: Condens. Matter 11, 10133 共1999兲. 15 T. B. Tan, A. J. Schultz, and D. A. Kofke, J. Chem. Phys. 132, 214103 共2010兲. 16 M. T. Dove, Introduction to Lattice Dynamics 共Cambridge University Press, New York, 2005兲. 17 J. M. Polson, E. Trizac, S. Pronk, and D. Frenkel, J. Chem. Phys. 112, 5339 共2000兲. 18 S. K. Kwak and D. A. Kofke, J. Chem. Phys. 122, 176101 共2005兲.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
Efficient calculation of temperature dependence of solid-phase free energies by overlap sampling coupled with harmonically targeted perturbation Tai Boon Tan, Andrew J. Schultz, and David A. Kofkea兲 Department of Chemical and Biological Engineering, University at Buffalo, The State University of New York, Buffalo, New York 14260-4200, USA
共Received 5 July 2010; accepted 9 August 2010; published online 4 October 2010兲 We examine a method for computing the change in free energy with temperature of a crystalline solid. In the method, the free-energy difference between nearby temperatures is calculated via overlap-sampling free-energy perturbation with 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. We test this harmonically targeted temperature perturbation 共HTTP兲 method through the application to the inverse-power soft potential, u共r兲 = 共 / r兲n, over a range of temperatures up to the melting condition. Three exponent values 共n = 12, 9, and 6兲 for the potential are studied with different crystal structures, specifically face-centered cubic 共fcc兲, body-centered cubic 共bcc兲, and hexagonal close packing. Absolute free energies 共classical only兲 for each system are obtained by implementing the series to near-zero temperature, where the harmonic model becomes very accurate. The HTTP method is shown to provide very precise results, with errors in the free energy smaller than two parts in 105. An analysis of the thermodynamic stability of the various structures in the infinite-system limit confirms previous findings. In particular, for n = 12 and 9, the fcc structure is stable for all temperatures up to melting, and for n = 6, the bcc crystal becomes stable relative to fcc for temperatures above kT / = 0.802⫾ 0.001. The effects of vacancies and other defects are not considered in the analysis. © 2010 American Institute of Physics. 关doi:10.1063/1.3483899兴 I. BACKGROUND AND INTRODUCTION
Free energy is one of the most fundamental quantities in thermodynamics. Knowledge of free energy is important in understanding the thermodynamic behavior and stability of a system, e.g., the study of phase equilibria and the prediction of crystal stability. Thermodynamic systems are driven to a state where the appropriate free energy is at its minimum with respect to any available degrees of freedom. This forms a foundation to determine the stability of a system at a particular thermodynamic state of interest, for example, in studies of protein folding.1 Very often, free-energy calculations performed in the context of molecular simulation involve great computational expense. Due to the significance of this thermodynamic property, considerable efforts have been devoted to developing new approaches to its calculation and to perfecting the existing methods. Typically, free-energy methods compute a relative rather than the absolute free energy. This is because the absolute free energy is directly related to the phase-space volume that is accessible to the system, and in general, there is no straightforward way to measure this quantity directly and reliably in a simulation. However, there are still exceptional cases where the absolute free energy of a system can be gauged directly, for example, the ideal gas and the harmonia兲
Electronic mail: [email protected].
0021-9606/2010/133共13兲/134104/11/$30.00
cally coupled system, for which there are analytical expression for the free energy. Still, most of the free-energy methods are designed to measure a free-energy difference instead of the absolute quantity itself. These calculation methods can be broadly categorized into two general approaches,2 which are density-of-state methods and work-based methods. Among the latter is thermodynamic integration3 共TI兲, which is a popular method because it is reliable and simple to apply. The efficiency of the method depends on the nature of the integration, and with a smooth, nondivergent integrand, good results can be obtained. Any inaccuracy can be made small by using a fine integration step, but this defeats some of the appeal of the method. Instead, we focus here on another work-based technique, specifically that based on instantaneous work processes and known as free-energy perturbation 共FEP兲.4 The FEP method has been reviewed and studied5–7 in great detail in recent years. The basic FEP working equation to calculate the free-energy difference between two systems 共labeled here system 0 and system 1兲 is given as exp关− 共A1 − A0兲兴 = 具exp关− 共U1 − U0兲兴典0 ,
共1兲
where  is the reciprocal of kT, with k being the Boltzmann constant and T the absolute temperature. Ui is the potential energy of system i. The angle bracket represents the canonical ensemble average, sampled according to the Boltzmann weight of the system specified by the subscript. The FEP
133, 134104-1
© 2010 American Institute of Physics
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-2
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
method is especially convenient to apply, even more so than TI. However, the method is highly prone to bias and inaccuracy, resulting from an inability to fully explore the 1 system when sampling the 0 system 共or vice versa兲.5,6 In this regard, it should be noted that for an unbiased result, it is not sufficient that the important configurations of the 0 and 1 system merely overlap; instead, it is necessary that the configurations of the perturbation or target system 共system 1兲 be a wholly contained subset of the sampling or reference system 共system 0兲, where “contained” is defined with respect to the amount of sampling contributing to the average.7 It is also important that the subset not be a vanishingly small fraction of the sampled configurations. The bias inherent in simple FEP can be severe and difficult to detect, so much so that it should in almost no case be applied in the single-stage form suggested by Eq. 共1兲. A very effective way to improve the accuracy of FEP is to use a two-stage technique, such as overlap or umbrella sampling. In a two-stage perturbation technique, one introduces an intermediate stage and calculates the free-energy difference between each system and the intermediate stage. With such approaches, mere overlap of the important configurations is sufficient to permit an accurate value to be obtained. The intermediate stage can be defined in such a way that the subset relation holds for both perturbation stages, so the overall calculation yields a precise and accurate result. In this work, we apply overlap sampling as optimized by Bennett6,8 to calculate the relative free energy between systems of different temperatures. Another way to improve the overlap of the sample and perturbation systems is to apply targeted perturbation.9 The idea is that an invertible mapping of configurations from one system to the other can be applied such that the important configurations of one can be made to coincide better with those of the other. This is difficult to achieve in general, but we propose here an approach that is highly effective for temperature perturbations of crystals, exploiting the harmonic character of these systems. The method presented here can be applied between any two temperatures to obtain the corresponding free-energy difference. If the solid remains stable to zero temperature 共while treating it entirely classically兲, the perturbation results can yield an absolute free energy for a solid crystalline structure for any temperature that connects reversibly to the lowtemperature system, where the 共analytically solvable兲 harmonic approximation becomes increasingly valid. Absolute free energies are usually needed when comparing phases of different symmetry, for which it is difficult to form a reversible path connecting them. In this work, we demonstrate how this information can be used to gauge the relative stability of different crystal structures. To demonstrate these schemes, we use the classical softsphere potential 共also known as the inverse-power potential兲10,11 as our test case. The interaction energy between two spherical-symmetric atoms given by the pairwiseadditive inverse-power potential is
u共r兲 =
冉冊 r
n
,
where and are the model energy and size parameters, respectively, while r is the atomic separation distance. n is the potential model parameter that dictates the “softness” of the potential. Increasing the value of exponent n makes the potential harder and more short-ranged, and as n → ⬁, the atomic interaction approaches the hard-sphere limit. An appealing feature of this potential is its thermodynamic scaling properties. Dimensionless thermodynamic properties of the inverse-power potential depend on only a single dimensionless state parameter, conventionally defined as ⌫ ⬅ 3共kT / 兲−3/n, where is the atomic density. The n values that are studied for the solid-phase free-energy calculation in this work are 12, 9, and 6. In addition, in the context of studying the relative stability between different crystal packings, we consider for each value of n free energies of face-centered cubic 共fcc兲, body-centered cubic 共bcc兲, and hexagonal close packed 共hcp兲 arrangements. The solid-fluid10–14 and the solid-solid12,13 phase transitions of inverse-power soft-sphere systems have been studied extensively in the past. The phase diagram of this model is well understood and characterized.13,14 For large n values, for example, n = 12 and 9, it is found that the bcc structure is unstable with respect to shear deformation and the fluid phase of such a system freezes into the fcc crystalline structure. However, the bcc phase becomes mechanically stable as n decreases. The structure is also reported to be thermodynamically stable for n less than 7.13 These previous works10–13 provide guidelines in our polymorphic transition study, especially for the n = 6 model. In addition to demonstrating the methods, the present work advances our understanding of these systems by providing free-energy values that are much more precise than previously obtained. The presentation of this paper is arranged such that we first review the formalism and the methods that we employed in this study 共Sec. II兲. Section III present discussion on our free-energy calculations results, and we will then end with concluding remarks for our findings in this work in Sec. IV.
II. FORMALISM AND METHOD
The key idea in this work is the application of targeted perturbation in a way that exploits the approximately harmonic character of the solid phase. This idea is further synthesized with the Bennett-optimized overlap-sampling FEP method5,6,8 to yield a very efficient technique for evaluating the difference in free energy of a crystalline solid at different temperatures. In cases where the solid is stable at very low temperatures, this method can further be used to evaluate absolute free-energy values. We present or review these various elements of the method in the following. It should be noted that in this development and application, we employ an entirely classical framework for the underlying mechanics. We also ignore momentum contributions to the partition function, as these cancel when computing free-energy differences.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-3
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
A. Harmonically targeted temperature perturbation 9
Targeted perturbation was introduced by Jarzynski as a means to improve the performance of FEP calculations by increasing the phase-space overlap of the two systems of interest via an invertible mapping of their coordinates. In applying the idea here, we are guided by the expectation that as the temperature is increased, each atom explores positions displaced further from its lattice site. Specifically, under the assumption that the potential well is quadratic, the magnitude of the displacement is proportional to the square-root of temperature. This suggests that when simulating a system at T0 and perturbing into the system at T1, we will find a configuration more typical of T1 if we accompany the temperature perturbation by a scaling of atomic displacements according to the equation ⌬x1 = ⌬x0
冉 冊 T1 T0
1/2
,
共2兲
where ⌬x is the atom’s Cartesian displacement from its nominal lattice site. Jarzynski points out that in a targeted FEP calculation, the Jacobian of the transformation must be included in the ensemble average. For the scaling described here, the Jacobian J = 兩x1 / x0兩 is independent of atom positions, J=
冉 冊 T1 T0
formulation given by Bennett. The details of overlapsampling method have been well studied and outlined in previous works,5,6,8 so we will provide only a brief review of the main equations here, particularly for this application. For two systems at temperature T0 and T1, respectively, we define an overlap system with weight function as given below,
共3兲
so it may be applied apart from the ensemble-averaging process; here, 3共N − 1兲 appears as the number of degrees of freedom 共not including the center-of-mass兲. Thus, for simple single-stage FEP, we have 共A兲1 = 共A兲0 + ⌬共AC兲 − 23 共N − 1兲ln共T1/T0兲.
Z␥ =
B. Overlap sampling
While this harmonically targeted temperature perturbation 共HTTP兲 method is formulated in the context of a singlestage FEP, it can just as well be applied in a two-stage implementation. In this study, we employ overlap sampling in the
冕
dx0␥共x0兲 =
so Z␥ 1 = Z0 Z0
冕
Z␥ 1 = Z1 Z1
冕
and
dx0
dx1
冕
1 dx1 ␥共x1兲, J
␥共x0兲 e0共x0兲 = e0共x0兲
共6兲
冓冔
1 ␥共x1兲 e1共x1兲 = J e1共x1兲
␥ e0
0
冓 冔 ␥ Je1
. 1
Combining we have
冉 冊
T0 Z0 = Z1 T1
共4兲
The quantity ⌬共AC兲 is the free-energy difference 共1 minus 0兲 between systems T0 and T1, which is measured in the presence of the transformation 关i.e., via Eq. 共1兲 but with the perturbation to U1 accompanied by the transformation in Eq. 共2兲兴, and the term added to it is the application of the Jacobian. This term is, not coincidentally, equal to the difference in free energies of a harmonic system at the two temperatures T1 and T0. In a recent work,15 we examined methods for computing the difference in free energy between a crystal and a harmonic reference system, and there we identified AC ⬅ A − Aharm as the focus of the calculation. We found that this term was much less dependent on system size 共N兲 than either A or Aharm, and we suggested that it could be computed for small systems 共where free-energy methods are easier to apply兲 and the result added to the free energy of a large-N harmonic system to yield a good estimate of the thermodynamic free energy of the target system. It is convenient then that in the proposed method, AC emerges as a direct result of the FEP calculation.
共5兲
where ei = exp共−U共xi兲 / kTi兲 is the Boltzmann factor for the system i. The 3N coordinates x1 are directly given by x0 via the invertible scaling relation in Eq. 共2兲, so ␥ is fully specified by either coordinate set. The parameter ␣ is the optimized Bennett’s parameter, where the optimum value of ␣ is given as ␣ ⬇ exp关−⌬共AC兲兴. The partition function for the overlap system can be expressed in terms of x0 or x1,
3共N−1兲/2
,
e0共x0兲e1共x1兲 , e1共x1兲 + ␣e0共x0兲
␥共x0兲 =
3共N−1兲/2
具␥/e1典1 , 具␥/e0典0
共7兲
so
−1 ⬅ exp关⌬共AC兲兴 =
具␥/e1典1 , 具␥/e0典0
共8兲
again showing, now for the two-stage implementation, that the anharmonic contribution to the free-energy difference is given directly from the targeted-perturbation averages. For notational convenience in the subsequent presentation, we have introduced the symbol defined as in this equation. In the overlap-sampling method, we would normally conduct two simulations to determine the free-energy difference for each adjacent pair of temperatures, one simulation conducted at each temperature and perturbing into their common intermediate region. In the HTTP technique examined in the present work, we are interested in doing this over a series of temperatures, so we can perturb the system simulated at one temperature T into the overlap regions on either side—one for the next lower temperature Tl and one for the next higher temperature, Th. We note that in doing so, the measured values of ␥OS,Tl / eT and ␥OS,Th / eT from the two perturbation averages taken in each direction from the simulation at temperature T are strongly correlated. In order to account for this issue in the error analysis, we measure the covariance and include its contribution to the uncertainty of
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-4
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
the free-energy difference. Other than this, the free energy for two temperatures separated by several overlap-sampling intervals is obtained by summing the differences for the individual overlap-sampling pairs. In each perturbation simulation, multiple initial guesses are made for the ␣ parameters 共11 ␣-values to be exact兲, each of which will yield a different estimate of 关defined in Eq. 共8兲兴. These ␣ values span a range from ␣C · exp共−1兲 to ␣C · exp共1兲, with ␣C being the center value. Bennett8 showed that the variance in i is minimized when i = ␣. The subscript i denotes the perturbation of adjacent systems Ti and Ti+1. We estimate the i value by interpolating ln共i兲 as a function of ln共␣兲 to locate i where ln共i兲 is equal to ln共␣兲. The so-determined i is then used as ␣C for any subsequent simulations having more simulation steps. The temperature interval for the perturbations should be chosen to balance two considerations. The first is that the interval is small enough that the perturbation between the two systems is not so difficult, or in other words, the important phase space of the two systems is not too different, as indicated by i being close to 1.0. The second consideration is that the temperature interval is not so small that spanning a given temperature range requires too many simulations, each of which must be equilibrated for some number of steps.
C. Absolute free energy
The foregoing development describes a procedure to determine the relative free-energy difference between two systems at different temperatures. If the crystalline form of the solid is stable down to zero temperature, we can use the method to evaluate absolute free energies as well. We assume that as we approach T = 0, a harmonic representation of the system becomes increasingly accurate, such that a temperature T sufficiently close to zero can be found to make arbitrarily small the difference in free energy ␦A between the 共classical兲 system of interest and a 共classical兲 harmonic approximation. In practice, ␦A is set by some accuracy goal for the calculation, typically established such that any inaccuracy is less than the standard error in the result. There are two ways to proceed from this point. First, we could imagine performing a temperature perturbation in the manner outlined above, from T to the lowest temperature T0 reached in the perturbation series. Given that the system at this temperature is by definition practically indistinguishable from a harmonic system, one could instead perturb a harmonic system at T to the system of interest at T0. However, if this is done using targeted perturbation as described above, the temperature of the harmonic reference becomes irrelevant. The method is in fact equivalent to a perturbation from a harmonic system to the system of interest, both at T0. This method was described in our recent work,15 and it should be expected to perform satisfactorily for this purpose. An alternative approach is to simply extrapolate the temperature perturbation data to T = 0. For this purpose, we choose a functional form for ⌬AC as an M th-order polyno-
mial in T. For the sake of generality, rather than just temperature, we present the fitting form in terms of the state parameter ⌼ = ⌫−n/3 = 共kT / 兲共3兲−n/3; thus M
ACfit共⌼兲
= 兺 c i⌼ i ,
共9兲
i=1
which enforces the requirement that AC = 0 at kT = 0 共Ref. 11兲 关we also verify that as T → 0, AC → 0, e.g., at kT = 0.001, AC is calculated to be 0.000 193共4兲 共previous work15兲 and 0.000 186 5共3兲 共this work, through fitting兲兴, and we minimize the weighted-least square objective function, fit ⍀ = 兺 共⌬AC,i − ⌬AC,i兲2/␦2i ,
共10兲
i
where ␦i is the uncertainty for ⌬AC,i. With AC expressed in terms of ⌼, the coefficients ci are independent of density and temperature. Results of the fit are presented in the Appendix. The fit may be used simply to establish an absolute free energy for the system at T0 and from this the values at all other temperatures connected to it via the perturbation series, or it may itself provide a surrogate for the free energies over the entire temperature range. Given values for AC at all temperatures of interest, the absolute free energy is obtained by adding the free energy of the harmonic system at each corresponding temperature. This free energy is obtained through the diagonalization of the 3N ⫻ 3N Hessian matrix, as detailed in many places.15,16 The result for the harmonic free energy can be expressed as
Aharm Ulat 1 = + N N 2N
3共N−1兲
兺i
冉 冊
ln
3 ln N 1 i + ln , − 2 N 2kT N 共11兲
where i are the eigenvalues of the Hessian matrix and Ulat is the lattice energy. The last two terms on the right-hand side of the equation are contributions from center-of-mass motion and vanish in the thermodynamic limit N → ⬁. D. Simulation details
We demonstrate the temperature perturbation method by computing free energies for three crystalline structures—fcc, bcc, and hcp—for inverse-power potentials of varying softness, namely, n = 12, 9, and 6. For the ABCABCABC. . .-planes stacking fcc structure, there are four basis atoms in each unit cell, with each basis atom located at scaled positions: 共0, 0, 0兲, 共0, 1/2, 1/2兲, 共1/2, 1/2, 0兲, and 共1/2, 0, 1/2兲, respectively. The bcc structure has two basis atoms at scaled positions: 共0, 0, 0兲 and 共1/2, 1/2, 1/2兲. Both of these structures retain a cubic packing with unit cell length, a given as a = 冑3 k / 3, where k and 3 are the basis number and atomic density respectively. The hcp has two basis atoms at 共0, 0, 0兲 and 共2/3, 1/3, 1/2兲 with the unit cell dimensions defined by vectors: 共a, 0, 0兲, 共−a cos共2 / 3兲, a cos共2 / 3兲, 0兲 and 共0, 0, c兲, where a = 冑2 / 3 is the unit cell length, while c = a冑8 / 3 is the height of the unit cell. Multiple system sizes are examined, specifically N = 500, 864, 1372, 2048, 2916, and 4000 for the fcc crystal 共with four basis atoms兲 and N = 686, 1024, 1458, 2662, 3456, and 4394 for both bcc and
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-5
Solid-phase free energies
hcp crystals, which have two basis atoms in each unit cell. The number particles are selected such that the number of unit cells in each direction are equal 共km3 with m as the number of unit cells in one direction兲. To gauge the effect of the shape of the simulation cell, we also simulated the n = 12 fcc system in a rhombohedral simulation box composed of rhombohedral unit cells, each containing one basis atom. The unit cell length on each edge is a = 冑2 / 3 with angle parameters ␣ =  = ␥ = 2 / 3. The system sizes studied with this simulation box geometry are N = 512, 1000, 1728, 2744, and 4096 particles, respectively. The potential parameters ⑀ and are each taken to be 1 in this study. Each simulation is run in the canonical 共NVT兲 ensemble at densities 3 = 1.1964, 1.3732, and 2.206 17 for n = 12, 9, and 6, respectively. Each density corresponds to melting condition10,11,17 at temperature kT / = 1.0. Each HTTP series starts at a low temperature 共kT / = 0.01兲 and perturbs into a higher temperature with temperature interval of 0.01, proceeding up to temperatures above the respective melting condition, ending at kT / = 1.2. In order to prevent the crystal structures from melting at temperatures above melting, we introduce a constraint potential on each atom, where the atom is always contained in the “enclosure” of its firstneighbor atoms’ shell, and there is no plane sliding allowed. This constraint potential is applied to simulations at all temperatures. We conduct 108 Monte Carlo 共MC兲 trial steps for each HTTP simulation. 50N MC trial moves are proposed as equilibration steps prior to any data collection. In each MC trial move, a pair of atoms selected at random is moved in directions opposite to one another so that the center-of-mass of the system is conserved. To improve the sampling time for systems with softness exponent n = 12 and 9, we apply a potential cutoff of 2.2 and 3.1, respectively. For n = 12, we test the truncation effect on the free-energy calculation for 864-particle system size at melting condition using different truncation radii of 2.2, 2.8, and 3.2, while 3.1, 3.3, and 3.6 are examined as the atomic truncation distance for the n = 9 soft potential. We find that the differences between the free energies calculated using these different atomic truncation radii are negligible, approximately 2 ⫻ 10−5 kT. However, for a longer-range soft potential 共n = 6兲, we apply halfa-simulation-box-length truncation 共=0.48ma兲 to sample the system. Also, in the context of minimizing the finite-size effect, the lattice sum used to calculate the lattice energy of each system and the Hessian for the harmonic free energy are computing using a cutoff of 58 and 15, respectively. The ␣C values are first estimated through a short simulation with 1 ⫻ 106 MC steps before computing the freeenergy difference with a longer simulation run. 108 MC trial steps are proposed in each simulation with the interval for data collection being the number of particles, N. The data are grouped into 100 blocks for error analysis. III. RESULTS AND DISCUSSION A. Results for correction term, AC calculation
In Fig. 1共a兲, we present the i values for the fcc crystal in a cubic box as a function of temperature kT / for different
J. Chem. Phys. 133, 134104 共2010兲
FIG. 1. 共a兲 Plot of i vs kT / for cubic fcc n = 12 soft potential at 3 = 1.1964. The symbols and lines represent the raw data and the fitted equation, respectively. The error bar in 共a兲 is smaller than the symbols. 共b兲 Plot of relative deviation of the simulation results and the fitted data, 共i − i,fitted兲 / i,fitted as a function of temperature, kT / , for N = 500 共black兲 and N = 4000 共violet兲 particles.
system sizes. The value of i is a measure of the free-energy difference between adjacent systems with temperatures Ti and Ti+1, and this does not necessarily relate to the degree of overlap of their important configuration spaces. However, the HTTP scheme developed here is designed to provide perfect overlap when applied to a harmonic system, so it is reasonable to expect that we will have a good degree of configuration-space overlap in the real systems to which the method is applied. In such a circumstance, we might look to i by itself as a rough measure of the difficulty of the HTTP calculation because a lot of the complications will be related to the asymmetry of the sizes of the two configuration spaces, and i can quantify this. A i value close to 1.0 indicates comparably sized configuration spaces, which coupled with an expectation of good overlap for the HTTP method should indicate an easier calculation. Figure 1共a兲 shows that i decreases with system size, indicating that, as expected, the perturbation between temperatures is harder for larger systems. The relative uncertainty, i / i, for each perturbation at any specific temperature 共with temperature inter-
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-6
Tan, Schultz, and Kofke
FIG. 2. AC / N vs kT / for cubic fcc exponent-12 soft potential at 3 = 1.1964. The inset is the blow-up of the curve at melting condition, kT / = 1.0.
val of 0.01兲 is approximately proportional to the system size. For example, for the perturbation between systems at kT / = 0.01 and 0.02, the relative uncertainties for i are about 0.3%, 0.8%, and 3% for N = 500 共smallest system兲, 1372, and 4000 共largest system size兲, respectively. Together with i from raw data, to demonstrate the quality of the fit described by Eq. 共9兲, we also plot i,fitted = exp共 −⌬AC兲 in Fig. 1共a兲, where AC comes from the fit. The simulation data are well described by the polynomial equation as can be seen for 500 and 4000 particles in Fig. 1共b兲. The relative deviation between the simulation data and the fitted line, 共i − i,fitted兲 / i,fitted, fluctuates about zero along the x-axis. The fluctuations are commensurate with the relative uncertainties. The root mean square of the deviation divided by the uncertainty ranges from 0.87 to 1.03. The results for the free-energy correction term, AC / N for different system sizes calculated using Eq. 共9兲, are presented in Fig. 2. We note that the correction term is nearly independent of system size, and in fact the results for different sizes cannot be distinguished on the plot. This behavior is consistent with what we found in our recent work.15 Nevertheless, an expanded scale reveals that AC / N decreases slightly with system size and appears to converge for the large system as shown in the inset of Fig. 2. To examine the system-size effect in more detail, we plot in Fig. 3 the correction term, AC / N, at the melting point as a function of 1 / N. We find the slopes to be 0.10共1兲 and 0.13共2兲 for the cubic and rhombohedral simulation box, respectively 共the last digit in parentheses represents the 67% uncertainty limit兲, which are about ten times smaller than the slopes when plotting the full free energy including the harmonic contribution.15,17 In Fig. 3, we observe that not only the correction term is barely dependent on the size of the system, the uncertainty for this term—as measured using the HTTP method—also does not vary much with N. We find that the values extrapolated to infinite-system size 共1 / N = 0兲 for both simulation box shapes agree well with each other, within their error bars.
J. Chem. Phys. 133, 134104 共2010兲
FIG. 3. n = 12. AC / N vs 1 / N at melting condition kT / = 1.0. The last digit in the parentheses is the 67% confidence limit.
B. Free-energy calculation results
In Fig. 4, we compare our results calculated using the HTTP method with previous work done by Polson et al.17 The plot presents the excess free energy per particle at the melting condition as a function of reciprocal of system size, where the excess free energy includes the harmonic contribution and is given with respect to an ideal gas as Aex / N = 共A − AIG兲 / N; the ideal gas free energy is AIG = ln − 1 + ln共2N兲 / 2N. The plot is constructed adding the term ln共N兲 / N to remove the leading-N dependence.17 We present results using our two box geometries, as well as the results reported by Polson et al., which used a rhombohedral box. Each result is then extrapolated to 1 / N = 0 to obtain its excess free energy for an infinite system. As we can see, the figure shows that the HTTP results for the two box shapes differ for finite-sized systems, but upon extrapolation, they are in excellent mutual agreement. Comparing our results with those of Polson et al., there is a slight disagreement of 0.002 observed. While this is a very small difference, it is significantly larger than 共more
FIG. 4. Excess free energy vs 1 / N at melting condition kT / = 1.0. The last digit in the parentheses and the error bars in the plot represent the 67% confidence limits; the error bars are not easily visible for the present results.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-7
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
TABLE I. Values and estimates of excess free energy 共Aex / N兲N→⬁ of hcp and fcc crystals for soft potential n = 12 at conditions 3 = 1.1964 and kT / = 1.0. The values reported in “This work” are results from hybrid method. Method This work Polson et al.a BenDBb HybridInfb
hcp
fcc
Difference, hcp− fcc
9.231 11共2兲 ¯ 9.235共2兲 9.231共2兲
9.228 31共2兲 9.226 08共7兲 9.231共2兲 9.226共2兲
0.002 80共3兲 0.0028共8兲 0.004共3兲 0.005共3兲
a
Reference 17. Reference 15.
b
than 30⫻兲 the confidence limits of either set of results. We have no definitive explanation for this discrepancy, though our suspicion is that it is due to the finite-step error inherent in the quadrature used for the Frenkel–Ladd method employed by Polson et al.—the difficulty may come, in particular, in handling the integrable singularity that arises in application to soft particles. Barring a mistake in our derivations or calculations and given the good phase-space overlap observed in our overlap-sampling calculations, there is no place in the HTTP method that could lead to an uncharacterized error that might cause this difference. Our preferred approach for calculating the infinitesystem free energy is the “hybrid method” that we presented in our previous study.15 In this approach, the free energy is calculated by performing separate extrapolations of the harmonic and anharmonic contributions to the free energy. The harmonic contribution has a greater system-size dependence, but we can calculate values for much larger systems, thereby improving the extrapolation. In contrast, it is much more difficult to calculate the anharmonic correction AC / N for large systems, but this contribution is much less system-size dependent 共at least for the soft-sphere systems examined so far; see Fig. 3兲. By separating the extrapolations, we can exploit the positive features of each to offset the respective difficulties. The extrapolated excess free energies 共1 / N = 0兲 from this hybrid method 共using HTTP for AC / N兲 are 9.228 32共2兲 and 9.228 30共2兲 for the cubic and rhombohedral simulation boxes, respectively. These results are consistent with the results presented in Fig. 4. We also calculate the free energy for n = 12 hcp solid at corresponding melting condition for the fcc crystal 共3 = 1.1964; kT / = 1.0兲. In similar fashion, we regress the data versus 1 / N as in Fig. 3 to yield the line 0.154 69共2兲 − 0.02共3兲 / N; adding in the extrapolated 共1 / N = 0兲 harmonic contribution yields the value given in Table I. Table I presents the infinite-system results from this work along with the findings15,17 that were previously reported. At melting condition, the fcc crystal is found to be more stable than the hcp crystal by an amount of 0.002 80共3兲. It is interesting to observe that this finding is in perfect agreement with the value reported by Polson et al.,17 even though their individual free energies differ from ours by an amount roughly equal to the hcp-fcc free-energy difference. We note that the results found using the HTTP method are more precise and conclusive than what we reported previously,15 especially in this application of determining the relative stability between different crystalline phases.
FIG. 5. The free-energy difference, 共⌬A / N兲N→⬁ between the competing crystal structures 共fcc, bcc, and hcp for different softness 共n兲 values. The Helmholtz free-energy difference, 共⌬A / N兲N→⬁, e.g., hcp− fcc, denotes 共A / N兲hcp,N→⬁ − 共A / N兲fcc,N→⬁. The subscript N → ⬁ indicates the extrapolated value to infinite-system size. The inset shows the bcc-fcc difference for n = 6 soft potential as a function of temperature. All results are collected at the same density for corresponding softness n values, which are 3 = 1.1964, 1.3732, and 2.20617 for n = 12, 9 and 6 respectively.
C. Relative stability of crystals
In this section, we use the HTTP method to examine the stability of different crystal structures 共fcc, bcc, and hcp兲 for different softness values, n. Lattice dynamics calculations10,11 show that the fcc and hcp crystals are stable 共or metastable兲 for all the softness values that we study here 共n = 12, 9, and 6兲, while the bcc structure is mechanically stable only for n = 6. For each crystalline phase, a constraint potential is imposed during simulation to the prevent the loss of its structure. To test the effect of the constraint potential on the free-energy calculation, we compare AC / N values calculated with and without the constraint at melting condition 共kT / = 1.0兲 using the 864-particle fcc n = 12 system. The effect is found to be negligible, with the calculated values being 0.157 10共3兲 共with兲 and 0.157 06共3兲 共without兲 for the two cases. The results for the crystal-stability calculation using the HTTP method are presented in Fig. 5. Also, see the Appendix for more calculations. Each curve reported in the plot is the infinite-system free-energy difference per particle between two crystals, where the absolute free energy of each crystal structure is determined as described above and the values at kT / = 0.0 indicate the difference in lattice energy between the two structures. We observe that the fcc crystal structure is more stable than the hcp structure for all values of the softness and temperatures, with fcc stability increasing with temperature. Not surprisingly, the free-energy difference between the fcc and hcp is small in comparison to the bcc-fcc difference. The free energy for the n = 6 system favors fcc at low temperatures, and the system undergoes a solid-solid phase transition to bcc at kT / = 0.802, in agreement with previous work, which found a transition temperature of approximately 0.822.13 Also, as can be seen from the inset of Fig. 5, the fcc-bcc difference, ⌬A = ⌬U − T⌬S, slowly decreases with temperature, and this indicates that the bcc
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-8
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
TABLE II. Values and estimates of excess free energy 共Aex / N兲N→⬁ for different crystalline structures 共fcc, bcc, and hcp兲. The excess free energy reported for n = 9 and 6 is at the melting condition 共3 = 1.3732 and kT / = 1.0兲 and the bcc-fcc phase transition 共3 = 2.206 17 and kT / = 0.802兲, respectively. Corrections AC are the results of this work 关Eq. 共9兲兴. n=9 Method Lattice energy Lattice energy+ harmonic Lattice energy+ harmonic+ correction共AC兲
n=6
hcp
fcc
Difference, hcp− fcc
fcc
bcc
Difference, fcc− bcc
6.401 648 46 11.631 837 06 11.689 03共2兲
6.401 293 69 11.627 759 84 11.683 72共2兲
0.000 354 77 0.004 077 22 0.005 31共3兲
22.612 353 77 28.332 440 12 28.314 05共2兲
22.716 955 46 28.303 564 27 28.314 07共4兲
⫺0.104 602 0.028 875 9 ⫺0.000 02共5兲
structure has more thermal fluctuation or higher entropy than the fcc structure over the temperature range. This effect becomes sufficiently important at high temperatures, where the transition occurs. With this HTTP method, we are capable of calculating the free energy for each crystalline phase with high precision and capturing the solid-solid transition with confidence. The largest uncertainty in the Helmholtz freeenergy difference between the crystalline structures shown in Fig. 5 is of the order of 10−5, which is very small on the scale of the plot. To better illustrate the entropic contribution to the relative stability between different crystalline structures, we tabulate in Table II the excess free energy for n = 9 and 6. For n = 9, we determine the relative stability between the hcp and fcc crystal at the fcc melting condition 共3 = 1.3732, kT / = 1.0兲, while the excess free energies of fcc and bcc crystals are calculated for the n = 6 system at the fcc-bcc phase transition 共3 = 2.206 17, kT / = 0.802兲. Table II decomposes each free-energy difference by considering contributions from the lattice energy alone, Ulat, and the lattice energy together with the harmonic contribution 关Eq. 共11兲兴. All the results shown in the table are infinite-limit system size and in excess of the ideal gas contribution, ln共兲 − 1. Corresponding results for n = 12 were presented in previous work.15 Table II shows for the n = 9 system that the lattice energy by itself favors the fcc structure by a very small amount 共0.005%兲. The addition of the entropic contribution to the free energy via lattice dynamics 关Eq. 共11兲兴 amplifies the difference significantly, increasing it by more than a factor of 10. The addition of the correction term further stabilizes the fcc structure but only marginally. The significance of these contributions becomes more prominent when we examine the n = 6 system at the fcc-bcc phase transition. The latticeenergy calculation, a crude approximation with complete neglect of the entropy, shows that the fcc crystal is more stable than the bcc structure by approximately 0.5%. However, the addition of the harmonic contribution switches the stability to bcc by 0.1%. The introduction of the correction term shifts stability back toward fcc, yielding the zero difference. The oscillation in polymorph stability with the introduction of various contributions demonstrates the potential importance of the entropic term in this application. D. Computational efficiency
In this part of the study, we investigate the computational efficiency of HTTP method as a function of temperature interval and number of simulation MC steps for the n = 12 system. The effect of simulation length is gauged by
collecting results using simulations of three lengths 共referred to here as short, medium, and long兲. Temperature intervals examined here are 0.002, 0.01, 0.05, and 0.2. For consistency, we pick the number of MC steps for each temperature interval such that the simulation runs for different temperature intervals employ equal amounts of computational effort by balancing the number of MC steps in each simulation and the total simulations required. Each series starts from the value of the temperature interval, e.g., 0.002, 0.01, 0.05, and 0.2, and proceeds in equal perturbation steps up to temperature, kT / = 1.2 共where again melting is at kT / = 1.0兲. Thus, the series for temperature intervals of 0.002 关0.2, 2, 20兴, 0.01 关1, 10, 100兴, 0.05 关5, 50, 500兴, and 0.2 关20, 200, 2000兴 comprises 600, 120, 24, and 6 simulations, respectively, across the temperature range. The numbers in the square brackets here indicate the total numbers of simulation steps 共in millions兲 considered for the short, medium, and long simulation runs. The calculations were run on an Intel Core2 Quad Q9550 Yorkfield 2.83 GHz processor. The computational efficiency results reported in Fig. 6 are for the 256- and 4000-particle system sizes at kT /
FIG. 6. Error measurement from Eq. 共12兲 vs computation time. The results are presented such that each curve represents results using different temperature intervals, e.g., 256-Int0.002 indicates the error measurement for 256particle system size with temperature interval of 0.002. The filled circle and triangle are for 256- and 4000-particle systems, respectively. Lines join points that apply the same method but use increasing 共from left to right兲 amounts of sampling at each temperature. TI and BenDB indicate the thermodynamic integration and harmonic perturbation methods 共Ref. 15兲, respectively. Each point within the TI curve differs by the number of MC steps at each temperature 共from left to right兲: 62.5, 125, 250, and 500⫻ 106. The solid and dashed lines are aids to the eye. The solid line with “1/2” indicates slope of magnitude 0.5.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-9
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
= 1.0. The accuracy of the HTTP method is calculated by comparing the correction term, AC, determined by the method under consideration, with what we take as the true value determined using temperature interval of 0.01 with 1 ⫻ 109 MC trial steps for each simulation. The error measure, ⌺, inclusive of accuracy and precision is given by
⌺=
冑 共AC − AC,true兲2 + 共SE兲2 , N
共12兲
where SE is the standard error of the mean in the final result. In Fig. 6, the three data points that are connected by the solid 共N = 256兲 or dashed lines 共N = 4000兲 are the error measurements determined from the short 共leftmost point兲, medium 共center兲, and long 共rightmost point兲 simulation runs, respectively, for different temperature intervals. As we can observe from the plot, the short, medium, and long runs 共of the same temperature interval兲 for both system sizes correspondingly amount to almost equal total computation hours of approximately 0.8, 8, and 80 h, respectively. However, a slight variation of about 7 h is observed in the total computation time between the 256- and 4000-particle system sizes for the long simulation run. We believe that the computation expense of the latter system has been greatly improved by the use of a potential cutoff of 2.2, for which the atomic interactions beyond the truncation radius are neglected. Figure 6 shows that 0.002, 0.01, and 0.05 are equally effective temperature intervals to determine the correction term, AC, using the HTTP method. The error measurement, ⌺, of these intervals decreases with the increasing number of MC steps, with a slope of 1/2. This suggests that the simulation results become more precise and converge to the true value 共more accurate兲, and the error for the long simulation comes primarily from the imprecision of the results. In contrast, for the temperature interval of 0.2, the computation efficiency does not improve with increasing calculation efforts. The error is mainly dominated by the inaccuracy, which records differences of 0.0007 and 0.004 from the true value, for N = 256 and 4000, respectively, for the long simulation run. The inaccuracies of these results outweigh their standard error by an order of magnitude, and the large error of this interval is well illustrated in Fig. 6. This large error is mainly caused by the poor perturbation between the two systems at different temperatures, where its temperature interval is apparently too large to have a good overlap in phase space between two systems. The problem worsens with increasing system size. Another reason is the extrapolation of the data to zero from temperature at kT / = 0.2, which seems to be a less appropriate procedure to follow in this case. One way to overcome the latter problem is to introduce an extra point at kT / = 0.05 in order to have a better regression of the data to zero. Judging from the plot, the temperature intervals of 0.01 and 0.05 are reasonable choices for this study because both intervals have a balance between 共1兲 a good amount of overlap in phase space between the perturbing systems and 共2兲 the number of simulations required. Also in the plot, we present the results calculated using
TI in temperature and the harmonic perturbation method 共BenDB兲, both of which were demonstrated in a previous work.15 The error measures at the melting condition from both methods are approximately two orders of magnitude larger than the one calculated using the HTTP method. In the BenDB method, for the 108-particle system size, the error is mainly contributed by its uncertainty, which is 0.0017. The result for the 256-particle system suffers from both its inaccuracy and imprecision; this is where we believe that the limit of the BenDB method is reached. TI shows an error of the order of 10−3. It improves slightly with increased computation, reflecting that it is prone to error introduced by both the quadrature method as well as stochastic errors in the integrated data. It should be noted that the BenDB method provides the free energy at only one point, in contrast to the TI and HTTP methods, which provide the free energy over a range of temperatures as a by-product of the calculation examined in Fig. 6. On the other hand, the TI and HTTP methods are unable provide absolute free energies without taking this path, while the BenDB method yields an absolute free energy without requiring the possibility of connecting to zero temperature. To demonstrate the importance of harmonic targeting as part of the overall HTTP method, we also employed the Bennett-optimized overlap-sampling FEP method to calculate AC / N between systems at different temperatures without performing scaling of the atomic positions. For both 500and 4000-particle system sizes with a temperature step of 0.01, we find that relative uncertainty of i 共or equivalently, the absolute error in the free energy per particle兲 is about an order of magnitude larger than that found using the HTTP method applied to the same pair of temperatures. This is true for perturbation steps over almost all of the temperature range and both system sizes. Further, the nontargeted overlap-sampling perturbation fails completely—with relative uncertainty of i being more than 70%—in the lowtemperature region 共e.g., perturbing from temperatures of 0.01–0.02兲, where the ratio of the phase spaces of the two perturbing systems is large. In addition to the poor precision, we also observe that the low-temperature perturbation results 共without harmonic targeting兲 are prone to inaccuracy. The comparison between targeted and nontargeted temperature perturbation is likely to be only worse for temperature steps larger than the 0.01 increment examined here. IV. CONCLUSIONS
We have demonstrated the effectiveness of harmonically targeted temperature perturbation, coupled with Bennettoptimized overlap sampling, to calculate the free-energy difference between systems at different temperatures, and further to compute the 共classical兲 absolute free energy by connecting to a harmonic system at zero temperature. The demonstration here is for a simple soft-sphere model, considering three values of the softness parameter. The method gives directly the free energy in excess of a harmonic system at the same temperature, and we show again that this quantity 共AC兲 is much less dependent on system size than the harmonic contribution is to the free energy and thus the total
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-10
J. Chem. Phys. 133, 134104 共2010兲
Tan, Schultz, and Kofke
TABLE III. Values for the lattice energy and vibrational contribution to the free energy.
n = 12
fcc hcp fcc hcp fcc bcc hcp
n=9 n=6
共Ulat / N兲N→⬁
共Avib,0 / N兲N→⬁
1.516 485 025 1.516 536 721 2.208 391 122 2.208 528 128 3.613 475 382 3.630 711 328 3.613 718 378
3.888 452 45 3.891 589 92 3.482 175 24 3.485 900 58 2.698 198 34 2.564 731 10 2.702 699 96
free energy itself. This HTTP method is capable of providing results of very high precision. The good precision obtained here allows us to determine the relative stability between different crystalline structures with high confidence. The results obtained here neglect the effects of thermodynamically stable defects, vacancies in particular, which have contributions greater than the errors in our calculated free energies.15 The effects of vacancy defects are however negligible compared to the harmonic and anharmonic contributions 共which are at least an order of magnitude larger兲. Also, a previous study18 has demonstrated that the relative stability of these crystal structures is not affected by the presence of defects. For soft potential with n = 12 and 9, we show that the fcc crystal is more stable than the hcp structure at all temperatures up to melting. The hcp-fcc free-energy difference increases substantially with temperature, indicating the significance of the entropic effects in determining the relative stability of crystal structures at nonzero temperature condition. This conclusion is further supported by the emergence of the bcc crystal structure at kT / = 0.802 for the n = 6 soft potential, where the fcc structure gives way to bcc crystal. The focus of this work is on the soft-sphere model, first to demonstrate the method, and second so that we can establish highly precise values for its free energy. Extension to more complex multiatomics is an important next step in demonstrating the utility of this approach.
ACKNOWLEDGMENTS
Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for support of this research. Computational resources were provided by the University at Buffalo Center for Computational Research.
APPENDIX: FREE ENERGY FORMULAS
For completeness, we present the expressions that we have determined for the free energy of the defect-free softsphere system of different exponent n values: 12, 9, and 6 with various crystal packings 共fcc, bcc, and hcp兲. In the thermodynamic limit, the configurational free energy 共apart from momentum contributions兲 can be expressed as
A Ulat Avib AC = + + . N N N N
共A1兲
The static-lattice-energy contribution is given in terms of density and temperature according to
冉 冊
Ulat 1 Ulat,0 = N ⌼ N
共A2兲
, N→⬁
where again ⌼ = 共kT / 兲共3兲−n/3 is the single independent state parameter for the model. The term in parentheses is the infinite-system lattice sum energy at unit density and is tabulated for the various cases in Table III. The vibrational contribution in Eq. 共A1兲 is the limiting value of the sum in Eq. 共11兲, and it depends on temperature and density according to
冉
Avib Avib,0 = N N
冊
3 − ln ⌼ + ln共3兲. 2 N→⬁
共A3兲
The temperature- and density-independent term Avib,0 / N also is presented in Table III. Finally, the anharmonic contribution Ac is given as a function of temperature and density via the group ⌼ in Eq. 共9兲. The state-independent constants ci are presented in Table IV. The infinite-system configurational free energy can be computed using these formulas for the corresponding crystal at any temperature and density up to the melting point 共values of ⌼ below about 0.5, 0.4, and 0.2 for n = 12, 9, and 6, respectively兲, with the coefficients given in Tables III and IV.
TABLE IV. Fitted coefficients for 共AC / N兲N→⬁ as defined in Eq. 共9兲.
n = 12 n=9 n=6
fcc hcp fcc hcp fcc bcc hcp
c1
c2
c3
c4
c5
c6
c7
0.464 148 0.462 065 0.249 526 0.250 099 ⫺0.024 270 0.283 913 ⫺0.023 687
⫺0.423 821 ⫺0.413 562 ⫺0.245 724 ⫺0.282 838 ⫺0.385 254 ⫺2.421 458 ⫺0.351 963
0.466 322 0.429 175 ⫺0.500 979 0.007 132 ⫺1.989 540 14.558 777 ⫺2.934 478
⫺0.526 094 ⫺0.461 716 4.258 323 0.922 672 24.366 184 ⫺77.105 689 37.070 278
0.148 125 0.111 491 ⫺16.542 027 ⫺5.218 832 ⫺208.856 577 229.711 416 ⫺280.028 273
0.401 893 0.368 595 30.811 960 10.890 988 860.698 810 ⫺391.161 702 953.287 096
⫺0.384 095 ⫺0.352 244 ⫺23.208 566 ⫺7.3695 08 ⫺1469.431 853 335.089 915 ⫺1140.820 226
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
134104-11 1
J. Chem. Phys. 133, 134104 共2010兲
Solid-phase free energies
P. A. Bash, U. C. Singh, F. K. Brown, R. Langridge, and P. A. Kollman, Science 235, 574 共1987兲; L. X. Dang, K. M. Merz, and P. A. Kollman, J. Am. Chem. Soc. 111, 8505 共1989兲. 2 D. A. Kofke, Fluid Phase Equilib. 228–229, 41 共2005兲. 3 D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. 共Academic, San Diego, 2002兲. 4 R. W. Zwanzig, J. Chem. Phys. 22, 1420 共1954兲. 5 N. D. Lu, D. A. Kofke, and T. B. Woolf, J. Comput. Chem. 25, 28 共2004兲. 6 N. D. Lu, J. K. Singh, and D. A. Kofke, J. Chem. Phys. 118, 2977 共2003兲; N. D. Lu, D. Wu, T. B. Woolf, and D. A. Kofke, Phys. Rev. E 69, 057702 共2004兲. 7 D. Wu and D. A. Kofke, J. Chem. Phys. 123, 054103 共2005兲; 123, 084109 共2005兲. 8 C. H. Bennett, J. Comput. Phys. 22, 245 共1976兲. 9 C. Jarzynski, Phys. Rev. E 65, 046122 共2002兲.
10
W. G. Hoover, S. G. Gray, and K. W. Johnson, J. Chem. Phys. 55, 1128 共1971兲. 11 W. G. Hoover, M. Ross, K. W. Johnson, D. Henderson, J. A. Barker, and B. C. Brown, J. Chem. Phys. 52, 4931 共1970兲. 12 R. Agrawal and D. A. Kofke, Mol. Phys. 85, 23 共1995兲. 13 S. Prestipino, F. Saija, and P. V. Giaquinta, J. Chem. Phys. 123, 144110 共2005兲. 14 D. C. Wang and A. P. Gast, J. Phys.: Condens. Matter 11, 10133 共1999兲. 15 T. B. Tan, A. J. Schultz, and D. A. Kofke, J. Chem. Phys. 132, 214103 共2010兲. 16 M. T. Dove, Introduction to Lattice Dynamics 共Cambridge University Press, New York, 2005兲. 17 J. M. Polson, E. Trizac, S. Pronk, and D. Frenkel, J. Chem. Phys. 112, 5339 共2000兲. 18 S. K. Kwak and D. A. Kofke, J. Chem. Phys. 122, 176101 共2005兲.
Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
