Skip to content

Double ufm #257

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
srmnitc merged 2 commits into from
Jul 8, 2026
Merged

Double ufm #257

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bff7b6e
add option for multilevel liquid switching
srmnitc Jun 25, 2026
303b130
update docs
srmnitc Jun 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 14 additions & 1 deletion calphy/input.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,7 +137,20 @@ def _extract_elements_from_pair_coeff(pair_coeff_string):

class UFMP(BaseModel, title="UFM potential input options"):
p: Annotated[float, Field(default=50.0)]
sigma: Annotated[float, Field(default=1.5)]
# sigma may be a scalar (single-component UFM reference, original behaviour) or a
# mapping of per-element-pair length scales for the two-leg UFM reference path, e.g.
# sigma: {H_H: 0.9, O_O: 2.8, H_O: 0.9}
# Keys are "<elementA>_<elementB>" (order-insensitive). Missing cross terms are
# filled by LAMMPS geometric mixing.
sigma: Annotated[Union[float, Dict[str, float]], Field(default=1.5)]
# When set, activates the two-leg UFM path: the system is first switched from the
# real potential to the (possibly multi-component) UFM at `sigma`, then from that
# UFM to a single-component UFM at `single_sigma` whose absolute free energy is
# known analytically (find_fe). Leaving this None reproduces the original
# single-leg behaviour exactly.
single_sigma: Annotated[Optional[float], Field(default=None)]
# p for the single-component endpoint; defaults to `p` when not given.
single_p: Annotated[Optional[float], Field(default=None)]


class MonteCarlo(
Expand Down
12 changes: 9 additions & 3 deletions calphy/integrators.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ def integrate_path(
return w, q, flambda


def find_w(mainfolder, calc, full=False, solid=True):
def find_w(mainfolder, calc, full=False, solid=True, prefix=""):
"""
Integrate the irreversible work and dissipation for independent simulations

Expand All @@ -154,6 +154,11 @@ def find_w(mainfolder, calc, full=False, solid=True):
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
Expand All @@ -168,11 +173,12 @@ def find_w(mainfolder, calc, full=False, solid=True):
ws = []
qs = []

tag = ("_%s" % prefix) if prefix else ""
for i in range(calc.n_iterations):
fwdfilestring = "forward_%d.dat" % (i + 1)
fwdfilestring = "forward%s_%d.dat" % (tag, i + 1)
fwdfilename = os.path.join(mainfolder, fwdfilestring)

bkdfilestring = "backward_%d.dat" % (i + 1)
bkdfilestring = "backward%s_%d.dat" % (tag, i + 1)
bkdfilename = os.path.join(mainfolder, bkdfilestring)

w, q, flambda = integrate_path(
Expand Down
218 changes: 194 additions & 24 deletions calphy/liquid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -306,6 +306,14 @@ def run_integration(self, iteration=1):
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
# ---------------------------------------------------------------
Expand All @@ -323,10 +331,13 @@ def run_integration(self, iteration=1):

for command in ph.hybrid_pair_coeff_commands(self.calc):
lmp.command(command)
lmp.command(
"pair_coeff * * ufm %f %f"
% (self.eps, self.calc.uhlenbeck_ford_model.sigma)
)
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
Expand Down Expand Up @@ -361,8 +372,8 @@ def run_integration(self, iteration=1):
lmp.command("fix_modify f2 temp Tcm")

lmp.command(
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" screen no file forward_%d.dat'
% iteration
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" screen no file %s'
% (leg1_fwd % iteration)
)
lmp.command("run %d" % self.calc._n_switching_steps)

Expand All @@ -378,10 +389,13 @@ def run_integration(self, iteration=1):
# ---------------------------------------------------------------

lmp.command("pair_style ufm %f" % self.ufm_cutoff)
lmp.command(
"pair_coeff * * %f %f"
% (self.eps, self.calc.uhlenbeck_ford_model.sigma)
)
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")
Expand Down Expand Up @@ -420,10 +434,13 @@ def run_integration(self, iteration=1):

for command in ph.hybrid_pair_coeff_commands(self.calc):
lmp.command(command)
lmp.command(
"pair_coeff * * ufm %f %f"
% (self.eps, self.calc.uhlenbeck_ford_model.sigma)
)
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
Expand Down Expand Up @@ -452,8 +469,8 @@ def run_integration(self, iteration=1):
lmp.command("fix_modify f2 temp Tcm")

lmp.command(
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" screen no file backward_%d.dat'
% iteration
'fix f3 all print 1 "${dU1} ${dU2} ${flambda}" screen no file %s'
% (leg1_bkd % iteration)
)
lmp.command("run %d" % self.calc._n_switching_steps)

Expand All @@ -464,6 +481,12 @@ def run_integration(self, iteration=1):
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)

if self.calc.script_mode:
file = os.path.join(self.simfolder, "integration.lmp")
lmp.write(file)
Expand All @@ -478,6 +501,124 @@ def run_integration(self, iteration=1):
except OSError as e:
self.logger.warning(f"Failed to rename log file: {e}")

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}" 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}" screen no file backward_leg2_%d.dat'
% iteration
)
lmp.command("run %d" % self.calc._n_switching_steps)
teardown()

def thermodynamic_integration(self):
"""
Calculate free energy after integration step
Expand All @@ -495,15 +636,44 @@ def thermodynamic_integration(self):
Calculates the final work, energy dissipation and free energy by
matching with UFM model
"""
w, q, qerr = find_w(self.simfolder, self.calc, full=True, solid=False)
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,
)
# 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(
Expand Down
Loading
Loading