make ramp more flexible by defining lambda through a variable

doc/src/fix_alchemy.rst

@@ -8,17 +8,18 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group-ID alchemy
+   fix ID group-ID alchemy v_name
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * alchemy = style name of this fix command
+* v_name = variable with name that determines the :math:`\lambda_p` value
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
-   fix trans all alchemy
+   fix trans all alchemy v_ramp
 
 Description
 """""""""""
@@ -30,14 +31,17 @@ simulation between two topologies (i.e. the same number and positions of
 atoms, but differences in atom parameters like type, charge, bonds,
 angles and so on).  For this a :ref:`multi-partition run <partition>` is
 required with exactly two partitions.  During the MD run, the fix will
-will determine a factor :math:`\lambda_p` that will be linearly ramped
+will determine a factor, :math:`\lambda_p`, for each partition *p* that
-*down* from 1.0 to 0.0 for the *first* partition (*p=0*) and ramped *up*
+will be taken from an equal style or equivalent :doc:`variable
-from 0.0 to 1.0 for the *second* partition (*p=1*).  The forces used for
+<variable>`.  Typically, this variable would be chose to linearly ramp
-the propagation of the atoms will be the sum of the forces of the two
+*down* from 1.0 to 0.0 for the *first* partition (*p=0*) and linearly
-systems combined and scaled with their respective :math:`\lambda_p`
+ramp *up* from 0.0 to 1.0 for the *second* partition (*p=1*).  The
-factor.  This allows to perform transformations that are not easily
+forces used for the propagation of the atoms will be the sum of the
-possible with :doc:`pair style hybrid/scaled <pair_hybrid>`, :doc:`fix
+forces of the two systems combined and scaled with their respective
-adapt <fix_adapt>` or :doc:`fix adapt/fep <fix_adapt_fep>`.
+:math:`\lambda_p` factor.  This allows to perform transformations that
+are not easily possible with :doc:`pair style hybrid/scaled
+<pair_hybrid>`, :doc:`fix adapt <fix_adapt>` or :doc:`fix adapt/fep
+<fix_adapt_fep>`.
 
 Due to the specifics of the implementation, the initial geometry and
 dimensions of the system must be exactly the same and the fix will
@@ -56,6 +60,8 @@ copper/aluminum bronze.
 
 .. code-block:: LAMMPS
 
+   variable name world pure alloy
+
    create_box 2 box
    create_atoms 1 box
    pair_style eam/alloy
@@ -66,6 +72,15 @@ copper/aluminum bronze.
    if "${name} == alloy" then &
      "set type 1 type/fraction 2 0.05 6745234"
 
+   # define ramp variable to combine the two different partitions
+   if "${name} == pure" then             &
+     "variable ramp equal ramp(1.0,0.0)"    &
+   else                                      &
+      "variable ramp equal ramp(0.0,1.0)"
+
+   fix 2 all alchemy v_ramp
+
+
 The ``examples/PACKAGES/alchemy`` folder contains complete example
 inputs for this command.
 
@@ -78,7 +93,7 @@ No information about this fix is written to :doc:`binary restart files <restart>
 None of the :doc:`fix_modify <fix_modify>` options are relevant to this fix.
 
 This fix stores a global scalar (the current value of :math:`\lambda_p`)
-and a global vector or length 3 which contains the potential energy of
+and a global vector of length 3 which contains the potential energy of
 the first partition, the second partition and the combined value,
 respectively. The global scalar is unitless and "intensive", the vector
 is in :doc:`energy units <units>` and "extensive".  These values can be
@@ -86,14 +101,6 @@ used by any command that uses a global value from a fix as input.  See
 the :doc:`Howto output <Howto_output>` doc page for an overview of
 LAMMPS output options.
 
-The value of :math:`\lambda_p` is influenced by the *start/stop* keywords
-of the :doc:`run <run>` command.  Without them it will be ramped
-linearly from 1.0 to 0.0 (partition 1) and 0.0 to 1.0 (partition 2) 
-during the steps of a run, with
-*start/stop* keywords the ramp is from the *start* time step to the
-*stop* timestep. This allows to break down a simulation over multiple
-*run* commands or to continue transparently from a restart.
-
 This fix is not invoked during :doc:`energy minimization <minimize>`.
 
 Restrictions

examples/PACKAGES/alchemy/in.alloy

@@ -19,9 +19,15 @@ pair_coeff * * AlCu.eam.alloy Cu Al
 if "${name} == alloy" then &
    "set type 1 type/fraction 2 0.05 6745234"
 
+# define ramp variable to combine the two different partitions
+if "${name} == pure" then             &
+   "variable ramp equal ramp(1.0,0.0)"    &
+else                                      &
+   "variable ramp equal ramp(0.0,1.0)"
+
 velocity all create 5000.0 6567345
 fix 1 all nvt temp 2500.0 500.0 0.002
-fix 2 all alchemy
+fix 2 all alchemy v_ramp
 
 compute pressure all pressure/alchemy 2
 

examples/PACKAGES/alchemy/in.twowater

@@ -73,11 +73,17 @@ else &
 velocity all create 300.0 5463576
 timestep 0.2
 
+# define ramp variable to combine the two different partitions
+if "${name} == twowater" then             &
+   "variable ramp equal ramp(1.0,0.0)"    &
+else                                      &
+   "variable ramp equal ramp(0.0,1.0)"
+
 # since the trajectory and forces are kept identical through fix alchemy,
 # we can do fix npt simulations, but we must use the "mixed" pressure
 
 fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0
-fix transform all alchemy
+fix transform all alchemy v_ramp
 compute pressure all pressure/alchemy transform
 fix_modify integrate press pressure
 

examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.0

@@ -19,7 +19,7 @@ timestep 0.002
 create_atoms 1 box
 Created 864 atoms
   using lattice units in orthogonal box = (0 0 0) to (22.5 22.5 22.5)
-  create_atoms CPU = 0.000 seconds
+  create_atoms CPU = 0.001 seconds
 displace_atoms all random 0.3 0.3 0.3 57845645
 Displacing atoms ...
 pair_style eam/alloy
@@ -29,9 +29,13 @@ Reading eam/alloy potential file AlCu.eam.alloy with DATE: 2008-10-01
 # replace 5% of copper with aluminium on the second partition only
 if "${name} == alloy" then    "set type 1 type/fraction 2 0.05 6745234"
 
+# define ramp variable to combine the two different partitions
+if "${name} == pure" then                "variable ramp equal ramp(1.0,0.0)"    else                                         "variable ramp equal ramp(0.0,1.0)"
+variable ramp equal ramp(1.0,0.0)
+
 velocity all create 5000.0 6567345
 fix 1 all nvt temp 2500.0 500.0 0.002
-fix 2 all alchemy
+fix 2 all alchemy v_ramp
 
 compute pressure all pressure/alchemy 2
 
@@ -159,7 +163,7 @@ Per MPI rank memory allocation (min/avg/max) = 3.474 | 3.474 | 3.474 Mbytes
       9800   520.37758     -31000.6       -2828.0252      8.0039334     -2886.0741      58.048903      0.02          -2918.4921    
       9900   506.49368     -32223.744     -2830.1705      8.0039334     -2886.6706      56.500132      0.01          -2920.4741    
      10000   507.13597     -33930.476     -2834.3236      8.0039334     -2890.8954      56.571781      0             -2924.5423    
-Loop time of 19.2208 on 2 procs for 10000 steps with 864 atoms
+Loop time of 17.7279 on 2 procs for 10000 steps with 864 atoms
 
 
-Total wall time: 0:00:19
+Total wall time: 0:00:17

examples/PACKAGES/alchemy/log.25Feb23.alloy.g++.4.1

@@ -32,9 +32,13 @@ set type 1 type/fraction 2 0.05 6745234
 Setting atom values ...
   37 settings made for type/fraction
 
+# define ramp variable to combine the two different partitions
+if "${name} == pure" then                "variable ramp equal ramp(1.0,0.0)"    else                                         "variable ramp equal ramp(0.0,1.0)"
+variable ramp equal ramp(0.0,1.0)
+
 velocity all create 5000.0 6567345
 fix 1 all nvt temp 2500.0 500.0 0.002
-fix 2 all alchemy
+fix 2 all alchemy v_ramp
 
 compute pressure all pressure/alchemy 2
 
@@ -162,7 +166,7 @@ Per MPI rank memory allocation (min/avg/max) = 3.474 | 3.474 | 3.474 Mbytes
       9800   524.81105     -31000.6       -2860.6103      7.8067107     -2919.1537      58.543462      0.98          -2918.4921    
       9900   507.60008     -32223.744     -2864.192       7.8067107     -2920.8155      56.623552      0.99          -2920.4741    
      10000   521.28866     -33930.476     -2866.3918      7.8067107     -2924.5423      58.150534      1             -2924.5423    
-Loop time of 19.2208 on 2 procs for 10000 steps with 864 atoms
+Loop time of 17.7277 on 2 procs for 10000 steps with 864 atoms
 
 
-Total wall time: 0:00:19
+Total wall time: 0:00:17

examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.0

@@ -1,9 +1,9 @@
 LAMMPS (8 Feb 2023)
 Processor partition = 0
   using 1 OpenMP thread(s) per MPI task
-# Example for alchemical transformation of two water molecules into a hydronium and hydroxyl ion
+# Example for an alchemical transformation of two water molecules into a hydronium and hydroxyl ion
-# WARNING: This input is intended for demonstrating the method, only
+# WARNING: This input is intended for demonstrating the method only,
-# The force field parameters are made up and NOT suitable for production simulations.
+# the force field parameters are mostly made up and NOT suitable for production simulations.
 
 # set up different names for two partitions
 variable name world twowater twoions
@@ -59,7 +59,7 @@ group transform id 1:6
 create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.33
 Created 96 atoms
   using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5)
-  create_atoms CPU = 0.001 seconds
+  create_atoms CPU = 0.000 seconds
 
 # change topology and settings for the two states
 # we cannot simply create a different topology directly or
@@ -156,11 +156,15 @@ Finding 1-2 1-3 1-4 neighbors ...
 velocity all create 300.0 5463576
 timestep 0.2
 
+# define ramp variable to combine the two different partitions
+if "${name} == twowater" then                "variable ramp equal ramp(1.0,0.0)"    else                                         "variable ramp equal ramp(0.0,1.0)"
+variable ramp equal ramp(1.0,0.0)
+
 # since the trajectory and forces are kept identical through fix alchemy,
 # we can do fix npt simulations, but we must use the "mixed" pressure
 
 fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0
-fix transform all alchemy
+fix transform all alchemy v_ramp
 compute pressure all pressure/alchemy transform
 fix_modify integrate press pressure
 
@@ -377,21 +381,21 @@ Per MPI rank memory allocation (min/avg/max) = 7.535 | 7.535 | 7.535 Mbytes
      19800   295.68946      81474.342      2308.3249      1.0240743      2219.3041      89.020804      0.01           1910.5545    
      19900   306.4947       43488.052      2330.4958      1.0222707      2238.222       92.27385       0.005          1957.4368    
      20000   313.31679     -25133.284     -1161.6979      1.0163289     -1256.0256      94.327722      0             -1649.7551    
-Loop time of 11.7903 on 2 procs for 20000 steps with 102 atoms
+Loop time of 8.92653 on 2 procs for 20000 steps with 102 atoms
 
-Performance: 29.312 ns/day, 0.819 hours/ns, 1696.311 timesteps/s, 173.024 katom-step/s
+Performance: 38.716 ns/day, 0.620 hours/ns, 2240.513 timesteps/s, 228.532 katom-step/s
-99.3% CPU use with 2 MPI tasks x 1 OpenMP threads
+96.3% CPU use with 2 MPI tasks x 1 OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 6.8605     | 8.6071     | 10.354     |  59.5 | 73.00
+Pair    | 4.7448     | 5.9908     | 7.2367     |  50.9 | 67.11
-Bond    | 0.044148   | 0.044761   | 0.045375   |   0.3 |  0.38
+Bond    | 0.036126   | 0.03641    | 0.036695   |   0.1 |  0.41
-Neigh   | 0.10568    | 0.10573    | 0.10578    |   0.0 |  0.90
+Neigh   | 0.082247   | 0.082294   | 0.082342   |   0.0 |  0.92
-Comm    | 0.74919    | 2.4939     | 4.2386     | 110.5 | 21.15
+Comm    | 0.68467    | 1.9303     | 3.1759     |  89.7 | 21.62
-Output  | 0.0033838  | 0.0036358  | 0.0038879  |   0.4 |  0.03
+Output  | 0.0023377  | 0.0052347  | 0.0081317  |   4.0 |  0.06
-Modify  | 0.46825    | 0.47526    | 0.48227    |   1.0 |  4.03
+Modify  | 0.8194     | 0.82122    | 0.82303    |   0.2 |  9.20
-Other   |            | 0.05993    |            |       |  0.51
+Other   |            | 0.06029    |            |       |  0.68
 
 Nlocal:             51 ave          51 max          51 min
 Histogram: 2 0 0 0 0 0 0 0 0 0
@@ -405,4 +409,4 @@ Ave neighs/atom = 359.21569
 Ave special neighs/atom = 2
 Neighbor list builds = 181
 Dangerous builds = 0
-Total wall time: 0:00:11
+Total wall time: 0:00:08

examples/PACKAGES/alchemy/log.25Feb23.twowater.g++.4.1

@@ -1,9 +1,9 @@
 LAMMPS (8 Feb 2023)
 Processor partition = 1
   using 1 OpenMP thread(s) per MPI task
-# Example for alchemical transformation of two water molecules into a hydronium and hydroxyl ion
+# Example for an alchemical transformation of two water molecules into a hydronium and hydroxyl ion
-# WARNING: This input is intended for demonstrating the method, only
+# WARNING: This input is intended for demonstrating the method only,
-# The force field parameters are made up and NOT suitable for production simulations.
+# the force field parameters are mostly made up and NOT suitable for production simulations.
 
 # set up different names for two partitions
 variable name world twowater twoions
@@ -106,7 +106,7 @@ Finding 1-2 1-3 1-4 neighbors ...
      1 = max # of 1-3 neighbors
      1 = max # of 1-4 neighbors
      7 = max # of special neighbors
-  special bonds CPU = 0.000 seconds
+  special bonds CPU = 0.002 seconds
 create_bonds single/bond 2 3 4
 Finding 1-2 1-3 1-4 neighbors ...
   special bond factors lj:    0        0        0       
@@ -180,11 +180,15 @@ Setting atom values ...
 velocity all create 300.0 5463576
 timestep 0.2
 
+# define ramp variable to combine the two different partitions
+if "${name} == twowater" then                "variable ramp equal ramp(1.0,0.0)"    else                                         "variable ramp equal ramp(0.0,1.0)"
+variable ramp equal ramp(0.0,1.0)
+
 # since the trajectory and forces are kept identical through fix alchemy,
 # we can do fix npt simulations, but we must use the "mixed" pressure
 
 fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0
-fix transform all alchemy
+fix transform all alchemy v_ramp
 compute pressure all pressure/alchemy transform
 fix_modify integrate press pressure
 
@@ -402,21 +406,21 @@ Per MPI rank memory allocation (min/avg/max) = 7.535 | 7.535 | 7.535 Mbytes
      19800   295.68946      81474.342      1996.4567      1.0240743      1907.4359      89.020804      0.99           1910.5545    
      19900   306.4947       43488.052      2048.2997      1.0222707      1956.0258      92.27385       0.995          1957.4368    
      20000   313.31679     -25133.284     -1555.4274      1.0163289     -1649.7551      94.327722      1             -1649.7551    
-Loop time of 11.7903 on 2 procs for 20000 steps with 102 atoms
+Loop time of 8.92592 on 2 procs for 20000 steps with 102 atoms
 
-Performance: 29.312 ns/day, 0.819 hours/ns, 1696.310 timesteps/s, 173.024 katom-step/s
+Performance: 38.719 ns/day, 0.620 hours/ns, 2240.664 timesteps/s, 228.548 katom-step/s
-99.3% CPU use with 2 MPI tasks x 1 OpenMP threads
+96.2% CPU use with 2 MPI tasks x 1 OpenMP threads
 
 MPI task timing breakdown:
 Section |  min time  |  avg time  |  max time  |%varavg| %total
 ---------------------------------------------------------------
-Pair    | 6.8748     | 8.631      | 10.387     |  59.8 | 73.20
+Pair    | 4.7753     | 6.0676     | 7.3599     |  52.5 | 67.98
-Bond    | 0.04667    | 0.047105   | 0.04754    |   0.2 |  0.40
+Bond    | 0.036352   | 0.036909   | 0.037466   |   0.3 |  0.41
-Neigh   | 0.10858    | 0.10862    | 0.10866    |   0.0 |  0.92
+Neigh   | 0.085317   | 0.085383   | 0.085449   |   0.0 |  0.96
-Comm    | 0.74367    | 2.4997     | 4.2558     | 111.1 | 21.20
+Comm    | 0.72794    | 2.022      | 3.3161     |  91.0 | 22.65
-Output  | 0.0034405  | 0.0037033  | 0.003966   |   0.4 |  0.03
+Output  | 0.0023255  | 0.005156   | 0.0079865  |   3.9 |  0.06
-Modify  | 0.4338     | 0.4388     | 0.4438     |   0.8 |  3.72
+Modify  | 0.64567    | 0.64689    | 0.64812    |   0.2 |  7.25
-Other   |            | 0.06136    |            |       |  0.52
+Other   |            | 0.06196    |            |       |  0.69
 
 Nlocal:             51 ave          51 max          51 min
 Histogram: 2 0 0 0 0 0 0 0 0 0
@@ -430,4 +434,4 @@ Ave neighs/atom = 359.20588
 Ave special neighs/atom = 2.0196078
 Neighbor list builds = 181
 Dangerous builds = 0
-Total wall time: 0:00:11
+Total wall time: 0:00:08

src/REPLICA/fix_alchemy.cpp

@@ -18,11 +18,13 @@
 #include "compute.h"
 #include "domain.h"
 #include "error.h"
+#include "input.h"
 #include "memory.h"
 #include "modify.h"
 #include "respa.h"
 #include "universe.h"
 #include "update.h"
+#include "variable.h"
 
 #include <cstring>
 
@@ -33,8 +35,12 @@ using namespace FixConst;
 
 FixAlchemy::FixAlchemy(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), commbuf(nullptr)
 {
-  if (narg != 3) error->all(FLERR, "Incorrect number of arguments for fix alchemy");
+  if (narg != 4) error->all(FLERR, "Incorrect number of arguments for fix alchemy");
   if (universe->nworlds != 2) error->all(FLERR, "Must use exactly two partitions");
+  if (utils::strmatch(arg[3], "^v_"))
+    id_lambda = arg[3] + 2;
+  else
+    error->all(FLERR, "Must use variable as lambda argument to fix alchemy");
 
   lambda = epot[0] = epot[1] = epot[2] = 0.0;
   progress = 0;
@@ -129,6 +135,14 @@ void FixAlchemy::init()
   if (modify->get_fix_by_style("^alchemy").size() > 1)
     error->all(FLERR, "There may only one fix alchemy at a time");
 
+  ivar = input->variable->find(id_lambda.c_str());
+  if (ivar < 0)
+    error->universe_one(FLERR, fmt::format("Variable {} for fix alchemy does not exist", id_lambda));
+  if (!input->variable->equalstyle(ivar))
+    error->universe_one(FLERR,
+                        fmt::format("Variable {} for fix alchemy is invalid style", id_lambda));
+  lambda = input->variable->compute_equal(ivar);
+
   // synchronize box dimensions, determine if resync during run will be needed.
 
   synchronize_box(domain, samerank);
@@ -152,7 +166,10 @@ void FixAlchemy::setup(int vflag)
   }
 
   if (universe->me == 0) {
-    progress = static_cast<int>(100.0 - lambda * 100.0);
+    double delta = update->ntimestep - update->beginstep;
+    if ((delta != 0.0) && (update->beginstep != update->endstep))
+          delta /= update->endstep - update->beginstep;
+    progress = static_cast<int>(delta*100.0);
     auto msg = fmt::format("Starting alchemical transformation at {:>3d}%\n", progress);
     if (universe->uscreen) fmt::print(universe->uscreen, msg);
     if (universe->ulogfile) fmt::print(universe->ulogfile, msg);
@@ -175,16 +192,6 @@ void FixAlchemy::post_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-static double get_lambda(const bigint &step, const bigint &begin, const bigint &end, int iworld)
-{
-  double lambda = step - begin;
-  if (lambda != 0.0) lambda /= end - begin;
-  if (iworld == 0) lambda = 1.0 - lambda;
-  return lambda;
-}
-
-/* ---------------------------------------------------------------------- */
-
 void FixAlchemy::post_force(int /*vflag*/)
 {
   if (3 * atom->nmax > nmax) {
@@ -194,7 +201,7 @@ void FixAlchemy::post_force(int /*vflag*/)
 
   const int nall = 3 * atom->nlocal;
   double *f = &atom->f[0][0];
-  lambda = get_lambda(update->ntimestep, update->beginstep, update->endstep, universe->iworld);
+  lambda = input->variable->compute_equal(ivar);
 
   for (int i = 0; i < nall; ++i) commbuf[i] = f[i] * lambda;
   MPI_Allreduce(commbuf, f, nall, MPI_DOUBLE, MPI_SUM, samerank);
@@ -218,7 +225,10 @@ void FixAlchemy::post_force(int /*vflag*/)
   // print progress info
 
   if (universe->me == 0) {
-    int status = static_cast<int>(100.0 - lambda * 100.0);
+    double delta = update->ntimestep - update->beginstep;
+    if ((delta != 0.0) && (update->beginstep != update->endstep))
+      delta /= update->endstep - update->beginstep;
+    int status = static_cast<int>(delta*100.0);
     if ((status / 10) > (progress / 10)) {
       progress = status;
       auto msg = fmt::format("  Alchemical transformation progress: {:>3d}%\n", progress);

src/REPLICA/fix_alchemy.h

@@ -42,14 +42,15 @@ class FixAlchemy : public Fix {
   MPI_Comm samerank;
   double *commbuf;
   class Compute *pe, *temp, *press;
-  std::string id_pe, id_temp, id_press;
+  std::string id_pe, id_temp, id_press, id_lambda;
-  double lambda;         // changes from 0 to 1 during run
+  double lambda;         // scaling prefactor for combining the partitions
   double epot[3];        // last (unscaled) potential energy from each replica and combined energy
   double pressure[6];    // joined pressure
   int progress;          // for progress indicator
   int sync_box;          // 1 of box dimensions need to be synchronized
   int ilevel_respa;
   int nmax;
+  int ivar;
 };
 }    // namespace LAMMPS_NS