- The paper implements nonequilibrium simulation techniques in LAMMPS to compute free energies of solids, reducing reliance on multiple equilibrium simulations.
- It combines Jarzynski's equality with forward-backward process results to mitigate statistical uncertainties and systematic errors.
- Application to iron polymorphism highlights phase stability challenges and underscores the need for validating interatomic potential models.
Free-Energy Calculation of Solids Using LAMMPS
This paper details the implementation of nonequilibrium techniques within the LAMMPS MD code for computing the free energies of solid materials. It emphasizes the efficiency gains achieved through these methods compared to traditional equilibrium thermodynamic integration approaches. The paper focuses on two nonequilibrium processes: calculating free-energy differences between systems with distinct Hamiltonians and determining the temperature dependence of a system's free energy from a single nonequilibrium simulation.
Theoretical Background
Traditional equilibrium free-energy calculations rely on thermodynamic integration, which involves constructing a series of equilibrium states along a path connecting two states of interest. The free-energy difference is then computed by averaging the thermodynamic driving force across these states through independent equilibrium simulations. Nonequilibrium approaches, however, utilize an explicitly time-dependent process, where the rate of execution determines the deviation from a quasistatic path. Jarzynski's equality provides the theoretical foundation for relating the work done in these nonequilibrium processes to the free-energy difference. While Jarzynski's equality is exact, its practical application can be limited by statistical uncertainties arising from the exponential average. An alternative method connects the mean value of the irreversible work distribution to the reversible work, reducing statistical uncertainties but introducing a systematic error in the form of dissipated heat. However, under linear response theory, this systematic error becomes the same for processes carried out in opposite directions, allowing its elimination by combining forward and backward process results. The paper underscores that, in contrast to equilibrium TI approaches, the nonequilibrium method samples the entire process during a single simulation, with the closeness to equilibrium assessed by monitoring the convergence as a function of process rate.
Implementation in LAMMPS
The paper outlines the implementation of the aforementioned nonequilibrium techniques within the LAMMPS MD code, demonstrating their application to calculating the temperature-dependent free energy of bcc iron using an EAM potential. The FL path is used to compute the reference free energy at a specific temperature, while the RS method extends the calculations across a broader temperature range. Details are provided regarding the preparation steps for the FL path, including timestep selection, thermostatting using the Langevin thermostat, and calculating the equilibrium volume using a barostat. The paper discusses the determination of the optimal spring constant for the harmonic reference system in the FL path, achieved by measuring the mean-squared displacement of atoms and applying the equipartition theorem. The implementation of the FL path in LAMMPS is facilitated by the fix ti/spring command, which performs the time-dependent switching between the EAM iron potential and the Einstein crystal Hamiltonian. The paper highlights the importance of accounting for the fixed center of mass constraint in the MD simulations by adding a correction term to the free-energy calculation. For the RS method, the scaling of the interatomic potential in LAMMPS is achieved using the fix adapt command, allowing for the time-dependent variation of the scaling parameter λ.
Application to Iron Polymorphism
The paper showcases the application of the nonequilibrium formalism to the study of polymorphism in iron, specifically examining the free energies of bcc, fcc, and hcp crystal structures using the EAM potential. System-size effects are carefully considered, with the analysis revealing that the free energy per atom converges with a leading term of $1/N$ after accounting for the center of mass term. The paper highlights the determination of minimum system sizes for each crystal structure to ensure convergence of the solid free energies to a desired accuracy. The results indicate that the EAM potential predicts a $#1{bcc} \rightarrow #1{fcc}$ transition temperature significantly lower than the experimental value and incorrectly predicts hcp as more stable than fcc across the analyzed temperature range. This underscores the importance of verifying the phase stability predicted by interatomic potentials before applying them to simulations of complex phenomena.
Conclusion
The paper successfully demonstrates the implementation and application of nonequilibrium simulation methods for calculating free energies of solids within LAMMPS. The efficiency and accuracy of the approach are highlighted through calculations on different polymorphs of iron. The paper concludes by discussing potential extensions of the technique to study more complex systems, including A-B compounds, alloy solid solutions, and systems with point and extended defects such as surfaces and interfaces. The authors also propose the use of this methodology for automating free-energy calculations, aiding in the development and benchmarking of classical interatomic potential models.