"""
calphy: a Python library and command line interface for automated free
energy calculations.
Copyright 2021 (c) Sarath Menon^1, Yury Lysogorskiy^2, Ralf Drautz^2
^1: Max Planck Institut für Eisenforschung, Dusseldorf, Germany
^2: Ruhr-University Bochum, Bochum, Germany
calphy is published and distributed under the Academic Software License 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.
calphy API is published and distributed under the BSD 3-Clause "New" or "Revised" License
See the LICENSE FILE for more details.
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/yury.lysogorskiy@icams.rub.de
"""
import copy
import numpy as np
import os
import time
from mendeleev import element
import yaml
from calphy.input import read_inputfile
# import calphy.queuekernel as cq
from calphy.errors import *
import calphy.helpers as ph
from calphy.liquid import Liquid
from calphy.solid import Solid
from calphy.composition_transformation import CompositionTransformation
[docs]class MeltingTemp:
"""
Class for automated melting temperature calculation.
Parameters
----------
options : dict
dict of input options
kernel : int
the index of the calculation that should be run from
the list of calculations in the input file
simfolder : string
base folder for running calculations
"""
[docs] def __init__(self, calculation=None, simfolder=None, log_to_screen=False):
self.calc = calculation
self.simfolder = simfolder
self.log_to_screen = log_to_screen
self.dtemp = self.calc.melting_temperature.step
self.maxattempts = self.calc.melting_temperature.attempts
self.attempts = 0
self.calculations = []
self.get_trange()
self.arg = None
logfile = os.path.join(os.getcwd(), f"{self.calc.create_identifier()}.log")
self.logger = ph.prepare_log(logfile, screen=log_to_screen)
# Values used when melting_temperature has to switch the structural
# phase-stability checks back on for its sub-calculations.
DETECTION_SOLID_FRACTION = 0.7
DETECTION_LIQUID_FRACTION = 0.05
def _enable_phase_detection(self, calc):
"""
Force the structural phase-stability checks on for a sub-calculation.
``run_jobs`` signals "solid melted" / "liquid froze" purely through
MeltedError / SolidifiedError, and ``start_calculation`` walks the
temperature bracket on those signals alone. The measured solid
fraction is bounded to [0, 1], so the shipped defaults
(``solid_fraction = 0``, ``liquid_fraction = 1``) make both checks
unreachable -- the bracket would never be corrected and Tm would be
reported without ever verifying that the solid stayed solid and the
liquid stayed liquid.
The checks are therefore enabled here for this mode only, leaving the
global defaults alone. An explicit user setting is always respected.
Parameters
----------
calc : dict
Raw sub-calculation dict, mutated in place before it is parsed.
"""
tolerance = calc.setdefault("tolerance", {})
if not isinstance(tolerance, dict):
return
if tolerance.get("solid_fraction") is None:
tolerance["solid_fraction"] = self.DETECTION_SOLID_FRACTION
self.logger.warning(
"mode melting_temperature: enabling melt detection with "
"tolerance.solid_fraction = %g (the default of 0 makes the "
"check unreachable, and the temperature bracket is advanced "
"only when it fires). Set tolerance.solid_fraction "
"explicitly to override." % self.DETECTION_SOLID_FRACTION
)
if tolerance.get("liquid_fraction") is None:
tolerance["liquid_fraction"] = self.DETECTION_LIQUID_FRACTION
self.logger.warning(
"mode melting_temperature: enabling solidification detection "
"with tolerance.liquid_fraction = %g (the default of 1 makes "
"the check unreachable, and the temperature bracket is "
"advanced only when it fires). Set tolerance.liquid_fraction "
"explicitly to override." % self.DETECTION_LIQUID_FRACTION
)
[docs] def prepare_calcs(self):
"""
Prepare calculations list from given object
Parameters
----------
None
Returns
-------
None
"""
# here, we need to prepare a new calculation
# protocol, read in, modify, write a output
# read input again
calculations = {"calculations": []}
with open(self.calc.inputfile, "r") as fin:
data = yaml.safe_load(fin)
calc = data["calculations"][int(self.calc.kernel)]
calc["mode"] = "ts"
calc["temperature"] = [int(self.tmin), int(self.tmax)]
calc["reference_phase"] = "solid"
# Preserve n_iterations from the original melting_temperature calculation
if "n_iterations" in data["calculations"][int(self.calc.kernel)]:
calc["n_iterations"] = data["calculations"][int(self.calc.kernel)][
"n_iterations"
]
self._enable_phase_detection(calc)
calculations["calculations"].append(calc)
with open(self.calc.inputfile, "r") as fin:
data = yaml.safe_load(fin)
calc = data["calculations"][int(self.calc.kernel)]
calc["mode"] = "ts"
calc["temperature"] = [int(self.tmin), int(self.tmax)]
calc["reference_phase"] = "liquid"
# Preserve n_iterations from the original melting_temperature calculation
if "n_iterations" in data["calculations"][int(self.calc.kernel)]:
calc["n_iterations"] = data["calculations"][int(self.calc.kernel)][
"n_iterations"
]
self._enable_phase_detection(calc)
calculations["calculations"].append(calc)
outfile = f"{self.calc.create_identifier()}.{self.attempts}.yaml"
with open(outfile, "w") as fout:
yaml.safe_dump(calculations, fout)
# now read in again, which would allow for checking and so on
# one could do this smartly, and simply create from here.
self.calculations = read_inputfile(outfile)
[docs] def get_trange(self):
"""
Get temperature range for calculations
Parameters
----------
None
Returns
-------
None
"""
tmin = self.calc._temperature - self.dtemp
if tmin < 0:
tmin = 10
tmax = self.calc._temperature + self.dtemp
self.tmax = tmax
self.tmin = tmin
[docs] def run_jobs(self):
"""
Run calculations
Parameters
----------
None
Returns
-------
None
"""
self.prepare_calcs()
self.soljob = Solid(
calculation=self.calculations[0],
simfolder=self.calculations[0].create_folders(),
)
self.lqdjob = Liquid(
calculation=self.calculations[1],
simfolder=self.calculations[1].create_folders(),
)
# Propagate MeltingTemp file handlers to sub-job loggers so that
# all sub-job output also appears in melting_temperature.log
for handler in self.logger.handlers:
self.soljob.logger.addHandler(handler)
self.lqdjob.logger.addHandler(handler)
self.logger.info(
"Free energy of %s and %s phases will be calculated"
% (self.soljob.calc.lattice, self.lqdjob.calc.lattice)
)
self.logger.info("Temperature range of %f-%f" % (self.tmin, self.tmax))
self.logger.info("STATE: Temperature range of %f-%f K" % (self.tmin, self.tmax))
self.logger.info("Starting solid fe calculation")
try:
self.soljob = routine_fe(self.soljob)
except MeltedError:
self.logger.info("Solid phase melted")
return 2
self.logger.info("Starting solid reversible scaling run")
for i in range(self.soljob.calc.n_iterations):
try:
self.soljob.reversible_scaling(iteration=(i + 1))
except MeltedError:
self.logger.info("Solid system melted during reversible scaling run")
return 2
self.solres = self.soljob.integrate_reversible_scaling(
scale_energy=True, return_values=True
)
self.logger.info("Starting liquid fe calculation")
try:
self.lqdjob = routine_fe(self.lqdjob)
except SolidifiedError:
self.logger.info("Liquid froze")
return 3
self.logger.info("Starting liquid reversible scaling calculation")
for i in range(self.lqdjob.calc.n_iterations):
try:
self.lqdjob.reversible_scaling(iteration=(i + 1))
except SolidifiedError:
self.logger.info("Liquid froze during reversible scaling calculation")
return 3
self.lqdres = self.lqdjob.integrate_reversible_scaling(
scale_energy=True, return_values=True
)
self._report_sweep_quality()
def _report_sweep_quality(self):
"""
Repeat any dissipation warning from the two sub-calculations here.
The solid and liquid sweeps log into their own simfolders, which nobody
reads when the point of the mode is a single number at the end. A
sweep that dissipated heavily produces a free energy that is wrong by
an amount no averaging removes, and the crossing of two such curves is
wrong with it -- so the warning has to reach the melting-temperature
log, where it will actually be seen.
"""
offenders = [
(job.calc.reference_phase, job.ediss)
for job in (self.soljob, self.lqdjob)
if getattr(job, "ediss_high", False)
]
if not offenders:
return
detail = ", ".join(
"%s %.3e eV/atom" % (phase, value) for phase, value in offenders
)
self.logger.warning(
"Melting temperature is being extrapolated from a sweep that did "
"not stay reversible (%s; tolerance.dissipation = %.3e). The phase "
"very likely changed during the sweep, so treat the reported Tm as "
"unreliable rather than as an answer -- see the warnings in the "
"sub-calculation logs." % (detail, self.calc.tolerance.dissipation)
)
self.logger.warning("STATE: Tm unreliable, sweep dissipation too high")
[docs] def start_calculation(self):
"""
Start calculation
Parameters
----------
None
Returns
-------
None
"""
for i in range(100):
returncode = self.run_jobs()
if returncode == 3:
self.tmin = self.tmin + self.dtemp
self.tmax = self.tmax + self.dtemp
elif returncode == 2:
self.tmin = self.tmin - self.dtemp
if self.tmin < 0:
self.tmin = 0
self.tmax = self.tmax - self.dtemp
else:
return True
self.attempts += 1
if self.attempts > self.maxattempts:
raise ValueError("Maximum number of tries reached")
# Half-width of the finite-difference stencil used to take the local
# dF/dT on either side of the crossing, in samples. Capped against the
# sweep length so a short sweep cannot index past either end.
CROSSING_STENCIL = 50
def _crossing_error(self, arg, suberr):
"""
Propagate the free-energy uncertainty at the crossing into an
uncertainty on Tm.
The two curves cross at index ``arg``; converting a free-energy error
into a temperature error needs the rate at which they separate there,
i.e. the difference of their local slopes. That slope is taken from a
symmetric finite difference around ``arg``.
The stencil is clamped to the array: ``find_tm`` only rejects
``arg == 0`` and ``arg == len - 1``, so a crossing near either end
would otherwise read ``arg + 50`` past the end (IndexError, raised in
postprocessing after all the MD has been paid for) or ``arg - 50`` as
a negative index, which numpy silently wraps to the far end of the
sweep and turns the local slope into a chord across the whole
temperature range.
Parameters
----------
arg : int
Index of the crossing.
suberr : float
Combined free-energy uncertainty at the crossing, in eV/atom.
Returns
-------
tmerr : float
Uncertainty on Tm in K, or ``np.nan`` if the local slopes are too
close to distinguish (parallel curves give no crossing scale).
"""
n = len(self.solres[1])
half = max(1, min(self.CROSSING_STENCIL, n // 20))
lo = max(arg - half, 0)
hi = min(arg + half, n - 1)
if hi <= lo:
self.logger.warning(
"Sweep has too few samples (%d) to estimate a slope at the "
"crossing; reporting Tm without an error estimate." % n
)
return np.nan
def _slope(res):
dt = res[0][hi] - res[0][lo]
if dt == 0:
return np.nan
return (res[1][hi] - res[1][lo]) / dt
slope_diff = _slope(self.solres) - _slope(self.lqdres)
if not np.isfinite(slope_diff) or slope_diff == 0:
self.logger.warning(
"Solid and liquid free-energy curves are parallel at the "
"crossing; reporting Tm without an error estimate."
)
return np.nan
return suberr / slope_diff
[docs] def find_tm(self):
"""
Find melting temperature
Parameters
----------
None
Returns
-------
None
"""
for i in range(100):
arg = np.argsort(np.abs(self.solres[1] - self.lqdres[1]))[0]
self.arg = arg
if (arg == 0) or (arg == len(self.solres[1]) - 1):
self.logger.info(
"From calculation, melting temperature is not within the selected range."
)
self.logger.info("STATE: From calculation, Tm is not within range.")
if arg == len(self.solres[1]) - 1:
arg = 999
# the above is just a trick to extrapolate
# now here we need to find a guess value;
tpred = self.extrapolate_tm(arg)
# now we have to run calcs again
self.tmin = tpred - self.dtemp
if self.tmin < 0:
self.tmin = 0
self.tmax = tpred + self.dtemp
self.logger.info(
"Restarting calculation with predicted melting temperature +/- %f"
% self.dtemp
)
# self.logger.info('STATE: Restarting calculation with predicted melting temperature +/- %f'%self.dtemp)
self.start_calculation()
else:
self.calc_tm = self.solres[0][arg]
# get errors
suberr = np.sqrt(self.solres[2][arg] ** 2 + self.lqdres[2][arg] ** 2)
self.tmerr = self._crossing_error(arg, suberr)
return self.calc_tm, self.tmerr
self.attempts += 1
self.logger.info("Attempt incremented to %d" % self.attempts)
if self.attempts > self.maxattempts:
raise ValueError("Maximum number of tries reached")
[docs] def calculate_tm(self):
# do a first round of calculation
self.start_calculation()
tm, tmerr = self.find_tm()
self.logger.info("Found melting temperature = %.2f +/- %.2f K " % (tm, tmerr))
if self.calc._melting_temperature is not None:
self.logger.info(
"Experimental melting temperature = %.2f K "
% (self.calc._melting_temperature)
)
self.logger.info("STATE: Tm = %.2f K +/- %.2f K" % (tm, tmerr))
[docs]def routine_fe(job):
"""
Perform an FE calculation routine
"""
ts = time.time()
job.run_averaging()
te = time.time() - ts
job.logger.info("Averaging routine finished in %f s" % te)
# now run integration loops
for i in range(job.calc.n_iterations):
ts = time.time()
job.run_integration(iteration=(i + 1))
te = time.time() - ts
job.logger.info("Integration cycle %d finished in %f s" % (i + 1, te))
job.thermodynamic_integration()
job.submit_report()
job.clean_up()
return job
[docs]def routine_ts(job):
"""
Perform ts routine
"""
routine_fe(job)
# Optional pre-flight temperature-range scan: run a fast real-thermostat
# ramp T0 -> Tf, detect the onset of a phase transition, and (for
# mode='adapt') reduce the upper temperature so the production sweep stays
# in the single-phase region. Disabled by default (mode='none').
if job.calc.phase_transition_detection.mode != "none":
ts = time.time()
job.scan_temperature_range()
te = time.time() - ts
job.logger.info("Pre-flight temperature-range scan finished in %f s" % te)
# now do rev scale steps
for i in range(job.calc.n_iterations):
ts = time.time()
job.reversible_scaling(iteration=(i + 1))
te = time.time() - ts
job.logger.info("TS integration cycle %d finished in %f s" % (i + 1, te))
job.integrate_reversible_scaling(scale_energy=True)
job.clean_up()
return job
[docs]def routine_tscale(job):
"""
Perform tscale routine
"""
routine_fe(job)
# now do rev scale steps
for i in range(job.calc.n_iterations):
ts = time.time()
job.temperature_scaling(iteration=(i + 1))
te = time.time() - ts
job.logger.info("Temperature scaling cycle %d finished in %f s" % (i + 1, te))
job.integrate_reversible_scaling(scale_energy=False)
job.clean_up()
return job
[docs]def routine_pscale(job):
"""
Perform pscale routine
"""
routine_fe(job)
# now do rev scale steps
for i in range(job.calc.n_iterations):
ts = time.time()
job.pressure_scaling(iteration=(i + 1))
te = time.time() - ts
job.logger.info("Pressure scaling cycle %d finished in %f s" % (i + 1, te))
job.integrate_pressure_scaling()
job.clean_up()
return job
[docs]def routine_alchemy(job):
"""
Perform an FE calculation routine
"""
ts = time.time()
job.run_averaging()
te = time.time() - ts
job.logger.info("Averaging routine finished in %f s" % te)
# now run integration loops
for i in range(job.calc.n_iterations):
ts = time.time()
job.run_integration(iteration=(i + 1))
te = time.time() - ts
job.logger.info("Alchemy integration cycle %d finished in %f s" % (i + 1, te))
job.thermodynamic_integration()
job.submit_report()
job.clean_up()
return job
[docs]def routine_composition_scaling(job):
"""
Perform a compositional scaling routine
"""
# we set up comp scaling first
job.logger.info("Calculating composition scaling")
comp = CompositionTransformation(job.calc)
forward_swap_types, reverse_swap_types = comp.get_swap_types(
allow_all_swaps=job.calc.monte_carlo.allow_all_swaps
)
job.calc.monte_carlo.forward_swap_types = forward_swap_types
job.calc.monte_carlo.reverse_swap_types = reverse_swap_types
job.logger.info(f"Forward swap types: {forward_swap_types}")
job.logger.info(f"Reverse swap types: {reverse_swap_types}")
# update pair styles
res = comp.update_pair_coeff(job.calc.pair_coeff[0])
job.calc.pair_style.append(job.calc.pair_style[0])
job.calc._pair_style_with_options.append(job.calc._pair_style_with_options[0])
job.calc.pair_coeff[0] = res[0]
job.calc.pair_coeff.append(res[1])
job.logger.info("Update pair coefficients")
job.logger.info(f"pair coeff 1: {job.calc.pair_coeff[0]}")
job.logger.info(f"pair coeff 2: {job.calc.pair_coeff[1]}")
job.calc._pair_style_names.append(job.calc._pair_style_names[0])
job.logger.info("Update pair styles")
job.logger.info(f"pair style 1: {job.calc._pair_style_names[0]}")
job.logger.info(f"pair style 2: {job.calc._pair_style_names[1]}")
backup_element = job.calc.element.copy()
job.calc.element = comp.pair_list_old
# job.calc._ghost_element_count = len(comp.new_atomtype) - len()
# write new file out and update lattice
outfilename = ".".join([job.calc.lattice, "comp", "data"])
comp.write_structure(outfilename)
job.calc.lattice = outfilename
job.logger.info(f"Modified lattice written to {outfilename}")
# prepare mass change methods
# update and backup mass
job.logger.info(f"Original mass: {job.calc.mass}")
backup_mass = job.calc.mass.copy()
mass_dict = {key: val for (key, val) in zip(backup_element, backup_mass)}
target_masses = []
target_counts = []
ref_mass_list = []
for mdict in comp.transformation_list:
ref_mass_list.append(mass_dict[mdict["primary_element"]])
target_masses.append(mass_dict[mdict["secondary_element"]])
target_counts.append(mdict["count"])
if len(backup_mass) > 2:
job.logger.warning("Composition scaling is untested for more than 2 elements!")
if len(np.unique(ref_mass_list)) > 1:
job.logger.warning("More than one kind of transformation found! Stopping")
raise RuntimeError("More than one kind of transformation found! Stopping")
ref_mass = ref_mass_list[0]
# now replace mass
job.calc.mass = [ref_mass for x in range(len(job.calc.element))]
job.logger.info(f"Temporarily replacing mass: {job.calc.mass}")
# update fict elements if needed
# job.calc._totalelements = comp.maxtype
# now start cycle
ts = time.time()
job.run_averaging()
te = time.time() - ts
job.logger.info("Averaging routine finished in %f s" % te)
# now run integration loops
for i in range(job.calc.n_iterations):
ts = time.time()
job.run_integration(iteration=(i + 1))
te = time.time() - ts
job.logger.info("Alchemy integration cycle %d finished in %f s" % (i + 1, te))
job.thermodynamic_integration()
job.logger.info("performing mass rescaling")
job.logger.info(f"Ref. mass is {ref_mass}")
job.logger.info(f"Target masses are {target_masses}")
# read the file
mcorsum = job.mass_integration(ref_mass, target_masses, target_counts)
job.fe = job.fe - mcorsum
job.submit_report(
extra_dict={
"results": {
"mass_correction": float(mcorsum),
"entropy_contribution": float(comp.entropy_contribution),
}
}
)
job.clean_up()
return job