"""
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 os
import numpy as np
import yaml
import matplotlib.pyplot as plt
import warnings
import pandas as pd
[docs]def read_report(folder):
"""
Read the finished calculation report
Parameters
----------
folder: string
folder from which calculation is to be read
Returns
-------
data: dict
dictionary with results
"""
repfile = os.path.join(folder, "report.yaml")
if not os.path.exists(repfile):
raise FileNotFoundError(f"file {repfile} not found")
with open(repfile, "r") as fin:
data = yaml.safe_load(fin)
return data
def _extract_error(errfile):
error_code = None
try:
if os.path.exists(errfile):
with open(errfile, "r") as fin:
for line in fin:
if "calphy.errors" in line:
break
error_code = line.split(":")[0].split(".")[-1]
except (OSError, IndexError):
pass
return error_code
def _read_sweep_data(folder, stride=1):
"""
Read the per-lambda switching data of every reversible-scaling replica.
Parameters
----------
folder : str
Calculation folder holding ``ts.forward_<i>.dat`` /
``ts.backward_<i>.dat``.
stride : int, optional
Keep every n-th sample.
Returns
-------
dict
``forward_energy_diff`` / ``backward_energy_diff`` (energy differential
per lambda) and ``forward_lambda`` / ``backward_lambda``, each a list
with one array per replica, or None where no replica was found.
Raises
------
Exception
Propagates any read/parse failure; the caller decides whether one bad
calculation should be fatal.
"""
f_ediffs, b_ediffs, f_lambdas, b_lambdas = [], [], [], []
i = 1
while True:
fwdfile = os.path.join(folder, f"ts.forward_{i}.dat")
bkdfile = os.path.join(folder, f"ts.backward_{i}.dat")
if not (os.path.exists(fwdfile) and os.path.exists(bkdfile)):
break
# ndmin=2 so a single-row sweep still unpacks as columns
fwd = np.loadtxt(fwdfile, comments="#", ndmin=2)[::stride]
bkd = np.loadtxt(bkdfile, comments="#", ndmin=2)[::stride]
fdx, flambda = fwd[:, 0], fwd[:, 3]
bdx, blambda = bkd[:, 0], bkd[:, 3]
f_ediffs.append(fdx / flambda)
b_ediffs.append(bdx / blambda)
f_lambdas.append(flambda)
b_lambdas.append(blambda)
i += 1
return {
"forward_energy_diff": f_ediffs or None,
"backward_energy_diff": b_ediffs or None,
"forward_lambda": f_lambdas or None,
"backward_lambda": b_lambdas or None,
}
[docs]def gather_results(
mainfolder,
reduce_composition=True,
extract_phase_prefix=False,
include_sweep_data=False,
sweep_data_stride=1,
):
"""
Gather results from all subfolders in a given folder into a Pandas DataFrame
Parameters
----------
mainfolder: string
folder where calculations are stored
reduce_composition: bool
If True, per species composition arrays are added.
Might be redundant.
extract_phase_prefix: bool
Should be used in conjuction with phase diagram mode.
Extracts the prefix and add it as a phase_name column.
include_sweep_data: bool, optional, default False
Load the raw per-lambda switching data of every reversible-scaling
replica into the frame. Off by default because it is expensive:
``ts.forward_*.dat`` / ``ts.backward_*.dat`` are written by
``fix print 1``, i.e. one row per MD step, so a sweep of
``n_switching_steps`` keeps roughly ``n_switching_steps * 32`` bytes
per replica per calculation resident in the frame (~1.6 MB per
calculation for a 50000-step sweep, scaling with ``n_iterations``).
Enable it when diagnosing *where* along a sweep forward and backward
diverge; leave it off to gather free energies.
sweep_data_stride: int, optional, default 1
Keep only every n-th sample of the sweep data. Only meaningful with
``include_sweep_data=True``; a stride of 10--100 preserves the shape
of the hysteresis while cutting the memory cost proportionally.
Returns
-------
df: pandas DataFrame
DataFrame with results. In addition to the columns produced
previously, this also includes:
- free_energy_error: array (ts/tscale, from temperature_sweep.dat)
or 0.0 (fe/alchemy/composition_scaling); the statistical
standard error of the mean free energy, not a hysteresis check.
- dissipation: mean switching dissipation (fe/alchemy) or NaN (ts/tscale)
- ts_dissipation: max hysteresis over a ts/tscale sweep, or NaN otherwise
- ts_dissipation_high: True where that hysteresis exceeded
tolerance.dissipation, i.e. the sweep did not stay reversible and the
free energy it produced should not be trusted; False when it was
within tolerance, NaN where no verdict was recorded
- forward_energy_diff / backward_energy_diff: list of arrays, one per
reversible-scaling replica (ts/tscale only, else None); the raw
per-lambda energy differential, useful to see *where* along the
sweep forward/backward diverge (phase-transition diagnostic).
Populated only when ``include_sweep_data`` is True, else None.
- forward_lambda / backward_lambda: matching lambda arrays for the
above (ts/tscale only, else None)
Notes
-----
A calculation whose sweep files are unreadable (a job killed mid-sweep
leaves a ragged final row) does not abort the gather: its sweep columns
stay None and the reason is recorded in ``error_code``.
"""
if sweep_data_stride < 1:
raise ValueError("sweep_data_stride must be >= 1")
try:
import pandas as pd
except ImportError:
raise ImportError("Please install pandas to use this function")
unique_elements = []
datadict = {}
datadict["calculation_mode"] = []
datadict["status"] = []
datadict["temperature"] = []
datadict["pressure"] = []
datadict["free_energy"] = []
datadict["free_energy_error"] = []
datadict["dissipation"] = []
datadict["ts_dissipation"] = []
datadict["ts_dissipation_high"] = []
datadict["reference_phase"] = []
datadict["error_code"] = []
datadict["composition"] = []
datadict["calculation"] = []
datadict["ideal_entropy"] = []
datadict["phase_name"] = []
datadict["reference_composition"] = []
datadict["forward_energy_diff"] = []
datadict["backward_energy_diff"] = []
datadict["forward_lambda"] = []
datadict["backward_lambda"] = []
folders = next(os.walk(mainfolder))[1]
for folder in folders:
# adjust for pyiron folder, see
if folder.split("_")[-1] == "hdf5":
# this could be a pyiron calc
withouthdf = folder.split("_hdf5")[0]
folder = f"{folder}/{withouthdf}"
inpfile = os.path.join(mainfolder, folder, "input_file.yaml")
# print(inpfile)
if not os.path.exists(inpfile):
continue
# ok, valid calculation, try to parse input file to get info
with open(inpfile, "r") as fin:
inp = yaml.safe_load(fin)
# grab the first calculation
inp = inp["calculations"][0]
# mode
mode = inp["mode"]
datadict["calculation_mode"].append(mode)
datadict["temperature"].append(inp["temperature"])
datadict["pressure"].append(inp["pressure"])
datadict["reference_phase"].append(inp["reference_phase"])
datadict["phase_name"].append(inp["phase_name"])
datadict["reference_composition"].append(inp["reference_composition"])
datadict["composition"].append(None)
datadict["ideal_entropy"].append(0)
datadict["calculation"].append(folder)
# check output file
outfile = os.path.join(mainfolder, folder, "report.yaml")
datadict["error_code"].append(None)
datadict["free_energy_error"].append(np.nan)
datadict["dissipation"].append(np.nan)
datadict["ts_dissipation"].append(np.nan)
datadict["ts_dissipation_high"].append(np.nan)
datadict["forward_energy_diff"].append(None)
datadict["backward_energy_diff"].append(None)
datadict["forward_lambda"].append(None)
datadict["backward_lambda"].append(None)
# print(inpfile)
if not os.path.exists(outfile):
datadict["status"].append("False")
datadict["free_energy"].append(np.nan)
# check if error file is found
errfile = os.path.join(os.getcwd(), mainfolder, folder + ".sub.err")
datadict["error_code"][-1] = _extract_error(errfile)
continue
if mode in ["fe", "alchemy", "composition_scaling"]:
datadict["status"].append("True")
datadict["free_energy_error"][-1] = 0.0
# ok, valid calculation, try to parse input file to get info
with open(outfile, "r") as fin:
out = yaml.safe_load(fin)
datadict["free_energy"].append(out["results"]["free_energy"])
# scalar quality metrics calphy already computes: qdiss (mean switching
# dissipation, fe/alchemy) and max hysteresis over a ts/tscale sweep
datadict["dissipation"][-1] = out["results"].get("dissipation", np.nan)
datadict["ts_dissipation"][-1] = out["results"].get("ts_dissipation", np.nan)
datadict["ts_dissipation_high"][-1] = out["results"].get(
"ts_dissipation_high", np.nan
)
# add normal composition
el_arr = np.array(out["input"]["element"].split(" ")).astype(str)
comp_arr = np.array(out["input"]["concentration"].split(" ")).astype(float)
composition = {x: y for x, y in zip(el_arr, comp_arr)}
datadict["composition"][-1] = composition
if mode == "composition_scaling":
# we need to update composition
compdict = inp["composition_scaling"]["output_chemical_composition"]
maxatoms = np.sum([val for key, val in compdict.items()])
for key, val in compdict.items():
compdict[key] = val / maxatoms
datadict["composition"][-1] = compdict
el_arr = list(compdict.keys())
# we also need to update entropy
if "entropy_contribution" in out["results"].keys():
datadict["ideal_entropy"][-1] = (
-1 * out["results"]["entropy_contribution"]
)
for el in el_arr:
if el not in unique_elements:
unique_elements.append(el)
# parse extra info
if mode in ["ts", "tscale"]:
datafile = os.path.join(
os.getcwd(), mainfolder, folder, "temperature_sweep.dat"
)
if os.path.exists(datafile):
datadict["status"].append("True")
t, f, ferr = np.loadtxt(datafile, unpack=True, usecols=(0, 1, 2))
datadict["temperature"][-1] = t
datadict["free_energy"][-1] = f
datadict["free_energy_error"][-1] = ferr
# Raw per-lambda energy differential of each forward/backward
# switching replica; diverging forward vs. backward curves
# signal hysteresis/a structural change during the sweep
# (ts_dissipation above is just the max of this difference).
#
# Opt-in: these files carry one row per MD step, so loading
# them for every calculation dominates both the runtime and
# the size of the returned frame.
if include_sweep_data:
try:
sweeps = _read_sweep_data(
os.path.join(mainfolder, folder), sweep_data_stride
)
except Exception as exc:
# a job killed mid-sweep leaves a ragged final row;
# that must not cost us every other calculation
warnings.warn(
"could not read switching data in %s: %s: %s"
% (folder, type(exc).__name__, exc),
RuntimeWarning,
)
datadict["error_code"][-1] = "unreadable switching data: %s" % (
type(exc).__name__,
)
else:
for key, val in sweeps.items():
datadict[key][-1] = val
else:
datadict["status"].append("False")
errfile = os.path.join(os.getcwd(), mainfolder, folder + ".sub.err")
datadict["error_code"][-1] = _extract_error(errfile)
if reduce_composition:
unique_element_dict = {x: [] for x in unique_elements}
for x in datadict["composition"]:
if x is None:
for key in unique_element_dict.keys():
unique_element_dict[key].append(0)
else:
for el in unique_elements:
if el in x.keys():
unique_element_dict[el].append(x[el])
else:
unique_element_dict[el].append(0)
# add the keys to datadict
for key, val in unique_element_dict.items():
datadict[key] = val
if not extract_phase_prefix:
del datadict["phase_name"]
df = pd.DataFrame(data=datadict)
return df
def _entropy(compC, comp_init=[1, 0]):
"""
l1: initial concentration [Au, Cu]
l2: final concentration [Au, Cu]
"""
compA = 1 - compC
def _log(val):
if val == 0:
return 0
else:
return np.log(val)
dA = compA * _log(compA) - comp_init[0] * _log(comp_init[0])
dC = compC * _log(compC) - comp_init[1] * _log(comp_init[1])
return kb * (dA + dC)
[docs]def clean_df(df, reference_element, combine_direct_calculations=False, smooth=False):
"""
Clean a parsed dataframe and drop unnecessary columns. This gets it ready for further processing
Note that `gather_results` should be run with `reduce_composition` and `extract_phase_name` for this to work.
Parameters
----------
df: DataFrame
dataframe parsed by `gather_results` with `reduce_composition=True`.
reference_element: str
reference element from the compositions, which will be renamed to `composition`
combine_direct_calculations: bool, optional
If True, combine direct calculations by fitting to produce temperature and free energy arrays
If used, an extra column `error` with RMSE of the fitting is also created
smooth : bool, optional
If True, smooth the F(T) data using the thermodynamic basis
``[1, T, T ln T, T²]``. If False (default), return the raw data
points without smoothing.
Returns
-------
df: DataFrame
combined, finished DataFrame
"""
if "phase_name" not in df.keys():
raise ValueError("phase_name key is not found, maybe add it?")
df = df.loc[df.status == "True"]
df = df.drop(
labels=[
"status",
"pressure",
"reference_phase",
"error_code",
"composition",
"calculation",
],
axis="columns",
)
phases = df.groupby(df.phase_name)
phases = [phases.get_group(x) for x in phases.groups]
df_dict = {}
for phase in phases:
if combine_direct_calculations:
gb = phase.groupby(by=reference_element)
gbs = [gb.get_group(x) for x in gb.groups]
fes = []
tes = []
errors = []
comps = []
mode_list = []
entropies = []
is_refs = []
for exdf in gbs:
temps = np.array(exdf.temperature.values)
fe = np.array(exdf.free_energy.values)
modes = np.array(exdf.calculation_mode.values)
entropy = np.array(exdf.ideal_entropy.values)
comp_ref = np.array(exdf.reference_composition.values)
unique_modes = np.unique(modes)
if len(unique_modes) > 1:
warnings.warn("mixing calculations from more than one mode!")
unique_mode = unique_modes[0]
# REMEMBER TO SORT EVERYTHING
args = np.argsort(temps)
temps = temps[args]
fe = fe[args]
# entropy = entropy[args]
# print(fe, entropy)
# print(len(fe), len(entropy))
if smooth:
# Thermodynamic basis: F(T) = a + b*T + c*T*ln(T) + d*T²
T_ = temps.astype(float)
basis = np.column_stack(
[np.ones_like(T_), T_, T_ * np.log(T_), T_**2]
)
coeffs_t, _, _, _ = np.linalg.lstsq(basis, fe, rcond=None)
fe_eval = basis @ coeffs_t
error = float(np.sqrt(np.mean((fe_eval - fe) ** 2)))
temp_arr = np.arange(temps.min(), temps.max() + 1, 1).astype(float)
basis_arr = np.column_stack(
[
np.ones_like(temp_arr),
temp_arr,
temp_arr * np.log(temp_arr),
temp_arr**2,
]
)
fe_arr = basis_arr @ coeffs_t
else:
fe_arr = fe
temp_arr = temps
error = 0
fes.append(fe_arr)
tes.append(temp_arr)
errors.append(error)
entropies.append(entropy[0])
comps.append(float(exdf[reference_element].values[0]))
# The fe job IS the reference by definition; composition_scaling jobs are not.
is_refs.append(unique_mode == "fe")
mode_list.append(unique_mode)
# replace df
df = pd.DataFrame(
data={
"temperature": tes,
"free_energy": fes,
"error": errors,
reference_element: comps,
"ideal_entropy": entropies,
"calculation_mode": mode_list,
"is_reference": is_refs,
}
)
df = df.rename(columns={reference_element: "composition"})
df_dict[phase.phase_name.values[0]] = df
return df_dict
[docs]def fix_composition_scaling(dfdict, correct_entropy=True, add_ideal_entropy=False):
"""
Correct composition-scaling free energies by adding the reference
free energy and (optionally) subtracting the ideal-entropy term.
Parameters
----------
dfdict : dict of DataFrame
Output from ``clean_df``.
correct_entropy : bool
If True **and** *add_ideal_entropy* is True, subtract
``T * S_ideal`` from the free energy.
add_ideal_entropy : bool
Controls whether the ideal-entropy correction is applied.
"""
for key, val in dfdict.items():
x = val
if "_entropy_corrected" not in x.columns:
x["_entropy_corrected"] = False
ref_row = x.loc[x.is_reference == True]
ref_fe = np.asarray(ref_row.free_energy.values[0], dtype=float)
ref_temp = np.asarray(ref_row.temperature.values[0], dtype=float)
# Fit reference F(T) with thermodynamic basis
basis_ref = np.column_stack(
[np.ones_like(ref_temp), ref_temp, ref_temp * np.log(ref_temp), ref_temp**2]
)
ref_coeffs, _, _, _ = np.linalg.lstsq(basis_ref, ref_fe, rcond=None)
for index, row in x.iterrows():
if (not row.is_reference) and (
row.calculation_mode == "composition_scaling"
):
T_row = np.asarray(row.temperature, dtype=float)
basis_row = np.column_stack(
[np.ones_like(T_row), T_row, T_row * np.log(T_row), T_row**2]
)
ref_at_T = basis_row @ ref_coeffs
if correct_entropy and add_ideal_entropy:
x.at[index, "free_energy"] = (
row.free_energy - row.temperature * row.ideal_entropy + ref_at_T
)
x.at[index, "_entropy_corrected"] = True
else:
x.at[index, "free_energy"] = row.free_energy + ref_at_T
dfdict[key] = x
return dfdict
[docs]def find_transition_temperature(folder1, folder2, fit_order=4, plot=True):
"""
Find transition temperature where free energy of two phases are equal.
Parameters
----------
folder1: string
directory with temperature scale calculation
folder2: string
directory with temperature scale calculation
fit_order: int, optional
default 4. Order for polynomial fit of temperature vs free energy
plot: bool, optional
default True. Plot the results.
"""
file1 = os.path.join(folder1, "temperature_sweep.dat")
file2 = os.path.join(folder2, "temperature_sweep.dat")
if not os.path.exists(file1):
raise FileNotFoundError(f"{file1} does not exist")
if not os.path.exists(file2):
raise FileNotFoundError(f"{file2} does not exist")
t1, f1 = np.loadtxt(file1, unpack=True, usecols=(0, 1))
t2, f2 = np.loadtxt(file2, unpack=True, usecols=(0, 1))
# do some fitting to determine temps
t1min = np.min(t1)
t2min = np.min(t2)
t1max = np.max(t1)
t2max = np.max(t2)
tmin = np.min([t1min, t2min])
tmax = np.max([t1max, t2max])
# warn about extrapolation
if not t1min == t2min:
warnings.warn(f"free energy is being extrapolated!")
if not t1max == t2max:
warnings.warn(f"free energy is being extrapolated!")
# now fit
f1fit = np.polyfit(t1, f1, fit_order)
f2fit = np.polyfit(t2, f2, fit_order)
# reevaluate over the new range
fit_t = np.arange(tmin, tmax + 1, 1)
fit_f1 = np.polyval(f1fit, fit_t)
fit_f2 = np.polyval(f2fit, fit_t)
# now evaluate the intersection temp
arg = np.argsort(np.abs(fit_f1 - fit_f2))[0]
transition_temp = fit_t[arg]
# warn if the temperature is shady
if np.abs(transition_temp - tmin) < 1e-3:
warnings.warn("It is likely there is no intersection of free energies")
elif np.abs(transition_temp - tmax) < 1e-3:
warnings.warn("It is likely there is no intersection of free energies")
# plot
if plot:
c1lo = "#ef9a9a"
c1hi = "#b71c1c"
c2lo = "#90caf9"
c2hi = "#0d47a1"
plt.plot(fit_t, fit_f1, color=c1lo, label=f"{folder1} fit")
plt.plot(fit_t, fit_f2, color=c2lo, label=f"{folder2} fit")
plt.plot(t1, f1, color=c1hi, label=folder1, ls="dashed")
plt.plot(t2, f2, color=c2hi, label=folder2, ls="dashed")
plt.axvline(transition_temp, ls="dashed", c="#37474f")
plt.ylabel("Free energy (eV/atom)")
plt.xlabel("Temperature (K)")
plt.legend(frameon=False)
return transition_temp