Nonequilibrium free-energy calculation of solids using LAMMPS

Rodrigo Freitas; Mark Asta; Maurice de Koning

Source-linked AI summary

Nonequilibrium free-energy calculation of solids using LAMMPS

Rodrigo Freitas, Mark Asta, Maurice de Koning

arXiv:2201.04622v1Jan 12, 2022cond-mat.mtrl-sci

TL;DR

The paper addresses the efficiency of computing free energies for crystalline solids compared with standard equilibrium thermodynamic integration. It implements two nonequilibrium routes in LAMMPS, applies them to iron polymorphs, and reports accurate estimates at practical system sizes and switching times.

Problem
Equilibrium free-energy calculations require path discretization and separate simulations for each state, motivating more efficient approaches for solids.
Method
The paper implements nonequilibrium Frenkel–Ladd and Reversible Scaling paths in LAMMPS for Hamiltonian free-energy differences and temperature-dependent free energies.
Results
Tenths of meV/atom or better were achieved for systems containing about 10,000 atoms using switching times on the order of 10^6 MD steps.
Takeaways & Limitations
The implementation supports efficient free-energy calculations and testing of phase stability for crystalline structures and interatomic potential models.

Abstract

from arXiv · show

This article describes nonequilibrium techniques for the calculation of free energies of solids using molecular dynamics (MD) simulations. These methods provide an alternative to standard equilibrium thermodynamic integration methods and often present superior efficiency. Here we describe the implementation in the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) code of two specific nonequilibrium processes that allow the calculation of the free-energy difference between two different system Hamiltonians as well as the free-energy temperature dependence of a given Hamiltonian, respectively. The theory behind the methods is summarized, and we describe (including fragments of LAMMPS scripts) how the process parameters should be selected to obtain the best-possible efficiency in the calculations of free energies using nonequilibrium MD simulations. As an example of the application of the methods we present results related to polymorphic transitions for a classical potential model of iron.

I. INTRODUCTION

The paper presents nonequilibrium free-energy methods implemented in LAMMPS as an efficient alternative to equilibrium thermodynamic integration. It develops the theoretical basis, practical considerations, and extensions for solids and related materials calculations.

Nonequilibrium processes provide an alternative to standard equilibrium free-energy calculations and can yield significant efficiency gains.Jarzynski’s equality connects nonequilibrium work distributions with equilibrium free-energy differences.
The implementation targets two routes: free-energy differences between different Hamiltonians and temperature dependence for a fixed Hamiltonian.The methods are implemented in the highly optimized LAMMPS molecular dynamics code.
The paper applies the methods to crystalline structures and discusses extensions to interfacial free energies and more complex crystalline systems.The implementation is presented with theoretical background, LAMMPS implementation guidance, and applications to iron polymorphism.
The exponential average in Jarzynski’s equality can produce substantial statistical uncertainties, motivating alternative estimators based on mean irreversible work.Forward and backward processes can eliminate systematic error when linear response theory applies.
Unlike equilibrium thermodynamic integration, nonequilibrium sampling can use a single simulation while convergence is monitored as the process rate changes.Equilibrium TI requires discretizing a path and running separate simulations for each state.

A. Nonequilibrium free-energy differences for a parameter-dependent Hamiltonian

For a parameter-dependent Hamiltonian, the free-energy difference between endpoint states can be estimated from a time-dependent nonequilibrium switching process. The method replaces pathwise equilibrium sampling with work accumulated along a trajectory, while requiring initial equilibration and control of dissipative error.

The target is ΔF = F(N,V,T;λf) − F(N,V,T;λi) between equilibrium states defined by two parameter values.
Equilibrium thermodynamic integration obtains this difference by integrating canonical averages over λ along a discretized path.Each grid point requires a separate equilibration and sampling simulation.
Nonequilibrium estimation varies λ explicitly with time from λi to λf during a switching time ts and evaluates the work integral along one trajectory.The continuous integral is approximated by an MD time-step sum using Δt and Δλk.
Dissipated heat introduces systematic error into averaged irreversible work, but forward and backward switching can remove it when the process is sufficiently slow for linear response theory.
The system must be equilibrated at the initial endpoint before the switching process begins.The required endpoint is λi or λf depending on the switching direction.

B. Free-energy difference between two systems: Frenkel–Ladd path

The Frenkel–Ladd path computes the free energy of an atomic solid by switching between its interaction Hamiltonian and an analytically tractable Einstein-crystal reference. Combining forward and backward switching estimates the reversible work and the solid’s free energy.

The interpolating Hamiltonian connects two different system Hamiltonians by setting λi = 0 and λf = 1.
The Frenkel–Ladd path uses the system’s interaction Hamiltonian as the initial state for calculating an atomic solid’s absolute free energy.
The reference state is an Einstein crystal of noninteracting particles attached to the equilibrium lattice by three-dimensional harmonic springs.Its Helmholtz free energy is known analytically.
The free energy of interest is estimated from the reversible work obtained by combining forward and backward switching processes.The two processes run between λ = 0 and λ = 1 in opposite directions.

C. Temperature dependence of the free energy: the Reversible Scaling path

The Reversible Scaling path maps a parameter-scaled Hamiltonian at reference temperature T0 onto the temperature dependence of a fixed Hamiltonian. Using irreversible work from one nonequilibrium simulation, it yields free energies across the mapped temperature interval.

Reversible Scaling assigns each value of λ to a temperature T = T0/λ for the original Hamiltonian H0.
The scaled Hamiltonian is constructed by introducing λ into H0’s potential-energy function, preserving a correspondence between configurational partition functions.
A single constant-temperature simulation varies λ from 1 to λf and estimates the irreversible work along that nonequilibrium path.
F0(T) can be calculated for all temperatures between T0 and T0/λf using the irreversible work from that simulation.
Applying the temperature-dependence relation requires prior knowledge of the free energy at the reference temperature T0.

III. CALCULATION OF BULK FREE ENERGY

The LAMMPS implementation applies nonequilibrium free-energy methods to a body-centered-cubic iron solid described by an EAM potential. The calculation combines a Frenkel–Ladd reference at 100 K with Reversible Scaling up to 1600 K at zero pressure.

LAMMPS computes temperature-dependent free energies for bcc iron using a Meyer–Entel EAM many-body potential.
The Frenkel–Ladd path establishes the reference free energy at T0 = 100 K, while Reversible Scaling extends the calculation to T = 1600 K.
The simulations calculate free energies at zero pressure over temperatures from 100 to 1600 K using the nonequilibrium approach.

A. Preparation for Frenkel–Ladd path

Preparation for the Frenkel–Ladd path sets stable simulation conditions, determines the equilibrium volume and an efficiency-oriented spring constant, and prevents center-of-mass thermostatting artifacts.

A timestep of Δt = 1 fs was selected after NVE energy-conservation benchmarks, and subsequent simulations use a thermostat.
The Frenkel–Ladd calculation uses fixed volume because the independent-particle Einstein crystal does not provide pressure through the virial pressure equation.
The equilibrium lattice parameter is first obtained at T0 = 100 K and zero hydrostatic pressure using a barostat.
The spring constant is chosen to match the Einstein crystal’s vibrational frequencies as closely as possible to the solid’s characteristic spectrum.
The spring constant must be measured at the same temperature and system size as the Frenkel–Ladd simulation because the phonon spectrum depends on both.
For 11,664 atoms at 100 K, the reported values are a(T0) = 2.8841(1) Å and k(T0) = 5.787(1) eV/Å2.
The simulations avoid the flying ice cube problem by measuring internal temperature, fixing the center of mass, and preventing thermostat-induced drift.

B. Frenkel–Ladd path in LAMMPS

The LAMMPS Frenkel–Ladd implementation switches between the EAM iron Hamiltonian and an Einstein crystal in forward and backward directions. It includes equilibration, configurable switching functions, and a finite-size correction for the fixed center of mass.

The fix ti/spring performs time-dependent switching between the EAM iron potential and the Einstein crystal along the Frenkel–Ladd path.
The fix syntax specifies the spring constant, full switching time, equilibration time, and switching function.
The protocol equilibrates under H0, switches λ from 0 to 1, equilibrates under HE, and switches λ back to 0.
A switching function with vanishing endpoint slope is implemented because such forms have been shown to produce less dissipative switching.
The fixed-center-of-mass constraint requires a first-order free-energy correction that vanishes per particle as ln N/N in the thermodynamic limit.
Forward and backward switches used 11,664 atoms at 100 K, k = 5.787 eV/Å2, teq = 0.1 ns, and ten independent realizations.

C. Reversible Scaling in LAMMPS

Reversible Scaling extends a reference free energy at T0 to a temperature-dependent free-energy curve by scaling the potential energy in LAMMPS. The implementation and parameter choices support efficient nonequilibrium calculations while controlling dissipation through switching-direction averaging.

Reference state: The RS path uses G(T0) at zero pressure as its reference point and extends the free energy to higher temperatures.The reference temperature is chosen at low temperature because Frenkel–Ladd calculations converge rapidly there.
Path implementation: During RS, the Hamiltonian uses a scaled potential energy λU(r), with λ varied from λi = 1 to λf.LAMMPS changes the interatomic potential directly using fix adapt while maintaining temperature and pressure with a thermostat and barostat.
Path implementation: LAMMPS implements potential scaling with fix adapt and a variable λ that multiplies all EAM pair interactions.The command is invoked every timestep and applies the scale factor to all atom-type pairs.
Parameter selection: The selected λ(τ) function produces a constant rate of temperature change and gives quickly converging RS calculations.The paper notes that narrower temperature ranges can reduce dissipation and improve precision.
Results: A switch from T0 = 100 K to 1600 K with λf = 0.0625 reproduced the free-energy curve with 1 meV/atom precision using 4 × 10^4 MD steps.At the endpoint, where cumulative dissipation is largest, averaging forward and backward simulations eliminates the systematic dissipation error.
Validation: Frenkel–Ladd checks at 400, 700, 1000, 1300, and 1600 K showed excellent agreement with the RS temperature-dependence curve.The RS and FL paths were compared using a longer switching time of ts = 2.0 ns.

IV. IRON POLYMORPHISM

Polymorphism is the temperature-dependent thermodynamic stability of crystal structures, making phase stability a central issue in modeling crystalline solids. Iron provides an example involving multiple crystal structures.

Definition: Polymorphism is the thermodynamic stability of phases with different crystal structures as a function of temperature.The paper identifies polymorphism as a central issue in modeling crystalline solids and notes that iron is an example.

A. Free-energy size dependence

The study examines finite-size effects in vibrational free energies of crystalline solids. After correcting for the fixed center of mass, the free energy per atom converges with leading term 1/N.

Finite-size effects: The calculations target vibrational free energy, with finite-size effects arising from the phonon-spectrum cutoff imposed by finite lattice size.A center-of-mass contribution is treated separately.
Finite-size effects: The bcc size-dependence study used cubic boxes from M = 1 to 30, including a largest system of N = 54,000 atoms.The calculations were performed for the bcc phase at 1600 K.
Convergence: After the fixed-center-of-mass correction, the free energy per atom converges with leading term 1/N.This behavior agrees with analytical finite-size calculations.
Implementation: The implementation includes a barostat setup using fix nph with anisotropic pressure control.The supplied script fragment specifies the corresponding LAMMPS command.

B. Results and discussion

The nonequilibrium calculations reproduce phase free-energy differences for an iron potential while revealing incorrect high-temperature phase stability. They also achieve high precision for systems of roughly 10,000 atoms.

Method: The same free-energy procedure was applied to bcc, fcc, and hcp iron structures.System sizes were approximately 10,000–11,700 atoms, with T0 = 100 K and switching time ts = 2 ns.
bcc–fcc transition: The predicted bcc → fcc transition occurs at 487 K, almost 2.5 times below the experimental value of 1183 K.Nonequilibrium results agree within error bars with independent equilibrium thermodynamic-integration results.
Finite-size convergence: At 1600 K, the free-energy convergence study targets accuracy within ±0.2 meV/atom of the thermodynamic-limit value.The thermodynamic-limit estimate uses the leading 1/N convergence behavior.
Validation: The reported free-energy differences relative to the independent comparison range from ±0.5 meV/atom to ±3.4 meV/atom.These differences are smaller than or comparable to the comparison results’ error bars.
fcc–hcp stability: The EAM potential predicts hcp to be more stable than fcc across the analyzed temperature range.Because hcp is not thermodynamically stable at zero pressure in observed iron, the potential gives incorrect high-temperature phase stability above 440 K.
Implications: The LAMMPS implementation provides a framework for testing phase stability in other interatomic potential models exhibiting polymorphic transitions below melting.The paper positions this as an application of the formalism to other solid models.

V. SUMMARY AND DISCUSSION

The implemented Frenkel–Ladd and Reversible Scaling nonequilibrium paths compute solid free energies efficiently and accurately, while supporting extensions to complex materials and potential-model benchmarking. Applications to iron polymorphs also reveal that the examined potential predicts incorrect phase stability above 440 K.

Summary and discussion: Tenths of meV/atom precision is achievable for systems containing about 10,000 atoms or more with switching times around 10^6 MD steps.Switching times as short as 4×10^4 steps for Reversible Scaling or 300 steps for Frenkel–Ladd achieve 1 meV/atom accuracy.
Summary and discussion: The LAMMPS implementation can extend nonequilibrium free-energy calculations from simple solids to compounds, alloys, surfaces, interfaces, and other defects.The methods support vibrational free energies for A–B compounds or alloys and excess free energies for surfaces and interfaces.
Summary and discussion: Different spring constants for atomic species with strongly differing vibrational frequencies may improve Frenkel–Ladd accuracy without additional computational cost.The recommended constants can be selected by monitoring the relevant vibrational behavior during the calculation.
Summary and discussion: For a surface, the excess free energy is obtained from γ = (Fsurf − Fbulk)/A, and the same framework can be applied to interfaces such as grain boundaries.The procedure involves subtracting the bulk free energy from the surface-system free energy and dividing by surface area.
Summary and discussion: Automated free-energy calculations can support development and benchmarking of classical interatomic potentials by evaluating predicted bulk thermodynamics and phase stability.Such calculations can be performed for typical crystal structures before applying a potential to more complex phenomena.