Configurational mapping significantly increases the efficiency of solid-solid phase coexistence calculations via molecular dynamics: Determining the FCC-HCP coexistence line of Lennard-Jones particles

Polymorphism, or the ability of a crystal to pack into multiple metastable states, is important in materials study and design. Polymorphism affects the properties of materials such as charge transport¹ and bioavailability.^2,3 When multiple metastable polymorphs with different properties are present, the calculation of solid-solid coexistence curves becomes important. Temperature and pressure transformations are present in materials such as pharmaceuticals^4,5 and metals.^6,7

Traditional phase-coexistence calculation methods, such as the Gibbs ensemble method,^8–10 are either not applicable to solid-solid systems or require a previously known coexistence point and suffer from increasing error due to the use of numerical integration.^11,12 We have previously introduced the Successive Interpolation of Multistate Reweighting (SIMR) method to predict solid-solid phase diagrams.¹³ This methodology does not rely on lattice dynamics and is thus applicable in systems that are far from harmonic. It calculates the phase diagram from direct calculation of the relative Gibbs free energy using a series of direct molecular dynamics or Monte Carlo simulations without any specialized sampling techniques. It can thus be wrapped around any molecular simulation code.

One drawback to this methodology is that it requires an overlap in the energy and volume phase space between adjacent temperature and pressure simulations. This presents a challenge in a number of situations, for example, when an extremely large pressure range is desired. Here, we present an extension of the SIMR method using a configurational mapping technique, inspired by the work of Tan and collaborators,^14–16 which reduces the number of simulations required and therefore the computational cost.

We have applied this method to the solid-solid phase diagram of Lennard-Jones (LJ) spheres, a common test systems in molecular simulation. The Lennard-Jones potential is often used to approximate the solid phase of noble gases such as argon, as well as the highly spherical methane, and to test methodologies for more chemically complex solids.^17–19 Many methods have been used to successfully and accurately calculate the melting and vaporization lines of Lennard-Jones system,^20–28 such as the Gibbs ensemble method^9,21,22 (for vapor liquid) and thermodynamic integration.²³ However, it is significantly harder to calculate solid-solid phase equilibria.

In this paper, we focus on solid infinite crystal systems. While the coexistence line of the hexagonal close-packed (HCP) and face-centered cubic (FCC) phases of Lennard-Jones spheres has been calculated using a variety of approximations and methods, there is substantial variation in the results from these studies. All studies agree that the two stable phases of the solid Lennard-Jones spheres are extremely close in free energy, on the range of ΔG = 1 − 10 × 10⁻⁴ per particle in reduced units throughout most of the phase diagram, meaning that very small uncertainties or errors result in large changes and uncertainties in the phase diagram.

A number of research groups has attempted to calculate the Lennard-Jones FCC/HCP phase diagram with a range of approximations, with generally inconsistent results. Choi et al.²⁹ performed an early calculation using perturbation theory of the Lennard-Jones potential around the hard-sphere close packing to determine the phase boundary between the FCC and HCP solids and the liquid phase. van der Hoef’s equations were then used to obtain the residual Helmholtz energy.³⁰ The configurational Helmholtz free energy can then be expanded as a perturbation series, and thermodynamic properties can be calculated from the first two terms.^31,32 However, as shown in Fig. 1, this approach was not consistent with the later, more comprehensive approaches. The coexistence line between FCC and HCP LJ structures has also been calculated using dynamic lattice theory (DLT).³³ However, this approach for obtaining the free energy uses a harmonic approximation, which is not valid at the level of the accuracy needed here.

Calculation of the fully anharmonic Helmholtz free energy has previously been performed using Monte Carlo simulations.^34–36 In the work of Adidharma and Tan,³⁵ canonical ensemble Monte Carlo simulations were performed at a variety of reduced temperatures and densities. The results from the simulations were then fit, using the energy and pressure from the simulations, and constants derived by Stillinger,³⁷ to the equation for the Helmholtz free energy of the Lennard-Jones solid derived by van der Hoef.³⁰ From the Helmholtz free energy, the coexistence line was then determined.

Lattice-switch Monte Carlo is another approach that has been used to calculate the free energy difference between phases. In the lattice-switch Monte Carlo method, a transformation between phases is proposed, which takes the molecules of one structure and converts the atomic positions and box vectors to those of the other phase. A range of multicanonical approaches must be applied in order to get sufficient exchange between the packings.^38,39 Once the free energy difference was calculated as a function of T^* and P^*, the phase boundary is easily determined. This method has been used to calculate the FCC-HCP coexistence line³⁶ and to study the instability of the BCC phase relative to FCC.⁴⁰ The results of Jackson are more consistent with later methods but are still temperature shifted, especially for smaller box sizes; for larger box sizes, the method was too inefficient to run at higher densities.

One recent comprehensive study of the LJ phase diagram was an extension of the earlier dynamic lattice theory approach of Travesset.^33,41 Rather than directly calculating the full free energy, the anharmonic contribution to the free energy was calculated using molecular dynamics with thermodynamic integration along a switching parameter λ, where λ = 0 corresponds to the harmonic potential energy U^DLT and λ = 1 is the full Lennard-Jones potential, as seen in Eq. (16) of Calero et al.⁴¹ Using the DLT energy and 20–50 simulations with a mixed potential λ value, the anharmonic contribution is calculated by thermodynamic integration. Once the harmonic (via DLT) and anharmonic contributions to the free energy have been calculated, conversion was performed to find the corresponding pressure for use in the calculation of the temperature-pressure coexistence curve. The addition of the anharmonic free energy term to the previous DLT results illustrates how a small change in free energy value results in a large change in the coexistence curve, as seen by the difference in the “Travesset” and “Calero” lines in Fig. 1. However, even more recently, Schultz and Kofke²⁸ published a comprehensive study of the Lennard-Jones phase diagram, using both exhaustively extensive direct simulation and analytic quasiharmonic approaches to explore the entire phase diagram of Lennard-Jones particles. Interestingly, the FCC-HCP equilibria results of Schultz and Kofke are in stronger agreement with Jackson’s results than the later results of Calero, et al.⁴¹ 

A comparison of many of these previous FCC-HCP coexistence lines over the full range of methods is shown in Fig. 1, which shows large differences between methods. All data were taken directly from temperature-pressure phase diagrams present in the papers using WebPlotDigitizer,⁴² with the exception of the Schultz et al. results, which were taken from correlations lines 7 and 9 of Table I of Schultz and Kofke.²⁸ We note that, as verified with the authors, the $v4$ which should be in the denominator of line 9 is incorrectly written as $v2$⁠. This was a transcription error in the table alone, not an error in the results of the study. However, none of the studies above published explicit error bars, making comparisons particularly difficult.

To obtain the phase diagrams of Lennard-Jones spheres using full molecular dynamics, we used the Successive Interpolation of Multistate Reweighting (SIMR) method.¹³ This method combines the reduced free energy difference values between temperature and pressure states within a polymorph, defined as βG, and a reference Gibbs free energy difference between polymorphs at the same temperature, which can then be combined to obtain the Gibbs free energy difference between polymorphs at all temperatures and pressures in the region of interest.

The reference Gibbs free energy difference value is determined using the pseudo-supercritical path (PSCP) method.^43–45 This method determines the free energy required to take each polymorph from a real crystal to an ideal gas; the free energy between the two polymorphs is then the difference of those values.

The reduced free energy between states within a polymorph is found using the multistate Bennett acceptance ratio (MBAR).⁴⁶ This method uses Eq. (1) to iteratively solve for the reduced free energy of each state f_i with respect to each other state f_k, where N_k is the number of configurations drawn from state k and u_k(x_jn) is the reduced energy of configuration n sampled in state j and evaluated in k. The reduced energy is defined as u_k(x_jn) = β_kU(x_nj) + P_kV(x_nj).

Using the definition of reduced free energy as given above, the Gibbs free energy difference between the two polymorphs at state i is given as Eq. (2), where Δf_ij is the difference in reduced free energy between states f_i and f_j and T_ref is some reference temperature where ΔG_ij is known. Linear interpolation is then used to find the points where the difference between polymorphs is zero, which is coexistence. The uncertainty in the coexistence points using this method is found in Eq. (3), where δd is the magnitude of the uncertainty in the coexistence line perpendicular to the line and δΔG is the uncertainty in the free energy at a point along the coexistence line. Full details of this method can be found in Schieber et al.¹³ In theory, reweighting can be used to refine estimates in between coexistence points rather than direct interpolation, but this is less reliable in the case of configuration mapping, as described below.

One requirement for simulations used for the SIMR method is that simulations adjacent in temperature or pressure have a non-negligible amount of phase space overlap,¹³ as defined in Eq. (4). Conceptually, this means that simulations have some set of configurations that they both sample. The overlap between states 1 and 2 is then dependent on the probability of all configurations, x, in each of the two distributions, P₁ and P₂. Due to this requirement, the number of simulations performed is dependent on the width of the energy and volume distributions of the simulations. Systems with wider potential energy and volume distributions are likely to still achieve phase space overlap with wider spacing. The width of these distributions, and thus the spacing in temperature and pressure that is allowable between simulations, depends on factors such as temperature, pressure, and the size and flexibility of the molecule. In order to decrease the number of simulations and therefore the computational resources required, it is desirable to increase the spacing between sampled states by increasing phase space overlap between states O_1,2, defined by

One potential way to increase phase space overlap between states is configuration mapping.^{14–16,47,48} Configuration mapping transforms the set of coordinates in one thermodynamic state into a set of coordinates that is more likely to have low energy in the other thermodynamic state of interest and evaluates the energy in the new state with the transformed configuration rather than the originally sampled configuration. We can then analytically calculate the free energy change for performing this mapping. This approach was used by Tan et al.¹⁴ to calculate the temperature dependence of the free energy of solids, by Paliwal and Shirts⁴⁷ to calculate the Gibbs free energy of transformation between different water models,⁴⁷ and in a differential form to calculate physical properties such as the heat capacity of HCP iron and the dielectric constant of the Stockmayer potential.¹⁵ Configuration mapping was shown to significantly improve the precision of these calculations. Here, we propose to use this methodology to improve the phase space overlap for SIMR.

Mathematically, we define a transformation T(x), with Jacobian J(x), which is allowed to depend on the current configuration x, and ΔU(x) is the difference in potential energy between the configuration when evaluated in the original and mapped states,¹⁶ ΔU(x) = U(T(x)) − U(x). In general, the potential can also change,⁴⁷ but in this study, we only change the temperature and pressure between states.

For a simple one-step transformation using the Zwanzig equation, the Helmholtz free energy difference in terms of mapping can be written as

More generally, when using configuration mapping to calculate free energy differences with multistate reweighting, we can derive equivalent formulas by replacing the reduced energy u(x) = βU(x) + βPV used in MBAR⁴⁶ or BAR⁴⁹ in the NPT ensemble with a “warped” reduced energy, defined in Eq. (6) by analogy with “warped bridge sampling,” a version of this technique used in statistics⁵⁰ to calculate the free energy difference between the states. In Eq. (6), i is the state the configuration was drawn from, j is the target state the energy is evaluated in, T_ij(x) is the transformed set of coordinates that were sampled from state i, and |J_ij(x)| is the determinant of the Jacobian of the transformation T_ij.

Equation (6) is applicable to all transformations that have a nonsingular Jacobian, but only a relatively small proportion transformations actually increase phase space overlap. A good transformation is one that is both relatively easy to implement and results in significant increase in the overlap. With a good transformation, the efficiency of the free energy calculation from state i to j can be made much more efficient.

In this study, we applied this coordinate mapping approach to the system of Lennard-Jones spheres to map between states defined by different temperature and pressures. In the case of point particles such as Lennard-Jones spheres, we only need to map the location of the particles themselves and not deal with any internal degrees of freedom. This scaling is applied between any two pairs of simulations, as the average box vectors change through both thermal expansion and compression/expansion due to changes in pressure.

To implement this type of transformation, first, the trajectories from states i and j are read. The desired transformation between the two states is determined, which usually requires some information from both trajectories to be efficient. This transformation is then applied to the coordinates in the trajectory i. In our case, these new coordinates are written to a new trajectory file and the energy of each frame of the trajectory is then reevaluated. The warped reduced energy is calculated using the warped energy and the new P and T, and the reduced free energy is calculated using Eq. (1). For the multistate reweighting process used in the SIMR, this transformation/reevaluation is performed for every set of pairs states in the (not necessarily regular) T, P grid.

The transformation we used is defined by the following:

where $r→iw$ is the new coordinate in the target ensemble, $r→i$ is the original coordinate, and B_i and B_j are the average box vector matrices of the original and target trajectories. If the two simulations differ in temperature, we also scale the deviation of the particle from its equilibrium position as derived by Tan, Schultz, and Kofke.¹⁴ This temperature scaling is carried out using the reduced coordinates $r→i,r$ or the fractional coordinates within the box, thus making this a three-step process. In Eq. (7), $Δr→i$ is the magnitude of the deviation of the molecule from its equilibrium position, $Δr→i=r→i−⟨r→i⟩$⁠, and T_i and T_j are the temperatures of the initial and target trajectories, respectively. The energy contribution of this Jacobian is $3N−32logTjTi+3N⁡log|Bj||Bi|$⁠, where N is the number of particles. The “−3” occurs because of the removal of translational center of mass motion in the simulation. Note in this case, a constant Jacobian is used for all configurations, though in theory it can be configuration dependent. Once all $uijw$ values have been calculated, Eq. (1) can be used directly to calculate the reduced free energy differences between the states.

An example of the effect of this mapping on the energy-volume distribution for two temperature and pressure states in this LJ system can be seen in Fig. 2. This figure shows the difference between the overlap in energy and volume achieved between the two states unmapped and using configuration mapping. The two unmapped trajectories show no overlap of their energy and volume distributions. The mapped distributions, however, show significant overlap in energy and density, which in most cases will translate directly to the overlap of configuration phase space.

The cost of energy reevaluations is very low compared with the cost of simulations. For example, in one test on a standard laptop with GROMACS,⁵⁴ the cost of mapping and reevaluating uncorrelated samples between the two states from a 9 ns simulation costs approximately 7.4 CPU-min. The cost of running one LJ particle mesh Ewald (PME) simulation itself for 9 ns is approximately 56 CPU-h. For a set of 187 states, as is examined here, the mapping cost between all pairs of states is then ∼4300 CPU-h, which is the cost of about 75–80 additional simulations. However, using mapping, we need at least 15–40× less simulation than we otherwise would have needed with SIMR, as discussed below, though the exact numbers will depend significantly on the specific code used.

When this particular mapping was applied to the system of Lennard-Jones spheres, the number of states required to achieve overlap decreased significantly. Without mapping, the minimum pressure spacing required to achieve sufficient overlap for the calculation of MBAR to converge with a finite value was approximately 1 P^*. With mapping, this could be increased to 25 P^* for the cutoff-based Lennard-Jones simulations and 12.5P^* for the higher precision PME calculations, which translates directly into an efficiency increase of 12–25 by decreasing the number of simulations required. Although overlap in the temperature dimension was already enough to achieve MBAR convergence with simulations every 0.066T^* (0.033T^* for PME simulations), mapping decreased the uncertainty of free energies between neighboring states in the T direction by roughly a factor of 1.31–1.48 (for example, at P^* = 127 with simulations spaced at T^* = 0.006 or T^* = 0.026, respectively), leading to an additional efficiency improvement in between ∼1.31² ≈ 1.71 and ∼1.48² ≈ 2.19, for an overall efficiency gain of 15–40×.

The Lennard-Jones phase diagram was produced using a system of 1200 LJ spheres and the standard Lennard-Jones 12-6 potential. The systems were set up with 10 layers of atoms in the x and y directions and 12 layers of atoms in the z direction, for a total of 300 FCC unit cells and 200 HCP unit cells. In a limited study by Jackson et al.,³⁶ a change in size from 216 to 1728 gave rise to a shift of about 0.1 T^* in the FCC-HCP coexistence line. Since our study is nearer the upper end, size dependence will be relatively small compared with the uncertainty, as discussed in more detail later when we look at the possible reasons for differences from previous results. The Lennard-Jones parameters for OPLS-UA methane were used in the simulations themselves (σ = 0.373 nm, ϵ = 1.2304 kJ/mol, m = 16.043 amu),⁵¹ though all results are reported in reduced units. The range of the phase diagram was chosen to correlate with previous lattice dynamics studies of Lennard-Jones spheres.³³ The temperature in reduced units was 0.066–0.466, and the pressure was between 0.003 and 508.9, which corresponded to 10–70 K and 1–200 001 bars in real units. This temperature range was chosen in order to include the region of predicted coexistence without including the region of melting. The pressure range was chosen to include the maximum HCP temperature stability point and the reentrant behavior and to include high enough pressures that the coexistence line is approximately linear and can be extrapolated. Simulations were initially spaced every 0.066 T^* and 25.44 P^*. For PME simulations, simulations were spaced every 0.033 T^* and 12.72P^*. The largest current limitation in the number of states that can be simulated is the limit on the memory that is available to run MBAR, which will fail with a memory error in the current implementation of pymbar⁴⁶ used to solve if the input matrix is too large.

All production molecular dynamics simulations of Lennard-Jones spheres were performed with GROMACS 5.1.2,^52,53 using a velocity Verlet integrator and Nosé-Hoover temperature control⁵⁴ with a time constant of 1.485 reduced units (2 ps). Isotropic Martyna-Tobias-Tuckerman-Klein (MTTK)⁵⁵ pressure control with a time constant of 7.24 reduced units (10 ps) was used. This combination of integrator and pressure control was chosen because it was shown to be the most stable for NPT simulation of small LJ systems in GROMACS at high pressure. For PSCP simulations, the sd integrator (Langevin dynamics) and thermostat and the Parrinello-Rahman barostat⁵⁶ were used to avoid nonergodicities at the fully restrained and noninteracting states of the PSCP. All particle mesh Ewald (PME) simulations were run for 9 × 10⁶ steps at a time step length of 0.00297 t^* (4 fs) for a total simulation time of 26 728 t^*. The cutoff used in these PME simulations was 3.5 σ. Potential switch simulations were run with the same time step for 11 879 t^*. Two types of simulations were run to determine the phase diagram: with cutoffs and using particle mesh Ewald (PME). For cutoff-based simulations, the van der Waals interactions were treated with a potential switch cutoff of 1.119 or 0.9325 nm, which corresponds to 3.0σ and 2.5σ, respectively (in units of Lennard-Jones radius). This method smoothly switched the potential over a range of 0.02 nm as a function of radius down to make the potential at the cutoff 0 using the vdw-modifier, potential-switch keyword. Because of the smaller relative energy difference between phases, all PME simulations were run for a factor of 2.5× longer than all simulations with the potential cutoff to add statistical accuracy.

Rather than increasing the cutoff size and extrapolating to infinity, we used particle mesh Ewald for dispersion interactions to incorporate long-range effects. This method works by only calculating the short-range interactions directly. Long-range interactions are calculated with a 3D fast Fourier transform on a grid.^57,58 The smooth version of this method, as implemented by Essmann et al.,⁵⁷ uses a B-spline interpolation on a grid to increase computational efficiency. This method was implemented for the dispersion term in GROMACS by Wennberg et al.^59,60 For all PME simulations, a Fourier grid of n_x = 36, n_y = 32, and n_z = 36, and a PME grid interpolation order of 6 were used. A cutoff of the direct-space summation was used and shifted from 3.47319σ (1.2955 nm) to 3.5σ (1.3055 nm) to guarantee smooth integration of the equations of motion, with tolerance of the direct space error at the cutoff (ewald-rtol-lj) equal to 1.0 × 10⁻⁶.

Two PSCP calculations were carried out for each calculation, one at 127 P^* and 0.33 T^* and the other at 152 P^* and 0.4 T^* in the NVT ensemble with a PV correction term of PΔ⟨V⟩ to convert between Helmholtz and Gibbs free energy. Simulations were carried out in NVT at the average volume corresponding to the equilibrated corresponding NPT simulation. The (127 P^*, 0.33 T^*) phase point was used for the phase diagram calculation, with the other point used to check cycle closure. In the PSCP, intermolecular interactions were turned off quartically, while simultaneously harmonic restraints were added to the average lattice positions. A set of 25 intermediates was used, with the force constant turned on quadratically to a maximum value of 113.076 in reduced units (1000 kJ mol⁻¹ nm⁻²), as per the protocol of Dybeck et al.⁶¹ 

Our more challenging phase diagram, generated using PME, is the result of one set of additional low pressure and temperature simulations, which doubled the number of simulations in the region below 200P^* and 0.3T^*, after the initial simulations. In this case, unsimulated states, where the average box vectors and displacement vectors are not known, cannot be incorporated into MBAR because the mapping procedure requires knowledge of the equilibrium P and T at each state, which must be obtained from the simulation. To calculate the uncertainty, 200 bootstrap samples of simulation configurations were generated. The uncertainty in the ΔG per particle between the FCC and HCP phases at each point was then determined from the standard deviation of the reduced free energy using the bootstrapped input. This was then converted to the uncertainty in the width of the coexistence line and plotted perpendicular to the line, as errors along the line do not actually affect the line, as described by Schieber et al.¹³ Alternately, using each of the bootstrapped free energies, a new phase equilibrium line for each bootstrap resampling line can be drawn. Error bounds can then be generated by taking the lines that represent a given confidence interval away from the median line. This is most accurately done in the perpendicular direction, though because in this case the phase line is mostly vertical and noise in the line makes tangent determination challenging, we look at the confidence interval in the temperature dimension. We use a linear spline for the phase equilibrium points obtained using SIMR to sort the points at any temperature of interest. A comparison of the two types of error analysis can be found in section II of the supplementary material (using a 1σ confidence interval). We note very good consistency between the methods above the 100P^* line; below the line, the bootstrap estimate becomes inconsistent because of low overlap, leading to large fluctuations in the bootstrap uncertainty.

Thermodynamic cycle closure was used to validate the process for the PME line, using two separate PSCP calculations at two different points. At 127 P^* and 0.40 T^*, the PSCP value is −0.000898(9). At 127 P^* and 0.27 T^*, the PSCP value is 0.000166(9). Using the second PSCP at the reference value in SIMR, the calculated ΔG at 127 P^* and 0.40 T^* is −0.00086(6), which is within the uncertainty of the PSCP value at that point, indicating cycle closure to high precision. After initial simulations and one round of additional simulations at areas of low overlap, at T^* = 0.20, the average uncertainty in Δμ was 0.0004 reduced units, which is similar to that shown in Fig. 6 of Calero et al.,⁴¹ which ranges from approximately 0.00025–0.0006, as seen in our Fig. 6. All free energies were calculated using the pymbar 3.0.3 implementation of MBAR.⁴⁶ 

We see in Fig. 4 that using mapping, the coexistence line was determined with good precision. The phase coexistence line is significantly affected by the treatment of the cutoff. Cutoffs should not affect the liquid-vapor transition significantly since those phases are essentially uncorrelated outside the cutoffs investigated here. However, in solids, cutoff effects of calculations using LJ spheres are important.^41,62 The coexistence line predicted by the SIMR and configuration mapping using a 2.5σ cutoff has a region of HCP coexistence that is higher in pressure and temperature and spans a wider range of pressures than the coexistence line predicted using a 3.0σ cutoff, as seen in Fig. 3. Both these results are a poor match for literature results extrapolated to a long cutoff in the high pressure region, as seen in Fig. 4. The reentrant behavior resulting from the SIMR method is much sharper than literature results with extrapolated large cutoffs, and the high-pressure and low-temperature region is not well approximated by this method. The poor match of the coexistence lines using a potential switch cutoff is due to the uneven increase in the contribution of the dispersion energy as the pressure is increased. The sudden inclusion of an entire shell of atoms under the cutoff value as pressure increases causes nonphysical behavior in the energy difference between phases. This is a known effect, analyzed, for example, by Jackson et al.,³⁶ and we analyze the reasons and resulting issues in more detail in Sec. I of the supplementary material (Fig. 5).

Using particle mesh Ewald for the dispersion term, we obtain a more accurate molecular dynamics phase diagram of FCC and HCP LJ spheres without the inconsistencies introduced by a potential cutoff. This phase diagram shows the same stability trends and reentrant behavior as literature results which in theory include the full statistical mechanics. The resulting PME SIMR phase diagram is approximately a standard deviation away from the results of Calero et al.,⁴¹ which at P^* = 250 is an average of 0.0234T^* and is within the uncertainty for most of the pressure range. The maximum HCP temperature stability is somewhat higher than their results. We note that although Calero et al. did not plot uncertainties in their phase diagram, their Fig. 6 included free energy differences between phases along the P^* = 290 isobar; this can be used to back-calculate an uncertainty of about 0.01–0.015 in T^* over this part of the range, slightly smaller than ours. The HCP temperature in our simulations stability is, however, lower than the results of Schultz and Kofke.²⁸ 

There are a number of factors that could potentially explain the discrepancies between other results that are statistically mechanically complete and our results; however, we find that most of them do not affect the phase diagram. The major methodological difference between previous anharmonic results and our PME SIMR results is the treatment of the cutoff, in particular, our use of PME to account for long-range interactions, rather than extrapolation to infinite cutoff. The use of PME and the parameters used for PSCP were validated to the results of Stillinger et al.,³⁷ with which multiple groups have got highly consistent results.²⁸ Our PME parameters gave energies at the lattice minimum of FCC and HCP lattice minima within 0.0001E^* of the results of Stillinger and ΔE^* between the phases of 7.3 × 10⁻⁷E^*. The average potential energy was changed significantly less than the statistical uncertainty by a change in the PME cutoff parameters. The difference in energy between HCP and FCC at the lattice minimum was ΔE^* = −8.70 × 10⁻⁴. Increasing the Ewald direct space cutoff by 0.25σ changed this energy difference by 8.5 × 10⁻⁶, less than 1% of the total relative energy and significantly below the statistical uncertainty. Increasing the number of Fourier points by 50% changed the energy by 1.34 × 10⁻⁵, approximately 1.5% relative error. Additionally, these terms cancel—when both increases in accuracy of the calculation are used, and the total change in energy is 3.4 × 10⁻⁷, less than 0.5% relative error. This degree of relative bias in the energy calculation should be roughly invariant, or even decrease, throughout the simulation, as the energy becomes dominated by the repulsive term that is entirely short range at higher pressures. Using a 1200 atom system and the PME settings described above, the initial structures were minimized for volumes ranging from 182.58 nm³ to 547.71 nm³. The minimum enthalpies were then compared against the values in the study by Stillinger³⁷ and the coexistence pressure was found. The T^* = 0 coexistence volume agreed to within 1.50 × 10⁻⁴V^*, and the enthalpy at that coexistence point agreed to within 0.136E^*, which is within 0.2%. Our results were also statistically tested for consistency by extrapolating the coexistence line at high pressures (where it is roughly linear) down to the 0 K point, though we did not simulate in this high-pressure region. We found that our extrapolated coexistence line crosses the T^* = 0 line between P^* = 845 and P^* = 958, which is within the uncertainty of the literature results of P^* = 878.

There are several other system differences that could, in theory, be causing the lower temperature coexistence in our results. One such effect is system size. Schultz and Kofke found that the system size effects are almost entirely harmonic.²⁸ To study these effects, we added a harmonic ΔG correction, equal to the difference in the harmonic free energy of each phase between a system of 1200 atoms and a system of infinite cutoff and size, with data graciously provided by the authors Schultz and Kofke.²⁸ However, this correction term only shifted the phase diagram up in temperature by an average value of 0.002T^*.

The effect of system size and cutoff was also tested at the T^* = 0 K line. A series of PME cutoffs between 3.2σ and 5.3σ was tested, and the T^* = 0 coexistence crossing was identified. The coexistence pressure was shifted by a value of 5P^* between the smallest and largest cutoffs; using the cutoff value of 3.2σ, the coexistence pressure is 878.26P^*, and that of 3.6σ, the pressure is 883.26P^*. Changing the system size from 1200 to 9600 atoms produced a change in energy at coexistence of 7.27 × 10⁻⁶E^*, and moving from the smallest to largest cutoff produced a change of 7.61 × 10⁻⁶E^*. Moving from a system size of 1200–9600 atoms produced a shift in the coexistence of 7P^*, which we considered to be sufficiently accurate, given that our main focus is the T^* > 0 portion of the phase boundary.

Another difference from the result of Schultz et al. is that they considered anisotropic expansion of the HCP box, which slightly shifts the phase coexistence curve in the direction of the HCP phase. However, their results (in Table III of Schultz and Kofke²⁸) show that this effect changes the location of the FCC-HCP-vapor triple point by ΔT^* = 0.0004 and the location of the T^* = 0 intersection of the coexistence curve by ΔP^* = 0.01 and is thus several orders of magnitude below the other differences considered in this paper.

After analyzing these other effects, we find that the maximum HCP stability temperature is almost entirely determined by the reference value obtained from the PSCP. If the PSCP reference value was shifted down by 0.0007E^*, our results would be within the uncertainty of the results of Schultz et al., with our maximum HCP stability raised to T^* = 0.402 compared with 0.40 in that study. Alternatively, if we use the coexistence point of Schultz and Kofke²⁸ of P^* = 127, T^* = 0.4 (obtained using the correlations found in Table I), as a reference, instead of the PSCP calculation, the maximum HCP stability moves to 0.42(2) and the T^* = 0 coexistence line still intersects between 853 and 1238 P^*, which are both within the uncertainty of literature results. The P^* = 0 intersection could not be calculated due to poor overlap in the low-pressure region with the simulated data set. Though this error in the reference ΔG generated using PSCP is small, it is several times the uncertainty in the calculation, indicating that some source of bias exists. This bias could be caused by several factors, including the lack of anisotropy in the PSCP simulations and finite-size effects in that calculation or very minor errors in the Parrinello-Rahman barostat implemented in GROMACS, which have been noted earlier.⁶³ We emphasize that any potential errors are so small so as to likely only be noticeable in calculations of this precision.

Other methods could be used to generate the PSCP, for example, lattice-switch Monte Carlo at a single (T^*, P^*) point, though the PSCP method is the most general for arbitrary crystal packings. In fact, the moves used in lattice-switch Monte Carlo themselves define a configuration mapping between the two phases, and the free energy between two phases can be determined directly from two simulations using Eq. (6), without ever actually implementing lattice-switch Monte Carlo. However, preliminary testing demonstrated that this approach was far too inaccurate for the LJ phase diagram because of the low overlap between the mapped FCC → HCP and HCP → FCC ensembles. This low overlap is not surprising given that lattice-switch Monte Carlo requires additional acceleration methods to yield reasonable results.³⁶ 

One challenge with studying the FCC and HCP structures using PME is the extremely small free energy differences between the polymorphs, as seen in Fig. 6, almost two orders of magnitude lower than the cutoff simulations. These small energy differences mean that longer (2.5× the potential switch simulations) and more closely spaced simulations are required to obtain sufficiently precise results. At pressures below 100P^*, the difference in free energy was small enough that sufficient precision to clearly resolve the phase diagram could not be achieved. This was due to consistently poor phase space overlap at these states, where the larger amount of movement in the atoms causes the mapping to work poorly.

This configuration mapping approach is likely to be useful for systems of point particles, whatever the potential may be, as the same mapping presented here will be valid. Based on previous results with configuration mapping on rigid water molecules in liquid water,⁴⁷ this approach is likely to work with small rigid molecules as well, but it is not clear how well it would work for more complex molecules with internal degrees of freedom that will change configuration ensemble between phases. There are a number of improvements that can be made as well. For example, adaptive choices of the simulation points, as discussed in Schieber et al.,¹³ can further decrease the clock time required at the expense of wall time. Many of the state-to-state mappings still lead to essentially negligible overlap and hence may be unnecessary; it may be possible to not map those states together, saving some amount of time spent mapping while losing negligible efficiency, though determining exactly which states to exclude or include is somewhat complicated.

The phase diagram of the solid phases of Lennard-Jones spheres has been predicted many times in the literature, with large discrepancies between methods. The Successive Interpolation of Multistate Reweighting (SIMR) is another method that can be used to determine the coexistence line of solid-solid transitions using full molecular dynamics simulations. However, this method is dependent on phase space overlap between adjacent simulations. Configuration mapping is a way to increase phase space overlap and therefore computational efficiency in phase diagram prediction.

Using a potential switch van der Waals cutoff, the coexistence curves were efficiently generated using this method to reasonable precision compared with the scale of Δμ between phases, and this demonstrates the efficiency of SIMR plus configuration mapping in a standard problem for point particles, such as that found in simulations of metals or inorganic materials. However, this cutoff approach introduces nonphysical behavior in the energy difference vs pressure curves, which can be understood by examining the radial distribution functions. As the pressure increases, more layers of atoms, and thus peaks in the radial distribution function (RDF), are brought under the value of the cutoff, which nonphysically affects the Gibbs free energy difference and thus the predicted coexistence line. Using particle mesh Ewald to calculate long-range interactions avoids nonphysical behavior in the energy differences between polymorphs, without extremely long cutoffs. The extremely small difference in potential energy and Gibbs free energy between phases and a small difference in slope between the free energy surfaces make the determination of this more accurate coexistence line challenging. In particular, at low pressures, the increased movement of the atoms decreases the effectiveness of the mapping, making overlap poor. Current limitations in analyzing larger numbers of states with MBAR prevents fully characterizing the phase diagram at lower pressures, though the uncertainty and bias are already very low (at or below 0.001E^*).

When using particle mesh Ewald, the phase diagram produced by the addition of configuration mapping to the SIMR method shows independence of cutoff and, despite the issues with statistical convergence, consistent trends in agreement with the current somewhat diverging literature results. The HCP phase is most stable at moderate pressures and low temperatures, with the FCC phase being more stable at high temperatures and extreme pressures. The coexistence curve displays reentrant behavior consistent with previous results. The maximum temperature of HCP stability is lower than the most likely most accurate results of Schultz et al.,¹⁵ likely due to the bias in the PSCP value used to generate our coexistence curve. We found that a change in the reference value on the order of 0.0007E^* brings out the curve to within the uncertainty of these literature results. These results demonstrate that the method, although not perfect for the present calculation, should be effective for most problems of applicable interest.

Included supplementary material consists of a discussion on the changes in free energy due to cutoffs, including analysis of the change in radial distribution functions as a function of pressure, and a comparison between uncertainty as determined by bootstrapping and the propagation of error. We also include the GROMACS input files used for simulations of Lennard-Jones particles for reference.

This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation Grant No. OCI-1053575. Specifically, it used the Bridges system, which is supported by the NSF Award No. ACI-1445606, at the Pittsburgh Supercomputing Center (PSC). This work was supported financially by the NSF through the Grant Nos. NSF-CBET 1351635 and NSF-DGE 1144083. We thank Nate Abraham (CU Boulder) for discussions and comparisons and Andrew Schultz (SUNY Buffalo) for helpful assistance in comparison with the results of Schultz et al.