"""
calphy: a Python library and command line interface for automated free
energy calculations.
Copyright 2021-2026 (c) Sarath Menon, Yury Lysogorskiy, Ralf Drautz
Interdisciplinary Centre for Advanced Materials Simulation (ICAMS),
Ruhr University Bochum, 44801 Bochum, Germany
calphy is published and distributed under the Academic Software Licence v1.0 (ASL).
calphy is distributed in the hope that it will be useful for non-commercial academic
research, but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the LICENSE file for details.
The ASL permits academic non-commercial use only. Contact
sarath.menon@ruhr-uni-bochum.de to enquire about commercial use rights.
More information about the program can be found in:
Menon, Sarath, Yury Lysogorskiy, Jutta Rogal, and Ralf Drautz.
"Automated Free Energy Calculation from Atomistic Simulations." Physical Review Materials 5(10), 2021
DOI: 10.1103/PhysRevMaterials.5.103801
For more information contact:
sarath.menon@ruhr-uni-bochum.de
"""
import numpy as np
import scipy.constants as const
import sys
import math
import os
import warnings
from calphy.splines import (
splines,
sum_spline1,
sum_spline25,
sum_spline50,
sum_spline75,
sum_spline100,
)
from scipy.integrate import cumulative_trapezoid as cumtrapz
from numpy import trapezoid as trapz
from tqdm import tqdm
import pyscal3.core as pc
from ase.io import read
# Constants
h = const.physical_constants["Planck constant in eV/Hz"][0]
hJ = const.physical_constants["Planck constant"][0]
hbar = h / (2 * np.pi)
kb = const.physical_constants["Boltzmann constant in eV/K"][0]
kbJ = const.physical_constants["Boltzmann constant"][0]
Na = const.physical_constants["Avogadro constant"][0]
eV2J = const.eV
J2eV = 6.242e18
# Pressure unit conversion: 1 eV/A^3 = e * 1e25 bar (exactly, given the SI
# definition of the electron volt). Used to convert LAMMPS pressures (bar)
# to the eV/A^3 units that the rest of the integrators work in.
# pressure_in_eV_per_A3 = pressure_in_bar / EV_A3_TO_BAR
EV_A3_TO_BAR = const.e * 1e25
# --------------------------------------------------------------------
# TI PATH INTEGRATION ROUTINES
# --------------------------------------------------------------------
[docs]def integrate_path(
calc,
fwdfilename,
bkdfilename,
solid=True,
):
"""
Get a filename with columns du and dlambda and integrate
Parameters
----------
fwdfilename: string
name of fwd integration file
bkdfilename: string
name of bkd integration file
usecols : list
column numbers to be used from input file
Returns
-------
w : float
irreversible work in switching the system
q : float
heat dissipation during switching of system
"""
natoms = np.array([calc._element_dict[x]["count"] for x in calc.element])
concentration = np.array(
[calc._element_dict[x]["composition"] for x in calc.element]
)
fdata = np.loadtxt(fwdfilename, unpack=True, comments="#")
bdata = np.loadtxt(bkdfilename, unpack=True, comments="#")
if solid:
fdui = fdata[0]
bdui = bdata[0]
fdur = np.zeros(len(fdui))
bdur = np.zeros(len(bdui))
for i in range(calc.n_elements):
if natoms[i] > 0:
fdur += concentration[i] * fdata[i + 1] / natoms[i]
bdur += concentration[i] * bdata[i + 1] / natoms[i]
flambda = fdata[calc.n_elements + 1]
blambda = bdata[calc.n_elements + 1]
else:
fdui = fdata[0]
bdui = bdata[0]
fdur = fdata[1]
bdur = bdata[1]
flambda = fdata[2]
blambda = bdata[2]
fdu = fdui - fdur
bdu = bdui - bdur
fw = trapz(fdu, flambda)
bw = trapz(bdu, blambda)
w = 0.5 * (fw - bw)
q = 0.5 * (fw + bw)
return w, q, flambda
[docs]def find_w(mainfolder, calc, full=False, solid=True, prefix=""):
"""
Integrate the irreversible work and dissipation for independent simulations
Parameters
----------
mainfolder: string
main simulation folder
nsims : int, optional
number of independent simulations, default 5
full : bool, optional
If True return error values, default False
usecols : tuple, optional
Columns to read in from data file. Default (0, 1)
prefix : str, optional
infix inserted into the switching-data filenames, e.g. "leg1" reads
forward_leg1_%d.dat / backward_leg1_%d.dat. Empty string (default) reads
the original forward_%d.dat / backward_%d.dat.
Returns
-------
ws : float
average irreversible work
qs : float
average energy dissipation, only returned if full is True
err : float
Error in free energy, only returned if full is True
"""
ws = []
qs = []
tag = ("_%s" % prefix) if prefix else ""
for i in range(calc.n_iterations):
fwdfilestring = "forward%s_%d.dat" % (tag, i + 1)
fwdfilename = os.path.join(mainfolder, fwdfilestring)
bkdfilestring = "backward%s_%d.dat" % (tag, i + 1)
bkdfilename = os.path.join(mainfolder, bkdfilestring)
w, q, flambda = integrate_path(
calc,
fwdfilename,
bkdfilename,
solid=solid,
)
ws.append(w)
qs.append(q)
wsmean = np.mean(ws)
qsmean = np.mean(qs)
# Uncertainty on the *mean* of n_iterations independent switching runs is the
# standard error of the mean (SEM = sample std / sqrt(N)), not the spread of
# the individual runs. A single run yields no error estimate.
wsstd = np.std(ws, ddof=1) / np.sqrt(len(ws)) if len(ws) > 1 else 0.0
if full:
return wsmean, qsmean, wsstd
else:
return wsmean
[docs]def integrate_rs(
simfolder, f0, t, natoms, p=0, nsims=5, scale_energy=False, return_values=False
):
"""
Carry out the reversible scaling integration
Parameters
----------
simfolder : string
main simulation folder
f0 : float
initial free energy for integration
t : float
initial temperature
nsims : int, optional
number of independent switching
scale_energy: bool, optional
if True, scale energy with switching parameter
Returns
-------
None
Notes
-----
Writes the output in a file reversible_scaling.dat
"""
ws = []
es = []
p = p / EV_A3_TO_BAR
for i in range(1, nsims + 1):
fdx, fp, fvol, flambda = np.loadtxt(
os.path.join(simfolder, "ts.forward_%d.dat" % i), unpack=True, comments="#"
)
bdx, bp, bvol, blambda = np.loadtxt(
os.path.join(simfolder, "ts.backward_%d.dat" % i), unpack=True, comments="#"
)
# Both sweeps start at lambda=1 (T0 end). Row-count mismatches can
# arise from integer rounding in the block planner or from a truncated
# forward sweep during phase-transition recovery. Trim from the
# lambda=lf end (the *second* temperature, T_stop) of the longer
# array so the two grids remain anchored at lambda=1 (T0).
# For a heating sweep (T0<Tf) this is the high-T tail; for a cooling
# sweep (T0>Tf) this is the low-T tail. Either way:
# forward rows: lambda[0]=1 … lambda[-1]=lf → trim [:n]
# backward rows: lambda[0]=lf … lambda[-1]=1 → trim [-n:]
nf = len(flambda)
nb = len(blambda)
if nf != nb:
n = min(nf, nb)
if nf > nb:
# forward is longer — drop its tail (lf / T_stop end)
fdx = fdx[:n]
fvol = fvol[:n] # slice before /natoms below
flambda = flambda[:n]
else:
# backward is longer — drop its first rows (lf / T_stop end,
# which is at the START of the backward file)
bdx = bdx[-n:]
bvol = bvol[-n:]
blambda = blambda[-n:]
if scale_energy:
fdx /= flambda
bdx /= blambda
# add pressure contribution
fvol = fvol / natoms
bvol = bvol / natoms
fdx = fdx + p * fvol
bdx = bdx + p * bvol
wf = cumtrapz(fdx, flambda, initial=0)
wb = cumtrapz(bdx[::-1], blambda[::-1], initial=0)
w = (wf + wb) / (2 * flambda)
e = np.max(np.abs((wf - wb) / (2 * flambda)))
ws.append(w)
es.append(e)
e_diss = np.min(es)
wmean = np.mean(ws, axis=0)
# standard error of the mean across the nsims independent sweeps (see find_w)
werr = np.std(ws, axis=0, ddof=1) / np.sqrt(len(ws)) if len(ws) > 1 else np.zeros_like(wmean)
temp = t / flambda
f = f0 / flambda + 1.5 * kb * temp * np.log(flambda) + wmean
outfile = os.path.join(simfolder, "temperature_sweep.dat")
np.savetxt(
outfile,
np.column_stack((temp, f, werr)),
header="reversible-scaling free energy vs temperature\n"
"temperature[K] free_energy[eV/atom] error[eV/atom]",
)
if not return_values:
return None, e_diss
else:
return (temp, f, werr), e_diss
[docs]def integrate_ps(simfolder, f0, natoms, pi, pf, nsims=1, return_values=False):
"""
Carry out the reversible scaling integration
Parameters
----------
simfolder : string
main simulation folder
f0 : float
initial free energy for integration
nsims : int, optional
number of independent switching
Returns
-------
None
Notes
-----
Writes the output in a file pressure_sweep.dat
"""
ws = []
for i in range(1, nsims + 1):
_, fp, fvol, _ = np.loadtxt(
os.path.join(simfolder, "ps.forward_%d.dat" % i), unpack=True, comments="#"
)
_, bp, bvol, _ = np.loadtxt(
os.path.join(simfolder, "ps.backward_%d.dat" % i), unpack=True, comments="#"
)
fvol = fvol / natoms
bvol = bvol / natoms
fp = fp / EV_A3_TO_BAR
bp = bp / EV_A3_TO_BAR
wf = cumtrapz(fvol, fp, initial=0)
wb = cumtrapz(bvol[::-1], bp[::-1], initial=0)
w = (wf + wb) / 2
ws.append(w)
wmean = np.mean(ws, axis=0)
# standard error of the mean across the nsims independent sweeps (see find_w)
werr = np.std(ws, axis=0, ddof=1) / np.sqrt(len(ws)) if len(ws) > 1 else np.zeros_like(wmean)
press = np.linspace(pi, pf, len(wmean))
f = f0 + wmean
if not return_values:
outfile = os.path.join(simfolder, "pressure_sweep.dat")
np.savetxt(
outfile,
np.column_stack((press, f, werr)),
header="pressure-scaling free energy vs pressure\n"
"pressure[bar] free_energy[eV/atom] error[eV/atom]",
)
else:
return (press, f, werr)
[docs]def integrate_mass(ref_mass, target_masses, target_counts, temperature, natoms):
mcorsum = 0
for i in range(len(target_masses)):
mcorsum += (
1.5
* kb
* temperature
* (1.0 * (target_counts[i] / natoms) * np.log(target_masses[i] / ref_mass))
)
return mcorsum
# --------------------------------------------------------------------
# REF. STATE ROUTINES: SOLID
# --------------------------------------------------------------------
[docs]def get_einstein_crystal_fe(
calc, vol, k, cm_correction=True, return_contributions=False, quantum=False
):
"""
Get the free energy of einstein crystal
Parameters
----------
calc : Calculation object
contains all input parameters
vol : float
converged volume per atom
k : spring constant, float
units - eV/Angstrom^2
cm_correction : bool, optional, default - True
add the centre of mass correction to free energy
return_contributions: bool, optional, default - True
If True, return individual contributions to the reference free energy.
quantum : bool, optional, default - False
If True, evaluate the *quantum* harmonic-oscillator free energy of
the Einstein crystal. This is the reference required when the
switching MD is driven by the Dammak quantum thermal bath, which
samples a quantum-statistical distribution rather than the
classical Boltzmann one. In the high-temperature limit
(k_B T >> hbar omega) this reduces to the classical expression.
Returns
-------
F_tot : float
total free energy of reference crystal
F_e : float
Free energy of Einstein crystal without centre of mass correction. Only if `return_contributions` is True.
F_cm : float
centre of mass correction. Only if `return_contributions` is True.
Notes
-----
The classical equations for free energy of Einstein crystal and centre
of mass correction are from https://doi.org/10.1063/5.0044833.
Quantum branch: each Cartesian mode i has angular frequency
omega_i = sqrt(k_i / m_i), and the quantum harmonic-oscillator free
energy per mode is
f_i = (1/2) hbar omega_i + k_B T ln(1 - exp(-hbar omega_i / k_B T)).
The reported F_e is sum_i f_i / N_atoms (3 modes per atom). The CM
correction is retained in its classical form: it is a phase-space
integration correction that does not change form under the QTB
sampling at leading order.
"""
# temperature
temp = calc._temperature
# natoms
natoms = np.sum([calc._element_dict[x]["count"] for x in calc.element])
# convert a to m3
vol = vol * 1e-30
# whats the beta
beta = 1 / (kbJ * temp)
# create an array of mass
mass = []
for x in calc.element:
for count in range(calc._element_dict[x]["count"]):
mass.append(calc._element_dict[x]["mass"])
mass = np.array(mass)
# convert mass to kg
mass = (mass / Na) * 1e-3
# create an array of k as well
karr = []
for c, x in enumerate(calc.element):
for count in range(calc._element_dict[x]["count"]):
karr.append(k[c])
k = np.array(karr)
# convert k from ev/A2 to J/m2
k = k * (eV2J / 1e-20)
if quantum:
# quantum harmonic oscillator reference (Dammak/QTB sampling)
hbar = hJ / (2 * np.pi)
omega = np.sqrt(k / mass) # rad/s, per atom
x = hbar * omega / (kbJ * temp)
# per-mode free energy in Joules; 3 modes per atom
f_mode = 0.5 * hbar * omega + kbJ * temp * np.log1p(-np.exp(-x))
F_e_J = 3.0 * np.sum(f_mode) / natoms # J per atom
F_e = F_e_J / eV2J # convert J -> eV
else:
# classical Einstein crystal
Z_e = ((beta**2 * k * hJ**2) / (4 * np.pi**2 * mass)) ** 1.5
F_e = np.log(Z_e)
F_e = kb * temp * np.sum(F_e) / natoms # *J2eV #convert back to eV
# now get the cm correction
if cm_correction:
mass_sum = np.sum(mass)
mu = mass / mass_sum
mu2_over_k = mu**2 / k
mu2_over_k_sum = np.sum(mu2_over_k)
prefactor = vol
F_cm = np.log(prefactor * (beta / (2 * np.pi * mu2_over_k_sum)) ** 1.5)
F_cm = kb * temp * F_cm / natoms # convert to eV
else:
F_cm = 0
F_tot = F_e - F_cm
if return_contributions:
return F_e, -F_cm
return F_tot
# --------------------------------------------------------------------
# REF. STATE ROUTINES: LIQUID
# --------------------------------------------------------------------
[docs]def get_ideal_gas_fe(temp, rho, natoms, mass, concentration):
"""
Get the free energy of an single/binary ideal gas
Parameters
----------
temp : temperature, float
the reference temperature in K
rho : number density, float
units - no of atoms/ angstrom^3
natoms: int
total number of atoms
mass : atomic mass, float
units - g/mol
xa : concentration of species a, float, optional
default 1
xb : concentration of species b, float, optional
default 0
Returns
-------
fe : float
free energy/atom of ideal gas system
"""
# find mass of one particle
mass = np.array(mass) / Na
beta = 1 / (kb * temp) # units - eV
# omega needs to be in m
omega = (beta * h * h / (2 * np.pi * mass)) ** 0.5
# convert omega
omega = omega * (const.eV / 1e-3) ** 0.5
# the above is in metres - change to Angstrom
omega = omega * 1e10
prefactor = 1 / beta
fe = 0
for count, conc in enumerate(concentration):
if concentration[count] > 0:
fe += conc * (3 * np.log(omega[count]) + np.log(rho) - 1 + np.log(conc))
# return prefactor*(ta + tb + (1/(2*natoms))*np.log(2*np.pi*natoms))
return prefactor * fe
[docs]def get_uhlenbeck_ford_fe(temp, rho, p, sigma):
"""
Get the excess free energy of Uhlenbeck-Ford model
Parameters
----------
temp : temperature, float
units - K
rho : density, float
units - no of atoms/ angstrom^3
p : uf scale, float
sigma : uf length scale, float
Returns
-------
fe : float
excess free energy/atom of uf system
"""
x = (0.5 * (np.pi * sigma * sigma) ** 1.5) * rho
_, fe = find_fe(p, x)
beta = 1 / (kb * temp)
fe = fe / beta
return fe
[docs]def press(x, coef):
"""
Find pressure of system
Parameters
----------
x : float
x value for UF system
coef : list of floats
coefficients
Returns
-------
result : float, optional
result pressure
"""
result = coef[0] * (x**3) + coef[1] * (x**2) + coef[2] * x + coef[3]
return result
[docs]def fe(x, coef, sum_spline, index):
"""
Fe inbuilt method
"""
if x < 0.0025:
result = coef[0] * (x**2) / 2.0 + coef[1] * x
return result
elif x < 0.1:
if x * 10000 % 25 == 0:
return sum_spline[index - 1]
else:
x_0 = 0.0025 * int(x * 400)
elif x < 1:
if x * 1000 % 25 == 0:
return sum_spline[index - 1]
else:
x_0 = 0.025 * int(x * 40)
elif x < 4:
if x * 100 % 10 == 0:
return sum_spline[index - 1]
else:
x_0 = 0.1 * int(x * 10)
else:
return sum_spline[index]
result = (
sum_spline[index - 1]
+ coef[0] * (x**2.0 - x_0**2.0) / 2.0
+ coef[1] * (x - x_0)
+ (coef[2] - 1.0) * math.log(x / x_0)
- coef[3] * (1.0 / x - 1.0 / x_0)
)
return result
[docs]def find_fe(p, x):
"""
Find free energy of UF system
Parameters
----------
x : float
x value of system
coef : list of floats
Coefficients of system
Returns
-------
fe : float
free energy of UF system
"""
if not p in splines:
raise ValueError("Invalid p. Valid numbers are: 1, 25, 50, 75, and 100.")
if (x <= 0.0) or (x > 4.0):
raise ValueError("Invalid x. Valid numbers are 0.0 < x <= 4.0")
table1 = splines[p]
if x < 0.1:
index = 0 + int(x * 400)
elif x < 1:
index = 40 + int((x * 40 - 4))
elif x < 4:
index = 76 + int((x * 10 - 10))
else:
index = 105
coef = table1[index]
pressure = press(x, coef)
if p == 1:
sum_spline = sum_spline1
elif p == 25:
sum_spline = sum_spline25
elif p == 50:
sum_spline = sum_spline50
elif p == 75:
sum_spline = sum_spline75
else:
sum_spline = sum_spline100
free_energy = fe(x, coef, sum_spline, index)
return pressure, free_energy
# --------------------------------------------------------------------
# PHASE DIAGRAM ROUTINES
# --------------------------------------------------------------------
[docs]def calculate_entropy_mix(conc):
"""
Calculate the entropy of mixing
Parameters
----------
conc : float
concentration
Returns
-------
s: float
entropy
"""
s = -kb * (conc * np.log(conc) + (1 - conc) * np.log(1 - conc))
return s
[docs]def calculate_fe_impurity(temp, natoms, fepure, feimpure):
"""
Calculate energy change of mixing, imput energies are
in eV/atom
Parameters
----------
temp : float
temperature
natoms : int
number of atoms
fepure : float
free energy of pure phase
feimpure : float
free energy of impure phase
Returns
-------
dg : float
entropy of mixing
"""
dg = feimpure * natoms - fepure * natoms + kb * temp * np.log(natoms)
return dg
[docs]def calculate_fe_mix(temp, fepure, feimpure, concs, natoms=4000):
"""
Calculate energy of mixing
Parameters
----------
temp : float
temperature
fepure : float
free energy of the pure phase
feimpure : float
energy due to impurity
concs : list of floats
concentration array
natoms : int
number of atoms
Returns
-------
fes : list of floats
free energy with concentration
"""
if concs[0] == 0:
print("zero is autodone")
concs = concs[1:]
fes = [fepure]
s = calculate_entropy_mix(concs)
dg = calculate_fe_impurity(temp, natoms, fepure, feimpure)
fe_conc = fepure + concs * dg - temp * s
for f in fe_conc:
fes.append(f)
return fes