这是用户在 2025-6-3 15:01 为 https://app.immersivetranslate.com/pdf-pro/d9fae50f-ac49-451e-8752-aa7536725749/ 保存的双语快照页面，由 沉浸式翻译 提供双语支持。了解如何保存？

Nonequilibrium free-energy calculation of solids using LAMMPS
使用 LAMMPS 计算固体的非平衡自由能

Rodrigo Freitas, 1 , 2 1 , 2 ^(1,2){ }^{1,2} Mark Asta, 1 1 ^(1){ }^{1} and Maurice de Koning 2 2 ^(2){ }^{2}
罗德里戈·弗雷塔斯、 1 , 2 1 , 2 ^(1,2){ }^{1,2} 马克·阿斯塔和 1 1 ^(1){ }^{1} 莫里斯·德·科宁 2 2 ^(2){ }^{2} 1 1 ^(1){ }^{1} Department of Materials Science and Engineering, University of California, Berkeley, CA 94720, USA
1 1 ^(1){ }^{1} 加州大学伯克利分校材料科学与工程系，美国 CA 94720 2 2 ^(2){ }^{2} Instituto de Física "Gleb Wataghin", Universidade Estadual de Campinas, UNICAMP, Campinas, São Paulo 13083-859, Brazil*
2 2 ^(2){ }^{2} Instituto de Física “Gleb Wataghin”， Universidade Estadual de Campinas， UNICAMP， 坎皮纳斯， 圣保罗 13083-859， 巴西*

(Dated: January 13, 2022)
（日期：2022 年 1 月 13 日）

Abstract  抽象

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.
本文介绍了使用分子动力学 （MD） 模拟计算固体自由能的非平衡技术。这些方法提供了标准平衡热力学积分方法的替代方案，并且通常具有卓越的效率。在这里，我们描述了在 LAMMPS（大规模原子/分子大规模并行模拟器）代码中两个特定非平衡过程的实现，这两个过程允许分别计算两个不同系统哈密顿量之间的自由能差以及给定哈密顿量的自由能温度依赖性。总结了这些方法背后的理论，并描述了（包括 LAMMPS 脚本的片段）应如何选择工艺参数，以便在使用非平衡 MD 模拟计算自由能时获得最佳效率。作为方法应用的一个例子，我们提出了与经典铁势模型的多态转变相关的结果。

I. INTRODUCTION  I. 引言

The calculation of free energies and derivative thermodynamic quantities for condensed-phase systems is a common and widespread application of atomistic simulation techniques. The continuous improvement of methods and the steady increase in computational power have enabled the efficient and reliable computation of free energies for complex systems of interest in materials science, including extended defects and interfaces. Freeenergy calculation methods have benefited substantially with the introduction of nonequilibrium approaches such as the adiabatic switching [1] (AS) method, which allow the calculations to be performed along explicitly timedependent processes that can lead to significant efficiency gains compared to standard equilibrium methods. Moreover, with the derivation of Jarzynski’s equality the connection between equilibrium free-energy differences and nonequilibrium processes has been placed on a firm theoretical basis [2].
凝聚相系统的自由能和导数热力学量的计算是原子仿真技术的常见且广泛的应用。方法的不断改进和计算能力的稳步提高，使得材料科学中感兴趣的复杂系统（包括扩展缺陷和界面）能够高效可靠地计算自由能。自由能计算方法随着非平衡方法的引入而受益匪浅，例如绝热开关 [1] （AS） 方法，该方法允许沿着显式的时间依赖过程执行计算，与标准平衡方法相比，这可以显着提高效率。此外，随着 Jarzynski 等式的推导，平衡自由能差异和非平衡过程之间的联系已经建立在坚实的理论基础上 [2]。
In this paper we describe the implementation of state-of-the-art nonequilibrium techniques for calculating free energies of solids in the highly optimized LAMMPS [3] (Large-scale Atomic/Molecular Massively Parallel Simulator) molecular dynamics code. Specifically, we focus on two particular kinds of nonequilibrium routes which, respectively, allow the calculation of free-energy differences between systems described by different Hamiltonians as well as the temperature dependence of the free energy of a given system Hamiltonian from a single nonequilibrium simulation. We demonstrate their application to the calculation of free energies of different crystalline structures for a classical interatomic potential model of iron [4] and discuss extensions of the approaches to the calculations of
在本文中，我们描述了在高度优化的 LAMMPS [3]（大规模原子/分子大规模并行模拟器）分子动力学代码中计算固体自由能的最先进的非平衡技术的实现。具体来说，我们专注于两种特殊类型的非平衡路线，它们分别允许计算由不同哈密顿量描述的系统之间的自由能差异，以及来自单个非平衡模拟的给定系统哈密顿量的自由能的温度依赖性。我们展示了它们在铁的经典原子间势模型 [4] 中不同晶体结构的自由能计算中的应用，并讨论了计算方法的扩展
interfacial free energies, which can be readily undertaken using the methods implemented in LAMMPS [5].
界面自由能，这可以使用 LAMMPS 中实现的方法轻松实现 [5]。
The paper has been organized as follows. In Sec. II we give a short and self-contained presentation of the theoretical framework underlying the nonequilibrium processes used to compute free-energy differences between equilibrium states. In Sec. III we describe the implementation of these methods into LAMMPS, and discuss how to optimize the parameters for efficiency, using bcc iron as an example system. In Sec. IV we describe the application of the methods to the study of polymorphism in iron and we end with a summary and a discussion of the extension of the methodology to more complex crystalline systems as well as surfaces and interfaces in Sec. V.
本文的组织结构如下。在第二节中，我们简短而独立地介绍了用于计算平衡状态之间自由能差异的非平衡过程的基础理论框架。在第 III 节中，我们描述了这些方法在 LAMMPS 中的实施，并讨论了如何以 bcc 铁作为示例系统优化参数以提高效率。在第四节中，我们描述了这些方法在铁多晶型研究中的应用，最后我们总结并讨论了该方法在第五节中扩展到更复杂的晶体系统以及表面和界面。

II. NONEQUILIBRIUM FREE-ENERGY ESTIMATION
II. 非平衡自由能估计

Standard equilibrium free-energy calculations are often based on thermodynamic integration (TI), [6,7] which is a general class of methods based on the construction of a sequence of equilibrium states on a path between two thermodynamic states of interest. The free-energy difference between two equilibrium states is then determined by computing ensemble averages of the relevant thermodynamic driving force for these states by means of a set of independent equilibrium simulations, followed by numerical integration. This approach embodies the thermodynamic equality between the free-energy difference between two equilibrium states and the reversible work W rev W rev  W_("rev ")W_{\text {rev }} that is done along a quasistatic path that connects them.
标准平衡自由能计算通常基于热力学积分 （TI） [6,7]，这是一类通用的方法，基于在两个感兴趣的热力学状态之间的路径上构建一系列平衡态。然后，通过一组独立的平衡模拟计算这些状态的相关热力学驱动力的集合平均值，然后进行数值积分，来确定两个平衡状态之间的自由能差。这种方法体现了两个平衡态之间的自由能差与沿连接它们的准静态路径完成的可逆功 W rev W rev  W_("rev ")W_{\text {rev }} 之间的热力学相等性。
Nonequilibrium approaches, on the other hand, envision this path in terms of an explicitly time-dependent process. The rate at which this process is executed then determines the amount by which it deviates from a quasistatic path. Jarzynski’s equality [2] relates the work
另一方面，非平衡方法根据明确的时间依赖过程来设想这条路径。然后，此过程的执行速率决定了它偏离准静态路径的量。Jarzynski 的等式 [2] 与这项工作相关
done along arbitrarily out-of-equilibrium processes, starting from an equilibrated initial state, to the free-energy difference Δ F Δ F Delta F\Delta F between this state and the equilibrium state at the end of the process. In contrast to a quasistatic path, however, the work done in a nonequilibrium process is a stochastic variable that differs for each realization and simulations provide a means for sampling its distribution function. Jarzynski’s equality then connects this distribution to the free-energy difference Δ F Δ F Delta F\Delta F between the equilibrium states defined by the parameters at the two ends of the path through
沿着任意的不平衡过程完成，从平衡的初始状态开始，到该过程结束时该状态与平衡状态之间的自由能差 Δ F Δ F Delta F\Delta F 。然而，与准静态路径相反，在非平衡过程中完成的工作是一个随机变量，每个实现都不同，模拟提供了一种对其分布函数进行采样的方法。然后，Jarzynski 等式将此分布与平衡状态之间的自由能差 Δ F Δ F Delta F\Delta F 联系起来，该平衡状态由路径两端的参数定义
exp ⁡ ( − β Δ F ) = exp ⁡ ( − β W irr ) ― exp ⁡ ( − β Δ F ) = exp ⁡ − β W irr ¯ exp(-beta Delta F)= bar(exp(-betaW_(irr)))\exp (-\beta \Delta F)=\overline{\exp \left(-\beta W_{\mathrm{irr}}\right)}
where W irr W irr  W_("irr ")W_{\text {irr }} is the irreversible work done along such a nonequilibrium process and the overline denotes averaging over an ensemble of different realizations. While the equality is exact, its practical application in free-energy calculations is often limited due the exponential average in Eq. (1), which may lead to very substantial statistical uncertainties in its evaluation [8].
其中 W irr W irr  W_("irr ")W_{\text {irr }} 是沿着这种非平衡过程完成的不可逆功，上划线表示在不同实现的集合上的平均。虽然相等式是精确的，但由于方程（1）中的指数平均值，它在自由能计算中的实际应用往往受到限制，这可能导致其评估中存在非常大的统计不确定性[8]。
Another way to relate Δ F Δ F Delta F\Delta F to the irreversible work distribution is by directly connecting its mean value to the true reversible work. One has
与不可逆功分布相关的 Δ F Δ F Delta F\Delta F 另一种方法是将其平均值直接连接到真正的可逆功。一个有
Δ F = W rev = W irr ― − E diss ― Δ F = W rev = W irr ¯ − E diss ¯ Delta F=W_(rev)= bar(W_(irr))- bar(E_(diss))\Delta F=W_{\mathrm{rev}}=\overline{W_{\mathrm{irr}}}-\overline{E_{\mathrm{diss}}}
where E diss ― E diss  ¯ bar(E_("diss "))\overline{E_{\text {diss }}} is the average dissipated heat generated for an ensemble of replicas of the nonequilibrium process. Because of the second law of thermodynamics we have E diss ≥ 0 E diss  ≥ 0 E_("diss ") >= 0E_{\text {diss }} \geq 0, with the equality being valid only in the limit of an infinitely slow, quasistatic process. Instead of the exponential average in Eq. (1), this relation involves a simple mean of W irr W irr  W_("irr ")W_{\text {irr }} values. This significantly reduces the statistical uncertainties but the price to pay is the presence of the a priori unknown systematic error in the form of the dissipated heat E diss ― E diss  ¯ bar(E_("diss "))\overline{E_{\text {diss }}}.
其中 E diss ― E diss  ¯ bar(E_("diss "))\overline{E_{\text {diss }}} 是非平衡过程的一组复制品产生的平均耗散热量。由于热力学第二定律，我们有 E diss ≥ 0 E diss  ≥ 0 E_("diss ") >= 0E_{\text {diss }} \geq 0 ，相等性仅在无限慢的准静态过程的极限下有效。与方程 （1） 中的指数平均值不同，这种关系涉及一个简单的 W irr W irr  W_("irr ")W_{\text {irr }} 值平均值。这大大降低了统计的不确定性，但付出的代价是存在以散热形式存在的先验未知系统误差 E diss ― E diss  ¯ bar(E_("diss "))\overline{E_{\text {diss }}} 。
However, provided that the nonequilibrium process is sufficiently “close” to the ideally quasistatic process for linear response theory to be valid, it can be shown that this systematic error becomes the same for two processes that are carried out in opposite directions [9]. In other words, for a linear-response nonequilibrium process connecting states 1 and 2 , we have E 1 → 2 diss ― = E 2 → 1 diss ― E 1 → 2 diss  ¯ = E 2 → 1 diss  ¯ bar(E_(1rarr2)^("diss "))= bar(E_(2rarr1)^("diss "))\overline{E_{1 \rightarrow 2}^{\text {diss }}}=\overline{E_{2 \rightarrow 1}^{\text {diss }}}, from which it follows
然而，如果非平衡过程足够“接近”理想准静态过程，线性响应理论才有效，就可以证明这种系统误差对于以相反方向进行的两个过程是相同的 [9]。换句话说，对于连接状态 1 和 2 的线性响应非平衡过程，我们有 E 1 → 2 diss ― = E 2 → 1 diss ― E 1 → 2 diss  ¯ = E 2 → 1 diss  ¯ bar(E_(1rarr2)^("diss "))= bar(E_(2rarr1)^("diss "))\overline{E_{1 \rightarrow 2}^{\text {diss }}}=\overline{E_{2 \rightarrow 1}^{\text {diss }}} ，它由此得出
Δ F ≡ F 2 − F 1 ≡ 1 2 [ W 1 → 2 rev − W 2 → 1 rev ] = 1 2 { [ W 1 → 2 irr ― − E 1 → 2 diss ― ] − [ W 2 → 1 irr ― − E 1 → 2 diss ― ] } = 1 2 [ W 1 → 2 irr ― − W 2 → 1 irr ― ] Δ F ≡ F 2 − F 1 ≡ 1 2 W 1 → 2 rev  − W 2 → 1 rev  = 1 2 W 1 → 2 irr  ¯ − E 1 → 2 diss  ¯ − W 2 → 1 irr  ¯ − E 1 → 2 diss  ¯ = 1 2 W 1 → 2 irr  ¯ − W 2 → 1 irr  ¯ {:[Delta F-=F_(2)-F_(1)],[-=(1)/(2)[W_(1rarr2)^("rev ")-W_(2rarr1)^("rev ")]],[=(1)/(2){[ bar(W_(1rarr2)^("irr "))- bar(E_(1rarr2)^("diss "))]-[ bar(W_(2rarr1)^("irr "))- bar(E_(1rarr2)^("diss "))]}],[=(1)/(2)[ bar(W_(1rarr2)^("irr "))- bar(W_(2rarr1)^("irr "))]]:}\begin{aligned} \Delta F & \equiv F_{2}-F_{1} \\ & \equiv \frac{1}{2}\left[W_{1 \rightarrow 2}^{\text {rev }}-W_{2 \rightarrow 1}^{\text {rev }}\right] \\ & =\frac{1}{2}\left\{\left[\overline{W_{1 \rightarrow 2}^{\text {irr }}}-\overline{E_{1 \rightarrow 2}^{\text {diss }}}\right]-\left[\overline{W_{2 \rightarrow 1}^{\text {irr }}}-\overline{E_{1 \rightarrow 2}^{\text {diss }}}\right]\right\} \\ & =\frac{1}{2}\left[\overline{W_{1 \rightarrow 2}^{\text {irr }}}-\overline{W_{2 \rightarrow 1}^{\text {irr }}}\right] \end{aligned}
eliminating the systematic error by combining the results of the forward and backward processes. Similarly, the magnitude of the dissipation in this regime can be estimated by
通过组合前进和后向过程的结果来消除系统误差。同样，这种状态下的耗散大小可以通过以下方式估算
E 1 → 2 diss ― = E 2 → 1 diss ― = 1 2 [ W 1 → 2 irr ― + W 2 → 1 irr ― ] E 1 → 2 diss  ¯ = E 2 → 1 diss  ¯ = 1 2 W 1 → 2 irr  ¯ + W 2 → 1 irr  ¯ bar(E_(1rarr2)^("diss "))= bar(E_(2rarr1)^("diss "))=(1)/(2)[ bar(W_(1rarr2)^("irr "))+ bar(W_(2rarr1)^("irr "))]\overline{E_{1 \rightarrow 2}^{\text {diss }}}=\overline{E_{2 \rightarrow 1}^{\text {diss }}}=\frac{1}{2}\left[\overline{W_{1 \rightarrow 2}^{\text {irr }}}+\overline{W_{2 \rightarrow 1}^{\text {irr }}}\right]
One of the challenges of standard equilibrium TI approaches is how to discretize the quasistatic path between the two states of interest. In addition to deciding on the number of states, it also requires choosing how to distribute them over the path. Furthermore, every state requires a separate simulation, each of which should allow sufficient equilibration as well as ensemble averaging time.
标准平衡 TI 方法的挑战之一是如何离散两种感兴趣状态之间的准静态路径。除了决定状态的数量外，还需要选择如何在路径上分配它们。此外，每个状态都需要一个单独的仿真，每个仿真都应该允许足够的平衡以及集成平均时间。
In the nonequilibrium approach discussed above, on the other hand, the entire process is sampled during a single simulation and the closeness to equilibrium can be systematically assessed by monitoring the convergence of Eq. (3) as a function of process rate. The characteristic that a desired free-energy value can be estimated from a few relatively short simulations and that its convergence can be systematically verified render it an attractive alternative to the standard equilibrium TI methodology and often gives substantially improved efficiency.
另一方面，在上面讨论的非平衡方法中，整个过程在单个模拟中采样，并且可以通过监测方程 （3） 作为过程速率的函数的收敛性来系统地评估接近平衡的程度。可以从几个相对较短的模拟中估计出所需的自由能值，并且可以系统地验证其收敛性，这使其成为标准平衡 TI 方法的有吸引力的替代方案，并且通常可以显着提高效率。

A. Nonequilibrium free-energy differences for a parameter-dependent Hamiltonian
A. 参数依赖性哈密顿量的非平衡自由能差值

Consider a system of N N NN particles confined to a volume V V VV and in thermal equilibrium with a heat reservoir at temperature T T TT. The Hamiltonian of the system is given by H ( Γ , λ ) H ( Γ , λ ) H(Gamma,lambda)H(\boldsymbol{\Gamma}, \lambda), where Γ = ( r , p ) Γ = ( r , p ) Gamma=(r,p)\boldsymbol{\Gamma}=(\mathbf{r}, \mathbf{p}) is a point in the phase space of the particles in the system (atoms in the present context) and λ λ lambda\lambda is a parameter, a specific example of which will be given below. The Helmholtz free energy of this system, for a particular value of λ λ lambda\lambda, is F ( N , V , T ; λ ) = − k B T ln ⁡ Z ( N , V , T ; λ ) F ( N , V , T ; λ ) = − k B T ln ⁡ Z ( N , V , T ; λ ) F(N,V,T;lambda)=-k_(B)T ln Z(N,V,T;lambda)F(N, V, T ; \lambda)=-k_{\mathrm{B}} T \ln Z(N, V, T ; \lambda) where k B k B k_(B)k_{\mathrm{B}} is the Boltzmann constant and Z Z ZZ is the system’s canonical partition function given as an integral over the entire phasespace volume:
考虑一个 N N NN 粒子系统，该系统被限制在一个体积 V V VV 内，并在温度下 T T TT 有一个热储存器。系统的哈密顿量由 H ( Γ , λ ) H ( Γ , λ ) H(Gamma,lambda)H(\boldsymbol{\Gamma}, \lambda) 给出，其中 Γ = ( r , p ) Γ = ( r , p ) Gamma=(r,p)\boldsymbol{\Gamma}=(\mathbf{r}, \mathbf{p}) 是系统中粒子（当前上下文中的原子）相空间中的一个点， λ λ lambda\lambda 是一个参数，下面将给出一个具体示例。对于特定值 ， λ λ lambda\lambda 该系统的亥姆霍兹自由能 是 F ( N , V , T ; λ ) = − k B T ln ⁡ Z ( N , V , T ; λ ) F ( N , V , T ; λ ) = − k B T ln ⁡ Z ( N , V , T ; λ ) F(N,V,T;lambda)=-k_(B)T ln Z(N,V,T;lambda)F(N, V, T ; \lambda)=-k_{\mathrm{B}} T \ln Z(N, V, T ; \lambda) 中的 k B k B k_(B)k_{\mathrm{B}} 玻尔兹曼常数 ，是 Z Z ZZ 系统的规范分配函数，作为整个相空间体积的积分给出：
Z ( N , V , T ; λ ) = ∫ d Γ h 3 N exp ⁡ [ − β H ( Γ , λ ) ] Z ( N , V , T ; λ ) = ∫ d Γ h 3 N exp ⁡ [ − β H ( Γ , λ ) ] Z(N,V,T;lambda)=int(dGamma)/(h^(3N))exp[-beta H(Gamma,lambda)]Z(N, V, T ; \lambda)=\int \frac{\mathrm{d} \boldsymbol{\Gamma}}{h^{3 N}} \exp [-\beta H(\boldsymbol{\Gamma}, \lambda)]
where β = 1 / k B T β = 1 / k B T beta=1//k_(B)T\beta=1 / k_{\mathrm{B}} T, and h h hh is Planck’s constant.
其中 β = 1 / k B T β = 1 / k B T beta=1//k_(B)T\beta=1 / k_{\mathrm{B}} T ，和 h h hh 是普朗克常数。
Let us consider the problem of calculating the freeenergy difference between two thermodynamic states characterized by different values of the parameter λ λ lambda\lambda (namely λ i λ i lambda_(i)\lambda_{\mathrm{i}} and λ f λ f lambda_(f)\lambda_{\mathrm{f}} ), i.e., our goal is to calculate Δ F ( N , V , T ) ≡ F ( N , V , T ; λ f ) − F ( N , V , T ; λ i ) Δ F ( N , V , T ) ≡ F N , V , T ; λ f − F N , V , T ; λ i Delta F(N,V,T)-=F(N,V,T;lambda_(f))-F(N,V,T;lambda_(i))\Delta F(N, V, T) \equiv F\left(N, V, T ; \lambda_{\mathrm{f}}\right)-F\left(N, V, T ; \lambda_{i}\right). From now on we will omit the dependence on N , V N , V N,VN, V, and T T TT when their values are clear from the context.
让我们考虑计算两个热力学状态之间的自由能差的问题，这些状态以参数的不同值 λ λ lambda\lambda （即 λ i λ i lambda_(i)\lambda_{\mathrm{i}} 和 λ f λ f lambda_(f)\lambda_{\mathrm{f}} ）为特征，即我们的目标是计算 Δ F ( N , V , T ) ≡ F ( N , V , T ; λ f ) − F ( N , V , T ; λ i ) Δ F ( N , V , T ) ≡ F N , V , T ; λ f − F N , V , T ; λ i Delta F(N,V,T)-=F(N,V,T;lambda_(f))-F(N,V,T;lambda_(i))\Delta F(N, V, T) \equiv F\left(N, V, T ; \lambda_{\mathrm{f}}\right)-F\left(N, V, T ; \lambda_{i}\right) 。从现在开始，我们将省略对 N , V N , V N,VN, V 的依赖性，并且 T T TT 当它们的值从上下文中明确时。
The desired free-energy difference can be computed by finding the derivative of the F F FF with respect to λ λ lambda\lambda,
所需的自由能差值可以通过求 关于 λ λ lambda\lambda 的 F F FF 导数来计算 。
∂ F ∂ λ = 1 Z ∫ d Γ h 3 N ∂ H ∂ λ exp ⁡ [ − β H ( Γ , λ ) ] = ⟨ ∂ H ∂ λ ⟩ λ ∂ F ∂ λ = 1 Z ∫ d Γ h 3 N ∂ H ∂ λ exp ⁡ [ − β H ( Γ , λ ) ] = ∂ H ∂ λ λ (del F)/(del lambda)=(1)/(Z)int((d)Gamma)/(h^(3N))(del H)/(del lambda)exp[-beta H(Gamma,lambda)]=(:(del H)/(del lambda):)_(lambda)\frac{\partial F}{\partial \lambda}=\frac{1}{Z} \int \frac{\mathrm{~d} \boldsymbol{\Gamma}}{h^{3 N}} \frac{\partial H}{\partial \lambda} \exp [-\beta H(\boldsymbol{\Gamma}, \lambda)]=\left\langle\frac{\partial H}{\partial \lambda}\right\rangle_{\lambda}
where ⟨ … ⟩ λ ⟨ … ⟩ λ (:dots:)_(lambda)\langle\ldots\rangle_{\lambda} is the canonical ensemble average for a specific value of the parameter λ λ lambda\lambda, and integrate it to obtain
其中 ⟨ … ⟩ λ ⟨ … ⟩ λ (:dots:)_(lambda)\langle\ldots\rangle_{\lambda} 是参数 λ λ lambda\lambda 的特定值的规范集成平均值，并将其积分得到
Δ F = F ( λ f ) − F ( λ i ) = ∫ λ i λ f d λ ⟨ ∂ H ∂ λ ⟩ λ ≡ W i → f rev Δ F = F λ f − F λ i = ∫ λ i λ f   d λ ∂ H ∂ λ λ ≡ W i → f rev {:[Delta F=F(lambda_(f))-F(lambda_(i))=int_(lambda_(i))^(lambda_(f))dlambda(:(del H)/(del lambda):)_(lambda)],[-=W_(irarrf)^(rev)]:}\begin{aligned} \Delta F=F\left(\lambda_{\mathrm{f}}\right)-F\left(\lambda_{\mathrm{i}}\right) & =\int_{\lambda_{\mathrm{i}}}^{\lambda_{\mathrm{f}}} \mathrm{~d} \lambda\left\langle\frac{\partial H}{\partial \lambda}\right\rangle_{\lambda} \\ & \equiv W_{\mathrm{i} \rightarrow \mathrm{f}}^{\mathrm{rev}} \end{aligned}
The integral on the right-hand side can be interpreted as the reversible work along a quasistatic process between the equilibrium states with the two λ λ lambda\lambda-values of interest. In the equilibrium TI approach, this integral is discretized on a grid of λ λ lambda\lambda-values and for each value a separate equilibration simulation is executed.
右侧的积分可以解释为平衡态与两个 λ λ lambda\lambda 感兴趣的值之间沿准静态过程的可逆功。在平衡 TI 方法中，此积分在 λ λ lambda\lambda -值网格上离散化，并且对于每个值执行单独的平衡仿真。
In the nonequilibrium approach the integral in Eq. (7) is estimated in terms of the irreversible work done along a single simulation in which λ = λ ( t ) λ = λ ( t ) lambda=lambda(t)\lambda=\lambda(t) is explicitly time dependent and varied from λ i λ i lambda_(i)\lambda_{\mathrm{i}} to λ f λ f lambda_(f)\lambda_{\mathrm{f}} in a switching time t s t s t_(s)t_{\mathrm{s}},
在非平衡方法中，方程 （7） 中的积分是根据在单个模拟中完成的不可逆功来估计的，其中 λ = λ ( t ) λ = λ ( t ) lambda=lambda(t)\lambda=\lambda(t) 显式地依赖于时间，并且在切换时间 t s t s t_(s)t_{\mathrm{s}} 中从 λ i λ i lambda_(i)\lambda_{\mathrm{i}} 到 λ f λ f lambda_(f)\lambda_{\mathrm{f}} 变化。
W i → f irr = ∫ 0 t s d t d λ d t ( ∂ H ∂ λ ) Γ ( t ) W i → f irr = ∫ 0 t s   d t d λ d t ∂ H ∂ λ Γ ( t ) W_(irarrf)^(irr)=int_(0)^(t_(s))dt((d)lambda)/((d)t)((del H)/(del lambda))_(Gamma(t))W_{\mathrm{i} \rightarrow \mathrm{f}}^{\mathrm{irr}}=\int_{0}^{t_{\mathrm{s}}} \mathrm{~d} t \frac{\mathrm{~d} \lambda}{\mathrm{~d} t}\left(\frac{\partial H}{\partial \lambda}\right)_{\boldsymbol{\Gamma}(t)}
where Γ ( t ) Γ ( t ) Gamma(t)\boldsymbol{\Gamma}(t) represents the phase-space trajectory of the system along the process. In practice, the integral in Eq. (8) is evaluated in terms of the sum
其中 Γ ( t ) Γ ( t ) Gamma(t)\boldsymbol{\Gamma}(t) 表示系统沿过程的相空间轨迹。在实践中，方程 （8） 中的积分是根据 和
W i → f irr = ∑ k = 0 N − 1 Δ λ k ( ∂ H ∂ λ ) Γ ( k Δ t ) W i → f irr = ∑ k = 0 N − 1   Δ λ k ∂ H ∂ λ Γ ( k Δ t ) W_(irarrf)^(irr)=sum_(k=0)^(N-1)Deltalambda_(k)((del H)/(del lambda))_(Gamma(k Delta t))W_{\mathrm{i} \rightarrow \mathrm{f}}^{\mathrm{irr}}=\sum_{k=0}^{N-1} \Delta \lambda_{k}\left(\frac{\partial H}{\partial \lambda}\right)_{\boldsymbol{\Gamma}(k \Delta t)}
where Δ t Δ t Delta t\Delta t is the MD time step, k k kk is the time step number, N N NN is the total number of time steps in which λ λ lambda\lambda varies from λ i λ i lambda_(i)\lambda_{\mathrm{i}} to λ f λ f lambda_(f)\lambda_{\mathrm{f}} and Δ λ k ≡ λ k + 1 − λ k Δ λ k ≡ λ k + 1 − λ k Deltalambda_(k)-=lambda_(k+1)-lambda_(k)\Delta \lambda_{k} \equiv \lambda_{k+1}-\lambda_{k} is the discretization step of the switching parameter λ λ lambda\lambda at time step k k kk.
其中 Δ t Δ t Delta t\Delta t 是 MD 时间步长， k k kk 是时间步长编号， N N NN 是 从 λ i λ i lambda_(i)\lambda_{\mathrm{i}} 到 λ f λ f lambda_(f)\lambda_{\mathrm{f}} λ λ lambda\lambda 变化的时间步长总数， Δ λ k ≡ λ k + 1 − λ k Δ λ k ≡ λ k + 1 − λ k Deltalambda_(k)-=lambda_(k+1)-lambda_(k)\Delta \lambda_{k} \equiv \lambda_{k+1}-\lambda_{k} 是时间步 k k kk 长 中开关参数 λ λ lambda\lambda 的离散化步长。
As mentioned previously, due to the nonequilibrium character of the process the average value of several realizations is subject to a systematic error due to the dissipated heat. But if the process is sufficiently slow for linear response theory to be accurate, it can be eliminated by combining the results of forward and backward switching processes as described by Eqs. (3) and (4).
如前所述，由于该过程的非平衡特性，由于散热，几个实现的平均值会受到系统误差的影响。但是，如果该过程足够慢，线性响应理论是准确的，则可以通过将 Eqs 描述的正向和反向切换过程的结果相结合来消除它。（3） 和 （4）。
Finally, it is important to emphasize that, before the nonequilibrium process is initiated, the system should be equilibrated at either λ i λ i lambda_(i)\lambda_{\mathrm{i}} or λ f λ f lambda_(f)\lambda_{\mathrm{f}}, depending on the sense of the switching process [2].
最后，需要强调的是，在非平衡过程开始之前，系统应该在 λ i λ i lambda_(i)\lambda_{\mathrm{i}} 或 λ f λ f lambda_(f)\lambda_{\mathrm{f}} 处平衡，具体取决于开关过程的意义 [2]。
Next we discuss the application of two specific thermodynamic paths that allow the calculation of the freeenergy difference between two systems described by different Hamiltonians as well as the computation of the temperature variation of the free energy for a given system Hamiltonian.
接下来，我们讨论了两个特定热力学路径的应用，这些路径允许计算由不同哈密顿量描述的两个系统之间的自由能差，以及计算给定系统哈密顿量的自由能温度变化。

B. Free-energy difference between two systems: Frenkel-Ladd path
B. 两个系统之间的自由能差：Frenkel-Ladd 路径

We choose the parametrical Hamiltonian H ( λ ) H ( λ ) H(lambda)H(\lambda) to be of the particular form
我们选择参数化哈密顿量 H ( λ ) H ( λ ) H(lambda)H(\lambda) 为特定形式
H ( λ ) = λ H f + ( 1 − λ ) H i H ( λ ) = λ H f + ( 1 − λ ) H i H(lambda)=lambdaH_(f)+(1-lambda)H_(i)H(\lambda)=\lambda H_{\mathrm{f}}+(1-\lambda) H_{\mathrm{i}}
where H i H i H_(i)H_{\mathrm{i}} and H f H f H_(f)H_{\mathrm{f}} represent two different system Hamiltonians. Setting λ i = 0 λ i = 0 lambda_(i)=0\lambda_{\mathrm{i}}=0 and λ f = 1 λ f = 1 lambda_(f)=1\lambda_{\mathrm{f}}=1, respectively Eq. (7) represents the free-energy difference between these systems, i.e.,
其中 H i H i H_(i)H_{\mathrm{i}} 和 H f H f H_(f)H_{\mathrm{f}} 表示两个不同的系统哈密顿量。设置 λ i = 0 λ i = 0 lambda_(i)=0\lambda_{\mathrm{i}}=0 和 λ f = 1 λ f = 1 lambda_(f)=1\lambda_{\mathrm{f}}=1 分别表示这些系统之间的自由能差，即
Δ F = F f − F i = ∫ 0 1 d λ ⟨ H f − H i ⟩ λ ≡ W i → f rev Δ F = F f − F i = ∫ 0 1   d λ H f − H i λ ≡ W i → f rev {:[Delta F=F_(f)-F_(i)=int_(0)^(1)dlambda(:H_(f)-H_(i):)_(lambda)],[-=W_(irarrf)^(rev)]:}\begin{aligned} \Delta F=F_{\mathrm{f}}-F_{\mathrm{i}} & =\int_{0}^{1} \mathrm{~d} \lambda\left\langle H_{\mathrm{f}}-H_{\mathrm{i}}\right\rangle_{\lambda} \\ & \equiv W_{\mathrm{i} \rightarrow \mathrm{f}}^{\mathrm{rev}} \end{aligned}
and the corresponding forward irreversible work estimator is determined as
，相应的前向不可逆功估算器确定为
W i → f irr = ∫ 0 t s d t d λ d t [ H f ( Γ ( t ) ) − H i ( Γ ( t ) ) ] W i → f irr = ∫ 0 t s   d t d λ d t H f ( Γ ( t ) ) − H i ( Γ ( t ) ) W_(irarrf)^(irr)=int_(0)^(t_(s))dt((d)lambda)/((d)t)[H_(f)(Gamma(t))-H_(i)(Gamma(t))]W_{\mathrm{i} \rightarrow \mathrm{f}}^{\mathrm{irr}}=\int_{0}^{t_{\mathrm{s}}} \mathrm{~d} t \frac{\mathrm{~d} \lambda}{\mathrm{~d} t}\left[H_{\mathrm{f}}(\boldsymbol{\Gamma}(t))-H_{\mathrm{i}}(\boldsymbol{\Gamma}(t))\right]
The Frenkel-Ladd (FL) path [10] uses this concept for computing the absolute free energy of atomic solids. The initial Hamiltonian H i H i H_(i)H_{\mathrm{i}} in Eq. (10) is chosen to be that of the system of interest for which we wish to compute the free energy, usually of the form
Frenkel-Ladd （FL） 路径 [10] 使用这个概念来计算原子固体的绝对自由能。选择方程（10）中的初始哈密顿量 H i H i H_(i)H_{\mathrm{i}} 是我们希望计算自由能的感兴趣系统的量程，通常为以下形式
H i ≡ H 0 = ∑ i = 1 N p i 2 2 m + U ( r ) H i ≡ H 0 = ∑ i = 1 N   p i 2 2 m + U ( r ) H_(i)-=H_(0)=sum_(i=1)^(N)(p_(i)^(2))/(2m)+U(r)H_{\mathrm{i}} \equiv H_{0}=\sum_{i=1}^{N} \frac{\mathbf{p}_{i}^{2}}{2 m}+U(\mathbf{r})
where U ( r ) U ( r ) U(r)U(\mathbf{r}) is some interaction potential and m m mm is the particle mass. We assume that at a certain temperature and volume we know that the stable phase of this system is a certain solid structure (e.g., face-centered cubic). As the second Hamiltonian in the path of Eq. (10) we consider that of a system of noninteracting particles of mass m m mm, each of which is attached to a lattice point by a 3dimensional harmonic spring. The crystallographic lattice to which these particles are connected corresponds precisely to that of the equilibrium phase of interest of system H 0 H 0 H_(0)H_{0}. The Hamiltonian of this harmonic reference system, known as an Einstein crystal, can be written as
其中 U ( r ) U ( r ) U(r)U(\mathbf{r}) 是一些相互作用势，是 m m mm 粒子质量。我们假设在一定的温度和体积下，我们知道这个系统的稳定相是某种固体结构（例如，面心立方）。作为方程（10）路径中的第二个哈密顿量，我们考虑一个质量不相互作用的粒子系统 m m mm ，每个粒子都通过三维谐波弹簧连接到晶格点。这些粒子所连接的晶格与系统 H 0 H 0 H_(0)H_{0} 的目标平衡相 的晶格精确对应。这个谐波参考系统的哈密顿量，称为爱因斯坦晶体，可以写成
H f ≡ H E = ∑ i = 1 N [ p i 2 2 m + 1 2 m ω 2 ( r i − r i 0 ) 2 ] H f ≡ H E = ∑ i = 1 N   p i 2 2 m + 1 2 m ω 2 r i − r i 0 2 H_(f)-=H_(E)=sum_(i=1)^(N)[(p_(i)^(2))/(2m)+(1)/(2)momega^(2)(r_(i)-r_(i)^(0))^(2)]H_{\mathrm{f}} \equiv H_{\mathrm{E}}=\sum_{i=1}^{N}\left[\frac{\mathbf{p}_{i}^{2}}{2 m}+\frac{1}{2} m \omega^{2}\left(\mathbf{r}_{i}-\mathbf{r}_{i}^{0}\right)^{2}\right]
where ω ω omega\omega is the oscillator frequency and r i 0 r i 0 r_(i)^(0)\mathbf{r}_{i}^{0} is the equilibrium position of particle i i ii in system H 0 H 0 H_(0)H_{0}. Its Helmholtz free energy is known analytically, namely,
其中 ω ω omega\omega 是振荡器频率， r i 0 r i 0 r_(i)^(0)\mathbf{r}_{i}^{0} 是粒子 i i ii 在系统中 H 0 H 0 H_(0)H_{0} 的平衡位置。它的亥姆霍兹自由能在解析上是已知的，即
F E ( N , V , T ) = 3 N k B T ln ⁡ ( ℏ ω k B T ) F E ( N , V , T ) = 3 N k B T ln ⁡ ℏ ω k B T F_(E)(N,V,T)=3Nk_(B)T ln((ℏomega)/(k_(B)T))F_{\mathrm{E}}(N, V, T)=3 N k_{\mathrm{B}} T \ln \left(\frac{\hbar \omega}{k_{\mathrm{B}} T}\right)
By estimating the reversible work between these two states, combining the results of forward and backward switching processes (with λ ( t = 0 ) = 0 λ ( t = 0 ) = 0 lambda(t=0)=0\lambda(t=0)=0 and λ ( t = t s ) = 1 λ t = t s = 1 lambda(t=t_(s))=1\lambda\left(t=t_{\mathrm{s}}\right)=1 and the opposite for the forward and backward processes, respectively), as described previously, the free energy of interest can be estimated as
如前所述，通过估计这两种状态之间的可逆功，结合正向和反向切换过程的结果（分别为正向和反向过程的 λ ( t = 0 ) = 0 λ ( t = 0 ) = 0 lambda(t=0)=0\lambda(t=0)=0 和 λ ( t = t s ) = 1 λ t = t s = 1 lambda(t=t_(s))=1\lambda\left(t=t_{\mathrm{s}}\right)=1 相反），感兴趣的自由能可以估计为
F 0 ( N , V , T ) = F E ( N , V , T ) + 1 2 ( W i → f irr ― − W f → i irr ― ) F 0 ( N , V , T ) = F E ( N , V , T ) + 1 2 W i → f irr ¯ − W f → i irr ¯ F_(0)(N,V,T)=F_(E)(N,V,T)+(1)/(2)( bar(W_(irarrf)^(irr))- bar(W_(frarri)^(irr)))F_{0}(N, V, T)=F_{\mathrm{E}}(N, V, T)+\frac{1}{2}\left(\overline{W_{\mathrm{i} \rightarrow \mathrm{f}}^{\mathrm{irr}}}-\overline{W_{\mathrm{f} \rightarrow \mathrm{i}}^{\mathrm{irr}}}\right)

C. Temperature dependence of the free energy: the Reversible Scaling path
C. 自由能的温度依赖性：可逆缩放路径

The Reversible Scaling (RS) path [11, 12] is a particular parametric form of the Hamiltonian H ( λ ) H ( λ ) H(lambda)H(\lambda) for which each value of λ λ lambda\lambda corresponds to a particular temperature of a given system Hamiltonian H 0 H 0 H_(0)H_{0}. In this way, applying the nonequilibrium free-energy approach allows the calculation of the temperature dependence of F 0 ( N , V , T ) F 0 ( N , V , T ) F_(0)(N,V,T)F_{0}(N, V, T) from a single constant temperature simulation. In the
可逆缩放 （RS） 路径 [11， 12] 是哈密顿量 H ( λ ) H ( λ ) H(lambda)H(\lambda) 的一种特殊参数形式，其中 的每个 λ λ lambda\lambda 值对应于给定系统 哈密顿量 H 0 H 0 H_(0)H_{0} 的特定温度。通过这种方式，应用非平衡自由能方法可以从单个恒温仿真中计算出 F 0 ( N , V , T ) F 0 ( N , V , T ) F_(0)(N,V,T)F_{0}(N, V, T) 的温度依赖性。在