Pcorr

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 2.2.0
+current_version = 2.2.1
 commit = True
 tag = True
 

CITATION.cff

@@ -61,4 +61,4 @@ abstract: >-
   examples illustrate the accuracy and efficiency of
   our implementation.
 license-url: 'https://github.com/ICAMS/calphy/blob/main/LICENSE'
-version: 2.2.0
+version: 2.2.1

calphy/__init__.py

@@ -28,7 +28,7 @@
 from calphy.alchemy import Alchemy
 from calphy.routines import MeltingTemp
 
-__version__ = "2.2.0"
+__version__ = "2.2.1"
 
 def addtest(a,b):
     return a+b

calphy/helpers.py

@@ -31,6 +31,26 @@
 from collections import Counter, defaultdict
 
 from ase.io import read, write
+from importlib.metadata import version as _dist_version
+
+#: pyscal 4 (still published as ``pyscal3``) is the minimum, as in
+#: pyproject.toml; conda environments and --no-deps installs bypass that pin.
+PYSCAL_MIN_MAJOR = 4
+
+
+def check_pyscal_version(installed=None):
+    """Raise ImportError unless the installed pyscal3 is at least version 4."""
+    if installed is None:
+        installed = _dist_version("pyscal3")
+    if int(installed.split(".")[0]) < PYSCAL_MIN_MAJOR:
+        raise ImportError(
+            "calphy needs pyscal3 >= %d.0.0 (the pyscal 4 API), found %s; "
+            "upgrade with: pip install 'pyscal3>=%d.0.0'"
+            % (PYSCAL_MIN_MAJOR, installed, PYSCAL_MIN_MAJOR)
+        )
+
+
+check_pyscal_version()
 
 import pyscal3
 from pyscal3 import Trajectory

calphy/input.py

@@ -58,7 +58,7 @@
     make_lattice,
 )
 
-__version__ = "2.2.0"
+__version__ = "2.2.1"
 
 
 def _check_equal(val):
@@ -156,6 +156,11 @@ def _extract_elements_from_pair_coeff(pair_coeff_string):
     "load_job": "job-state pickling was removed in calphy v2; rerun from the input file",
     "seed": "use md.seed -- one master seed now controls every stochastic step "
             "(the old quantum_thermal_bath seed was never actually applied)",
+    "lambda_schedule": "was removed in calphy 2.2.1; reversible scaling always "
+                       "ramps lambda linearly, as in de Koning et al., since only "
+                       "then does the barostat follow the scaled pressure lambda*P "
+                       "exactly (uniform_temperature only changed the sampling "
+                       "density, not the free energy)",
 }
 
 #: lazily built: ({field name -> [locations]}, {model class -> location label})
@@ -588,7 +593,6 @@ class Calculation(_StrictInput, title="Main input class"):
     n_print_steps_equilibration: Annotated[int, Field(default=0)]
     alchemy_coupling: Annotated[bool, Field(default=False)]
     n_iterations: Annotated[int, Field(default=1)]
-    lambda_schedule: Annotated[str, Field(default="linear")]
     equilibration_control: Annotated[Union[str, None], Field(default=None)]
     folder_prefix: Annotated[Union[str, None], Field(default=None)]
 
@@ -658,12 +662,6 @@ def _validate_all(self) -> "Input":
             if self.pair_mode not in ["overlay"]:
                 raise ValueError("pair_mode should be one of: overlay")
 
-        self.lambda_schedule = self.lambda_schedule.lower()
-        if self.lambda_schedule not in ["linear", "uniform_temperature"]:
-            raise ValueError(
-                "lambda_schedule must be one of: 'linear', 'uniform_temperature'"
-            )
-
         if self.pair_mode == "overlay":
             if self.pair_style is None or self.pair_coeff is None:
                 raise ValueError("pair_mode overlay requires pair_style and pair_coeff")

calphy/integrators.py

@@ -292,7 +292,9 @@ def integrate_rs(
         ws.append(w)
         es.append(e)
 
-    e_diss = np.min(es)
+    # worst iteration: one sweep that crossed a transition must not be hidden
+    # behind a clean one
+    e_diss = np.max(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)

calphy/phase.py

@@ -1435,13 +1435,47 @@ def _run_sweep(
         lmp.command("run               %d" % n_sweep)
         lmp.command("unfix             f3")
 
+    def _install_sweep_barostat(self, lmp, p_start, p_stop):
+        """
+        Replace the COM-constrained NPT fix ``f1`` by one whose pressure target
+        ramps ``p_start`` → ``p_stop`` over the run that follows.
+
+        The scaled system at (T0, λ) samples the real system at (T0/λ, P) only
+        if the barostat imposes the *scaled* pressure P_RS = λ·P (de Koning,
+        Antonelli and Yip, J. Chem. Phys. 115, 11025 (2001), Eq. 11), and
+        ``integrate_rs`` relies on exactly that when it integrates U + P·V at
+        the constant real pressure P.  ``fix npt`` can only ramp its target
+        linearly over a single run, so the constant-pressure equilibration fix
+        has to be swapped for a ramping one right before the sweep; the
+        ``fix_modify`` binding the COM-corrected temperature has to be reissued
+        for the new fix.  At ``pressure: 0`` the ramp is flat and the only
+        effect is a restart of the Nose-Hoover chain, as in ``pressure_scaling``.
+
+        Parameters
+        ----------
+        lmp : BaseRunner
+            Active LAMMPS runner with fix ``f1`` and compute ``tcm`` defined.
+        p_start, p_stop : float
+            Scaled pressure at the start and at the end of the sweep, in bar.
+        """
+        t0 = self.calc._temperature
+        lmp.command("unfix             f1")
+        lmp.command(
+            "fix               f1 all npt temp %f %f %f %s %f %f %f "
+            "fixedpoint ${xcm} ${ycm} ${zcm}"
+            % (t0, t0, self.calc.md.thermostat_damping[1],
+               self.iso, p_start, p_stop, self.calc.md.barostat_damping[1])
+        )
+        lmp.command("fix_modify        f1 temp tcm")
+
     def _reversible_scaling_forward(self, iteration: int = 1) -> None:
         """
         Perform the forward sweep of a reversible-scaling calculation.
 
         1. Short NPT warm start at T0 (see :meth:`_warm_start_steps`).
         2. COM-constrained equilibration at T0.
-        3. Forward sweep: λ 1 → T0/Tf.
+        3. Forward sweep: λ 1 → T0/Tf under the scaled pressure P_RS = λ·P,
+           i.e. the barostat ramps P → (T0/Tf)·P.
         4. Write ``conf.ts.forward_{iteration}.data`` for the backward sweep.
 
         Parameters
@@ -1542,36 +1576,15 @@ def _reversible_scaling_forward(self, iteration: int = 1) -> None:
             "forward sweep (iteration %d): COM-constrained equilibration done", iteration
         )
 
-        # ----------------------------------------------------------------
-        # Lambda schedule for the forward sweep.
-        #
-        # "linear" (default): lambda = ramp(li, lf) — simple linear
-        #   interpolation; LAMMPS ramp() resets automatically each run.
-        #
-        # "uniform_temperature": T_eq(s) = T0/lambda is linear in step
-        #   so every Kelvin bin gets the same number of MD samples.
-        #   Requires explicit step0 capture before each sweep.
-        # ----------------------------------------------------------------
-        lmp.command("variable         T0_rs equal %f" % t0)
-        if self.calc.lambda_schedule == "uniform_temperature":
-            lmp.command("variable         Nsweep equal %d" % self.calc._n_sweep_steps)
-            lmp.command("variable         Tf_rs equal %f" % tf)
-            # Capture the step at the START of the sweep so the formula is
-            # independent of any prior MD steps (no reset_timestep needed).
-            lmp.command("variable         step0 equal $(step)")
-            lmp.command(
-                "variable         flambda equal "
-                "v_T0_rs/(v_T0_rs+(v_Tf_rs-v_T0_rs)*(step-v_step0)/v_Nsweep)"
-            )
-            lmp.command(
-                "variable         blambda equal "
-                "v_T0_rs/(v_Tf_rs-(v_Tf_rs-v_T0_rs)*(step-v_step0)/v_Nsweep)"
-            )
-        else:  # "linear" (default)
-            lmp.command("variable         flambda equal ramp(${li},${lf})")
-            lmp.command("variable         blambda equal ramp(${lf},${li})")
-        lmp.command("variable         ftemp equal v_T0_rs/v_flambda")
-        lmp.command("variable         btemp equal v_T0_rs/v_blambda")
+        # ── Sweep barostat: P_RS ramps pi → pf = lf·pi ─────────────────────
+        if self.calc.npt:
+            self._install_sweep_barostat(lmp, pi, pf)
+
+        # λ ramps linearly in the step, like the target of the sweep
+        # barostat, so the barostat follows P_RS = λ·P exactly.  LAMMPS
+        # ramp() resets automatically each run.
+        lmp.command("variable         flambda equal ramp(${li},${lf})")
+        lmp.command("variable         blambda equal ramp(${lf},${li})")
 
         # Scaled Hamiltonian lambda*U from a single copy of the potential.
         # The earlier two-copy form, 1*U + (lambda-1)*U, gives the identical
@@ -1595,12 +1608,17 @@ def _reversible_scaling_forward(self, iteration: int = 1) -> None:
             )
             for combo in swap_combos:
                 self.logger.info("  swapping types %s ↔ %s", combo[0], combo[1])
+            # atom/swap accepts on the energy of the scaled potential, λ·ΔU,
+            # and the scaled system is canonical at the thermostat temperature
+            # T0: exp(-λ·ΔU / kT0) is the Boltzmann factor of the real system
+            # at T0/λ.  So the swap temperature is T0 throughout the sweep, not
+            # T0/λ, which would sample exp(-λ²·ΔU / kT0).
             for idx, (type1, type2) in enumerate(swap_combos):
                 lmp.command(
-                    "fix  swap%d all atom/swap %d %d %d ${ftemp} ke yes types %s %s"
+                    "fix  swap%d all atom/swap %d %d %d %f ke yes types %s %s"
                     % (idx, self.calc.monte_carlo.n_steps,
                        self.calc.monte_carlo.n_swaps,
-                       np.random.randint(1, 10000), type1, type2)
+                       np.random.randint(1, 10000), t0, type1, type2)
                 )
 
         if self.calc.n_print_steps > 0:
@@ -1664,8 +1682,9 @@ def _reversible_scaling_backward(self, iteration: int = 1) -> None:
 
         1. Load ``conf.ts.forward_{iteration}.data`` written by the forward
            sweep.
-        2. Middle equilibration at Tf.
-        3. Backward sweep: λ T0/Tf → 1.
+        2. Middle equilibration at Tf, i.e. at λ = T0/Tf and the scaled
+           pressure (T0/Tf)·P.
+        3. Backward sweep: λ T0/Tf → 1, the scaled pressure ramping back to P.
 
         Parameters
         ----------
@@ -1716,7 +1735,9 @@ def _reversible_scaling_backward(self, iteration: int = 1) -> None:
         # scaled potential.  Using a constant scaling variable (rather
         # than the ramp) keeps λ frozen at lf during this run.  A single
         # copy of the potential scaled by lambda is used throughout the
-        # reversible-scaling stage (see _reversible_scaling_forward).
+        # reversible-scaling stage (see _reversible_scaling_forward).  The
+        # barostat likewise holds the scaled pressure pf = lf·pi that the
+        # forward sweep ended at, see _install_sweep_barostat.
         lmp.command("variable          lambda_eq equal %f" % lf)
         lmp.command(ph.scaled_pair_style_command(self.calc, ["v_lambda_eq"]))
         for cmd in ph.hybrid_pair_coeff_commands(self.calc):
@@ -1731,7 +1752,7 @@ def _reversible_scaling_backward(self, iteration: int = 1) -> None:
                 "fix               f1 all npt temp %f %f %f %s %f %f %f "
                 "fixedpoint ${xcm} ${ycm} ${zcm}"
                 % (t0, t0, self.calc.md.thermostat_damping[1],
-                   self.iso, pi, pi, self.calc.md.barostat_damping[1])
+                   self.iso, pf, pf, self.calc.md.barostat_damping[1])
             )
         else:
             lmp.command(
@@ -1763,28 +1784,15 @@ def _reversible_scaling_backward(self, iteration: int = 1) -> None:
         else:
             self.check_if_solidfied(lmp, "traj.temp.dat")
 
+        # ── Sweep barostat: P_RS ramps pf = lf·pi → pi ─────────────────────
+        if self.calc.npt:
+            self._install_sweep_barostat(lmp, pf, pi)
+
         # ── Switch from the constant-λ scaled potential to the ramping
         # scaled potential for the backward sweep: define the lambda
         # variables, then re-install hybrid/scaled driven by blambda.
-        # T0_rs is needed by both schedules for ftemp/btemp.
-        lmp.command("variable         T0_rs equal %f" % t0)
-        if self.calc.lambda_schedule == "uniform_temperature":
-            lmp.command("variable         Nsweep equal %d" % self.calc._n_sweep_steps)
-            lmp.command("variable         Tf_rs equal %f" % tf)
-            lmp.command("variable         step0 equal $(step)")
-            lmp.command(
-                "variable         flambda equal "
-                "v_T0_rs/(v_T0_rs+(v_Tf_rs-v_T0_rs)*(step-v_step0)/v_Nsweep)"
-            )
-            lmp.command(
-                "variable         blambda equal "
-                "v_T0_rs/(v_Tf_rs-(v_Tf_rs-v_T0_rs)*(step-v_step0)/v_Nsweep)"
-            )
-        else:  # "linear"
-            lmp.command("variable         flambda equal ramp(${li},${lf})")
-            lmp.command("variable         blambda equal ramp(${lf},${li})")
-        lmp.command("variable         ftemp equal v_T0_rs/v_flambda")
-        lmp.command("variable         btemp equal v_T0_rs/v_blambda")
+        lmp.command("variable         flambda equal ramp(${li},${lf})")
+        lmp.command("variable         blambda equal ramp(${lf},${li})")
 
         lmp.command(ph.scaled_pair_style_command(self.calc, ["v_blambda"]))
         for cmd in ph.hybrid_pair_coeff_commands(self.calc):
@@ -1805,12 +1813,13 @@ def _reversible_scaling_backward(self, iteration: int = 1) -> None:
             )
             for combo in swap_combos:
                 self.logger.info("  swapping types %s ↔ %s", combo[0], combo[1])
+            # Swap temperature T0, as in the forward sweep.
             for idx, (type1, type2) in enumerate(swap_combos):
                 lmp.command(
-                    "fix  swap%d all atom/swap %d %d %d ${btemp} ke yes types %s %s"
+                    "fix  swap%d all atom/swap %d %d %d %f ke yes types %s %s"
                     % (idx, self.calc.monte_carlo.n_steps,
                        self.calc.monte_carlo.n_swaps,
-                       np.random.randint(1, 10000), type1, type2)
+                       np.random.randint(1, 10000), t0, type1, type2)
                 )
 
         if self.calc.n_print_steps > 0:
@@ -1996,12 +2005,14 @@ def scan_temperature_range(self):
         lmp.command("unfix             1")
 
         # ── Real-thermostat ramp T0 -> Tf, recording every step ─────────────
-        pf = (t0 / tf) * p0
+        # The real system is heated, so the barostat holds the real pressure
+        # p0 throughout; the scaled-pressure ramp belongs to the λ-scaled
+        # sweep of _reversible_scaling_forward only.
         lmp.command("variable          dU      equal pe/atoms")
         lmp.command(
             "fix               f2 all npt temp %f %f %f %s %f %f %f"
             % (t0, tf, self.calc.md.thermostat_damping[1],
-               self.iso, p0, pf, self.calc.md.barostat_damping[1])
+               self.iso, p0, p0, self.calc.md.barostat_damping[1])
         )
         scan_file = "prescan.forward.dat"
         lmp.command(
@@ -2155,8 +2166,11 @@ def temperature_scaling(self, iteration=1):
         tf = self.calc._temperature_stop
         li = 1
         lf = t0 / tf
+        # The real temperature is ramped by the thermostat and λ = T0/T is only
+        # the recorded label, so the real pressure p0 is held fixed throughout;
+        # a scaled-pressure ramp would sample the wrong isobar (it belongs to
+        # the λ-scaled Hamiltonian of reversible_scaling, not here).
         p0 = self.calc._pressure
-        pf = lf * p0
 
         # create lammps object
         lmp = ph.create_object(self.calc, self.simfolder)
@@ -2208,7 +2222,7 @@ def temperature_scaling(self, iteration=1):
                 self.calc.md.thermostat_damping[1],
                 self.iso,
                 p0,
-                pf,
+                p0,
                 self.calc.md.barostat_damping[1],
             )
         )
@@ -2234,8 +2248,8 @@ def temperature_scaling(self, iteration=1):
                 tf,
                 self.calc.md.thermostat_damping[1],
                 self.iso,
-                pf,
-                pf,
+                p0,
+                p0,
                 self.calc.md.barostat_damping[1],
             )
         )
@@ -2255,18 +2269,19 @@ def temperature_scaling(self, iteration=1):
         else:
             self.check_if_solidfied(lmp, "traj.temp.dat")
 
-        # start reverse loop
+        # start reverse loop: the thermostat ramps Tf -> T0, mirroring the
+        # forward sweep
         lmp.command("variable          lambda equal ramp(${lf},${li})")
 
         lmp.command(
             "fix               f2 all npt temp %f %f %f %s %f %f %f"
             % (
-                t0,
+                tf,
                 t0,
                 self.calc.md.thermostat_damping[1],
                 self.iso,
                 p0,
-                pf,
+                p0,
                 self.calc.md.barostat_damping[1],
             )
         )