"""
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 numpy as np
import yaml
import os
from calphy.integrators import *
import calphy.helpers as ph
import calphy.phase as cph
from calphy.errors import *
[docs]class Liquid(cph.Phase):
"""
Class for free energy calculation with liquid as the reference state
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):
"""
Set up class
"""
# call base class
super().__init__(
calculation=calculation,
simfolder=simfolder,
log_to_screen=log_to_screen,
)
[docs] def rattle_structure(self, lmp):
"""
Disorder the structure using random displacements followed by a
controlled NVT cool-down before liquid equilibration.
This is a lightweight alternative to :meth:`melt_structure`:
"""
self.logger.info(
"Rattling structure: velocities at 2*thigh=%f, NVT cool-down to T=%f",
self.calc._temperature_high,
self.calc._temperature,
)
lmp.command(
"displace_atoms all random 0.1 0.1 0.1 %d" % np.random.randint(1, 10000)
)
lmp.command(
"velocity all create %f %d"
% (2.0 * self.calc._temperature_high, np.random.randint(1, 10000))
)
lmp.command(
"fix nh_rattle all nvt temp %f %f %f"
% (
self.calc._temperature_high,
self.calc._temperature,
self.calc.md.thermostat_damping[1],
)
)
lmp.command("run %d" % int(self.calc.md.n_small_steps))
lmp.command("unfix nh_rattle")
[docs] def melt_structure(self, lmp):
""" """
if self.calc._fix_lattice and self.calc.melting_cycle:
raise ValueError(
"Cannot fix lattice and melt structure (set to False) at the same time"
)
melted = False
# this is the multiplier for thigh to try melting routines
for thmult in np.arange(1.0, 2.0, 0.1):
trajfile = os.path.join(self.simfolder, "traj.melt")
if os.path.exists(trajfile):
os.remove(trajfile)
self.logger.info(
"Starting melting cycle with thigh temp %f, factor %f"
% (self.calc._temperature_high, thmult)
)
factor = (self.calc._temperature_high / self.calc._temperature) * thmult
lmp.command(
"velocity all create %f %d"
% (self.calc._temperature * factor, np.random.randint(1, 10000))
)
self.fix_nose_hoover(lmp, temp_start_factor=factor, temp_end_factor=factor)
lmp.command("run %d" % int(self.calc.md.n_small_steps))
self.unfix_nose_hoover(lmp)
self.dump_current_snapshot(lmp, "traj.melt")
# we have to check if the structure melted
solids = ph.find_solid_fraction(os.path.join(self.simfolder, "traj.melt"))
self.logger.info("fraction of solids found: %f", solids / self.natoms)
if solids / self.natoms < self.calc.tolerance.liquid_fraction:
melted = True
break
# if melting cycle is over and still not melted, raise error
if not melted:
self.lammps_close(lmp=lmp)
lmp.rotate_logs("melting")
raise SolidifiedError(
"Liquid system did not melt, maybe try a higher thigh temperature."
)
[docs] def run_averaging(self):
"""
Run averaging routine
Parameters
----------
None
Returns
-------
None
Notes
-----
Run averaging routine using LAMMPS. Starting from the initial lattice two different routines can
be followed:
If pressure is specified, MD simulations are run until the pressure converges within the given
threshold value.
If `fix_lattice` option is True, then the input structure is used as it is and the corresponding pressure
is calculated.
At the end of the run, the averaged box dimensions are calculated.
"""
# create lammps object
lmp = ph.create_object(self.calc, self.simfolder)
lmp = ph.set_pair_style(lmp, self.calc)
# set up structure
lmp = ph.create_structure(lmp, self.calc)
# set up potential
lmp = ph.set_pair_coeff(lmp, self.calc)
lmp = ph.set_mass(lmp, self.calc)
# Melt regime for the liquid
lmp.command(
"velocity all create %f %d"
% (self.calc._temperature_high, np.random.randint(1, 10000))
)
# add some computes
lmp.command("variable mvol equal vol")
lmp.command("variable mlx equal lx")
lmp.command("variable mly equal ly")
lmp.command("variable mlz equal lz")
lmp.command("variable mpress equal press")
lmp.command("variable mpe equal pe/atoms")
lmp.command("variable metotal equal etotal/atoms")
lmp.command("variable mtemp equal temp")
# Disorder the structure before equilibration
if self.calc.melting_cycle:
self.melt_structure(lmp)
else:
self.rattle_structure(lmp)
if not self.calc._fix_lattice:
# now assign correct temperature and equilibrate
self.run_zero_pressure_equilibration(lmp)
# equilibration-frame dump (post warm-up; no-op unless
# n_print_steps_equilibration > 0)
self.start_equilibration_dump(lmp)
# converge pressure
self.run_pressure_convergence(lmp)
else:
self.start_equilibration_dump(lmp)
self.run_constrained_pressure_convergence(lmp)
# check melted error
self.stop_equilibration_dump(lmp)
self.dump_current_snapshot(lmp, "traj.equilibration_stage1.dat")
self.check_if_solidfied(lmp, "traj.equilibration_stage1.dat")
self.dump_current_snapshot(lmp, "traj.equilibration_stage2.dat")
lmp = ph.write_data(lmp, "conf.equilibration.data")
self.lammps_close(lmp=lmp)
lmp.rotate_logs("averaging")
[docs] def run_integration(self, iteration=1):
"""
Run integration routine
Parameters
----------
iteration : int, optional
iteration number for running independent iterations
Returns
-------
None
Notes
-----
Run the integration routine where the initial and final systems are connected using
the lambda parameter. See algorithm 4 in publication.
"""
lmp = ph.create_object(self.calc, self.simfolder)
# Adiabatic switching parameters.
lmp.command("variable li equal 1.0")
lmp.command("variable lf equal 0.0")
lmp = ph.set_pair_style(lmp, self.calc)
# read in the conf file
# conf = os.path.join(self.simfolder, "conf.equilibration.dump")
conf = os.path.join(self.simfolder, "conf.equilibration.data")
lmp = ph.read_data(lmp, conf)
# set hybrid ufm and normal potential
# lmp = ph.set_hybrid_potential(lmp, self.options, self.eps)
lmp = ph.set_pair_coeff(lmp, self.calc)
lmp = ph.set_mass(lmp, self.calc)
# remap the box to get the correct pressure
lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)
lmp.command("fix f1 all nve")
lmp.command(
"fix f2 all langevin %f %f %f %d zero yes"
% (
self.calc._temperature,
self.calc._temperature,
self.calc.md.thermostat_damping[1],
np.random.randint(1, 10000),
)
)
lmp.command("run %d" % self.calc.n_equilibration_steps)
lmp.command("unfix f1")
lmp.command("unfix f2")
# Output file naming. For the original single-leg path the files are
# forward_%d.dat / backward_%d.dat (unchanged). For the two-leg path this
# first leg (real -> multi-component UFM) writes *_leg1_%d.dat and the
# second leg (multi-component UFM -> single-component UFM) writes
# *_leg2_%d.dat.
leg1_fwd = "forward_leg1_%d.dat" if self._is_two_leg else "forward_%d.dat"
leg1_bkd = "backward_leg1_%d.dat" if self._is_two_leg else "backward_%d.dat"
# ---------------------------------------------------------------
# FWD cycle
# ---------------------------------------------------------------
lmp.command("variable flambda equal ramp(${li},${lf})")
lmp.command("variable blambda equal 1.0-v_flambda")
lmp.command(
ph.scaled_pair_style_command(
self.calc,
["v_flambda"],
extra_terms=["v_blambda ufm %f" % self.ufm_cutoff],
)
)
for command in ph.hybrid_pair_coeff_commands(self.calc):
lmp.command(command)
for command in self.ufm_pair_coeff_commands(
self.eps,
self.calc.uhlenbeck_ford_model.sigma,
self._ufm_sigma_by_type,
substyle="ufm",
):
lmp.command(command)
compute_commands, real_energy, compute_ids = ph.real_pair_compute_commands(
self.calc
)
for command in compute_commands:
lmp.command(command)
lmp.command("compute c2 all pair ufm")
lmp.command("variable step equal step")
lmp.command("variable dU1 equal (%s)/atoms" % real_energy)
lmp.command("variable dU2 equal c_c2/atoms")
lmp.command("thermo_style custom step v_dU1 v_dU2")
lmp.command("thermo 1000")
lmp.command(
"velocity all create %f %d mom yes rot yes dist gaussian"
% (self.calc._temperature, np.random.randint(1, 10000))
)
lmp.command("fix f1 all nve")
lmp.command(
"fix f2 all langevin %f %f %f %d zero yes"
% (
self.calc._temperature,
self.calc._temperature,
self.calc.md.thermostat_damping[1],
np.random.randint(1, 10000),
)
)
lmp.command("compute Tcm all temp/com")
lmp.command("fix_modify f2 temp Tcm")
lmp.command(
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" '
'title "# dU_sys[eV/atom] dU_ref[eV/atom] lambda" '
"screen no file %s"
% (leg1_fwd % iteration)
)
lmp.command("run %d" % self.calc._n_switching_steps)
lmp.command("unfix f1")
lmp.command("unfix f2")
lmp.command("unfix f3")
for compute_id in compute_ids:
lmp.command("uncompute %s" % compute_id)
lmp.command("uncompute c2")
# ---------------------------------------------------------------
# EQBRM cycle
# ---------------------------------------------------------------
lmp.command("pair_style ufm %f" % self.ufm_cutoff)
for command in self.ufm_pair_coeff_commands(
self.eps,
self.calc.uhlenbeck_ford_model.sigma,
self._ufm_sigma_by_type,
substyle="",
):
lmp.command(command)
lmp.command("thermo_style custom step pe")
lmp.command("thermo 1000")
lmp.command("fix f1 all nve")
lmp.command(
"fix f2 all langevin %f %f %f %d zero yes"
% (
self.calc._temperature,
self.calc._temperature,
self.calc.md.thermostat_damping[1],
np.random.randint(1, 10000),
)
)
lmp.command("fix_modify f2 temp Tcm")
lmp.command("run %d" % self.calc.n_equilibration_steps)
lmp.command("unfix f1")
lmp.command("unfix f2")
# ---------------------------------------------------------------
# BKD cycle
# ---------------------------------------------------------------
lmp.command("variable flambda equal ramp(${lf},${li})")
lmp.command("variable blambda equal 1.0-v_flambda")
lmp.command(
ph.scaled_pair_style_command(
self.calc,
["v_flambda"],
extra_terms=["v_blambda ufm %f" % self.ufm_cutoff],
)
)
for command in ph.hybrid_pair_coeff_commands(self.calc):
lmp.command(command)
for command in self.ufm_pair_coeff_commands(
self.eps,
self.calc.uhlenbeck_ford_model.sigma,
self._ufm_sigma_by_type,
substyle="ufm",
):
lmp.command(command)
compute_commands, real_energy, compute_ids = ph.real_pair_compute_commands(
self.calc
)
for command in compute_commands:
lmp.command(command)
lmp.command("compute c2 all pair ufm")
lmp.command("variable step equal step")
lmp.command("variable dU1 equal (%s)/atoms" % real_energy)
lmp.command("variable dU2 equal c_c2/atoms")
lmp.command("thermo_style custom step v_dU1 v_dU2")
lmp.command("thermo 1000")
lmp.command("fix f1 all nve")
lmp.command(
"fix f2 all langevin %f %f %f %d zero yes"
% (
self.calc._temperature,
self.calc._temperature,
self.calc.md.thermostat_damping[1],
np.random.randint(1, 10000),
)
)
lmp.command("fix_modify f2 temp Tcm")
lmp.command(
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" '
'title "# dU_sys[eV/atom] dU_ref[eV/atom] lambda" '
"screen no file %s"
% (leg1_bkd % iteration)
)
lmp.command("run %d" % self.calc._n_switching_steps)
lmp.command("unfix f1")
lmp.command("unfix f2")
lmp.command("unfix f3")
for compute_id in compute_ids:
lmp.command("uncompute %s" % compute_id)
lmp.command("uncompute c2")
# ---------------------------------------------------------------
# LEG 2 (two-leg path only): multi-component UFM -> single-component UFM
# ---------------------------------------------------------------
if self._is_two_leg:
self._run_leg2(lmp, iteration)
# close object
self.lammps_close(lmp=lmp)
lmp.rotate_logs("integration")
def _run_leg2(self, lmp, iteration):
"""
Second leg of the two-leg UFM reference path.
Switches the system from the (multi-component) UFM reference used in leg 1
to a single-component UFM reference at ``single_sigma`` whose absolute free
energy is known analytically. Writes forward_leg2_%d.dat / backward_leg2_%d.dat.
At entry the simulation may be in any state (leg-1 backward leaves it on the
real potential); this method re-establishes pure multi-component UFM, briefly
equilibrates, then performs the forward and backward switching runs. The two
UFM end states are distinguished inside ``hybrid/scaled`` as substyles
``ufm 1`` (multi-component, leg-1 sigmas) and ``ufm 2`` (single-component,
single_sigma).
Both substyles use the SAME eps so that only the length scale changes along
the leg. The single-component eps used for the analytic free energy
(``self.single_eps``) is applied at the find_fe stage, not here.
"""
T = self.calc._temperature
tdamp = self.calc.md.thermostat_damping[1]
# re-establish pure multi-component UFM and equilibrate
lmp.command("pair_style ufm %f" % self.ufm_cutoff)
for command in self.ufm_pair_coeff_commands(
self.eps, self.calc.uhlenbeck_ford_model.sigma,
self._ufm_sigma_by_type, substyle="",
):
lmp.command(command)
lmp.command("thermo_style custom step pe")
lmp.command("thermo 1000")
lmp.command("fix f1 all nve")
lmp.command(
"fix f2 all langevin %f %f %f %d zero yes"
% (T, T, tdamp, np.random.randint(1, 10000))
)
lmp.command("fix_modify f2 temp Tcm")
lmp.command("run %d" % self.calc.n_equilibration_steps)
lmp.command("unfix f1")
lmp.command("unfix f2")
# helper to set up the dual-ufm hybrid/scaled state for one direction
def setup_dual_ufm(forward):
if forward:
lmp.command("variable flambda equal ramp(${li},${lf})")
else:
lmp.command("variable flambda equal ramp(${lf},${li})")
lmp.command("variable blambda equal 1.0-v_flambda")
# flambda scales the multi-component UFM (ufm 1); blambda the single (ufm 2)
lmp.command(
"pair_style hybrid/scaled v_flambda ufm %f v_blambda ufm %f"
% (self.ufm_cutoff, self.single_ufm_cutoff)
)
for command in self.ufm_pair_coeff_commands(
self.eps, self.calc.uhlenbeck_ford_model.sigma,
self._ufm_sigma_by_type, substyle="ufm 1",
):
lmp.command(command)
for command in self.ufm_pair_coeff_commands(
self.eps, self.single_sigma, None, substyle="ufm 2",
):
lmp.command(command)
lmp.command("compute c1 all pair ufm 1")
lmp.command("compute c2 all pair ufm 2")
lmp.command("variable dU1 equal c_c1/atoms")
lmp.command("variable dU2 equal c_c2/atoms")
lmp.command("thermo_style custom step v_dU1 v_dU2")
lmp.command("thermo 1000")
lmp.command("fix f1 all nve")
lmp.command(
"fix f2 all langevin %f %f %f %d zero yes"
% (T, T, tdamp, np.random.randint(1, 10000))
)
lmp.command("fix_modify f2 temp Tcm")
def teardown():
lmp.command("unfix f1")
lmp.command("unfix f2")
lmp.command("unfix f3")
lmp.command("uncompute c1")
lmp.command("uncompute c2")
# FWD: multi -> single
setup_dual_ufm(forward=True)
lmp.command(
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" '
'title "# dU_sys[eV/atom] dU_ref[eV/atom] lambda" '
"screen no file forward_leg2_%d.dat"
% iteration
)
lmp.command("run %d" % self.calc._n_switching_steps)
teardown()
# EQBRM at pure single-component UFM
lmp.command("pair_style ufm %f" % self.single_ufm_cutoff)
for command in self.ufm_pair_coeff_commands(
self.eps, self.single_sigma, None, substyle="",
):
lmp.command(command)
lmp.command("thermo_style custom step pe")
lmp.command("thermo 1000")
lmp.command("fix f1 all nve")
lmp.command(
"fix f2 all langevin %f %f %f %d zero yes"
% (T, T, tdamp, np.random.randint(1, 10000))
)
lmp.command("fix_modify f2 temp Tcm")
lmp.command("run %d" % self.calc.n_equilibration_steps)
lmp.command("unfix f1")
lmp.command("unfix f2")
# BKD: single -> multi
setup_dual_ufm(forward=False)
lmp.command(
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" '
'title "# dU_sys[eV/atom] dU_ref[eV/atom] lambda" '
"screen no file backward_leg2_%d.dat"
% iteration
)
lmp.command("run %d" % self.calc._n_switching_steps)
teardown()
[docs] def thermodynamic_integration(self):
"""
Calculate free energy after integration step
Parameters
----------
None
Returns
-------
None
Notes
-----
Calculates the final work, energy dissipation and free energy by
matching with UFM model
"""
if self._is_two_leg:
# Two-leg UFM path:
# leg 1: real potential -> multi-component UFM (work w1)
# leg 2: multi-component UFM -> single-component UFM at single_sigma (work w2)
# The analytic reference is the single-component UFM at single_sigma/single_p,
# and the total switching work is w1 + w2.
w1, q1, qerr1 = find_w(
self.simfolder, self.calc, full=True, solid=False, prefix="leg1"
)
w2, q2, qerr2 = find_w(
self.simfolder, self.calc, full=True, solid=False, prefix="leg2"
)
w = w1 + w2
q = q1 + q2
qerr = np.sqrt(qerr1**2 + qerr2**2)
# store legs for diagnostics
self.w_leg1 = w1
self.w_leg2 = w2
single_p = self.calc.uhlenbeck_ford_model.single_p
if single_p is None:
single_p = self.calc.uhlenbeck_ford_model.p
f1 = get_uhlenbeck_ford_fe(
self.calc._temperature,
self.rho,
single_p,
self.calc.uhlenbeck_ford_model.single_sigma,
)
else:
w, q, qerr = find_w(self.simfolder, self.calc, full=True, solid=False)
# TODO: Hardcoded UFM parameters - enable option to change
f1 = get_uhlenbeck_ford_fe(
self.calc._temperature,
self.rho,
self.calc.uhlenbeck_ford_model.p,
self.calc.uhlenbeck_ford_model.sigma,
)
# Get ideal gas fe
f2 = get_ideal_gas_fe(
self.calc._temperature,
self.rho,
self.natoms,
[val["mass"] for key, val in self.calc._element_dict.items()],
[val["composition"] for key, val in self.calc._element_dict.items()],
)
self.ferr = qerr
self.fref = f1
self.fideal = f2
self.w = w
self.qdiss = q
# add pressure contribution if required
if self.calc._pressure != 0:
p = self.calc._pressure / EV_A3_TO_BAR
v = self.vol / self.natoms
self.pv = p * v
else:
self.pv = 0
# calculate final free energy
self.fe = self.fideal + self.fref - self.w + self.pv