add documentation

doc/src/Commands_compute.rst

@@ -104,6 +104,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`pe/tally <compute_tally>`
    * :doc:`plasticity/atom <compute_plasticity_atom>`
    * :doc:`pressure <compute_pressure>`
+   * :doc:`pressure/alchemy <compute_pressure_alchemy>`
    * :doc:`pressure/uef <compute_pressure_uef>`
    * :doc:`property/atom <compute_property_atom>`
    * :doc:`property/chunk <compute_property_chunk>`

doc/src/Commands_fix.rst

@@ -29,6 +29,7 @@ OPT.
    * :doc:`adapt/fep <fix_adapt_fep>`
    * :doc:`addforce <fix_addforce>`
    * :doc:`addtorque <fix_addtorque>`
+   * :doc:`alchemy <fix_alchemy>`
    * :doc:`amoeba/bitorsion <fix_amoeba_bitorsion>`
    * :doc:`amoeba/pitorsion <fix_amoeba_pitorsion>`
    * :doc:`append/atoms <fix_append_atoms>`

doc/src/compute.rst

@@ -258,6 +258,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`pe/tally <compute_tally>` - potential energy between two groups of atoms via the tally callback mechanism
 * :doc:`plasticity/atom <compute_plasticity_atom>` - Peridynamic plasticity for each atom
 * :doc:`pressure <compute_pressure>` - total pressure and pressure tensor
+* :doc:`pressure/alchemy <compute_pressure_alchemy>` - mixed system total pressure and pressure tensor for :doc:`fix alchemy <fix_alchemy>` runs
 * :doc:`pressure/uef <compute_pressure_uef>` - pressure tensor in the reference frame of an applied flow field
 * :doc:`property/atom <compute_property_atom>` - convert atom attributes to per-atom vectors/arrays
 * :doc:`property/chunk <compute_property_chunk>` - extract various per-chunk attributes

doc/src/compute_pressure_alchemy.rst

@@ -0,0 +1,80 @@
+.. index:: compute pressure/alchemy
+
+compute pressure/alchemy command
+================================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID pressure/alchemy fix-ID
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* pressure/alchemy = style name of this compute command
+* fix-ID = ID of :doc:`fix alchemy <fix_alchemy>` command
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix trans all alchemy
+   compute mixed all pressure/alchemy trans
+   thermo_modify press mixed
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+Define a compute style that makes the "mixed" system pressure available
+for a system that uses the :doc:`fix alchemy <fix_alchemy>` command to
+transform one topology to another.  This can be used in combination with
+either :doc:`thermo_modify press <thermo_modify>` or :doc:`fix_modify
+press <fix_modify>` to output and access a pressure consistent with the
+simulated combined two topology system.
+
+The actual pressure is determined with :doc:`compute pressure
+<compute_pressure>` commands that are internally used by :doc:`fix
+alchemy <fix_alchemy>` for each topology individually and then combined.
+This command just extracts the information from the fix.
+
+The ``examples/PACKAGES/alchemy`` folder contains an example input for this command.
+
+----------
+
+Output info
+"""""""""""
+
+This compute calculates a global scalar (the pressure) and a global
+vector of length 6 (the pressure tensor), which can be accessed by
+indices 1--6.  These values can be used by any command that uses global
+scalar or vector values from a compute as input.  See the :doc:`Howto
+output <Howto_output>` page for an overview of LAMMPS output options.
+
+The ordering of values in the symmetric pressure tensor is as follows:
+:math:`p_{xx},` :math:`p_{yy},` :math:`p_{zz},` :math:`p_{xy},`
+:math:`p_{xz},` :math:`p_{yz}.`
+
+The scalar and vector values calculated by this compute are "intensive".
+The scalar and vector values will be in pressure :doc:`units <units>`.
+
+Restrictions
+""""""""""""
+
+This compute is part of the REPLICA package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+
+Related commands
+""""""""""""""""
+
+:doc:`fix alchemy <fix_alchemy>`, :doc:`compute pressure <compute_pressure>`,
+:doc:`thermo_modify <thermo_modify>`, :doc:`fix_modify <fix_modify>`
+
+Default
+"""""""
+
+none

doc/src/fix.rst

@@ -181,6 +181,7 @@ accelerated styles exist.
 * :doc:`adapt/fep <fix_adapt_fep>` - enhanced version of fix adapt
 * :doc:`addforce <fix_addforce>` - add a force to each atom
 * :doc:`addtorque <fix_addtorque>` - add a torque to a group of atoms
+* :doc:`alchemy <fix_alchemy>` - perform an "alchemical transformation" between two partitions
 * :doc:`amoeba/bitorsion <fix_amoeba_bitorsion>` - torsion/torsion terms in AMOEBA force field
 * :doc:`amoeba/pitorsion <fix_amoeba_pitorsion>` - 6-body terms in AMOEBA force field
 * :doc:`append/atoms <fix_append_atoms>` - append atoms to a running simulation

doc/src/fix_alchemy.rst

@@ -0,0 +1,104 @@
+.. index:: fix alchemy
+
+fix alchemy command
+===================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fix ID group-ID alchemy
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* alchemy = style name of this fix command
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix trans all alchemy
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+This fix command enables running an "alchemical transformation" MD
+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
+*down* from 1.0 to 0.0 for the *first* partition (*p=0*) and ramped *up*
+from 0.0 to 1.0 for the *second* partition (*p=1*).  The forces used for
+the propagation of the atoms will be the sum of the forces of the two
+systems combined and scaled with their respective :math:`\lambda_p`
+factor.  This allows to perform transformations that are not easily
+possible with :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
+synchronize them during the run.  It is thus not possible to initialize
+the two partitions by reading different data files or creating different
+systems from scratch, but rather they have to be started from the same
+system and then the desired modifications need to be applied to the
+system of the second partition.  The commands :doc:`fix adapt <fix_adapt>`
+or :doc:`fix adapt/fep <fix_adapt_fep>` could be used for simulations
+where the requirements for fix alchemy are not given.
+
+The ``examples/PACKAGES/alchemy`` folder contains example inputs for
+this command.
+
+----------
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+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
+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".  This values can be
+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 or 0.0 to 1.0 during the steps of a run, with
+*start/stop* keywords the ramp us 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
+""""""""""""
+
+This fix is part of the REPLICA package.  It is only enabled if LAMMPS
+was built with that package.  See the :doc:`Build package
+<Build_package>` page for more info.
+
+There may be only one instance of this fix in use at any time.
+
+This fix requires to perform a :ref:`multi-partition run <partition>`
+with *exactly* two partitions.
+
+This fix is *not* compatible with :doc:`load balancing <fix_balance>`.
+
+Related commands
+""""""""""""""""
+
+:doc:`compute pressure/alchemy <compute_pressure_alchemy>` command,
+:doc:`fix adapt <fix_adapt>` command
+
+Default
+"""""""
+
+none

src/REPLICA/fix_alchemy.cpp

@@ -51,12 +51,14 @@ FixAlchemy::FixAlchemy(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
   nmax = 6;
   sync_box = 0;
 
-  // set up rank-to-rank communicator
+  // set up rank-to-rank communicator for inter-partition communication
+
   int color = comm->me;
   int key = universe->iworld;
   MPI_Comm_split(universe->uworld, color, key, &samerank);
 
   // check that we have the same domain decomposition on all ranks
+
   int my_nlocal[2] = {0, 0};
   int all_nlocal[2] = {0, 0};
   my_nlocal[universe->iworld] = atom->nlocal;
@@ -124,7 +126,11 @@ void FixAlchemy::init()
   if (modify->get_fix_by_style("^balance").size() > 0)
     error->all(FLERR, "Fix alchemy is not compatible with load balancing");
 
+  if (modify->get_fix_by_style("^alchemy").size() > 1)
+    error->all(FLERR, "There may only one fix alchemy at a time");
+
   // synchronize box dimensions, determine if resync during run will be needed.
+
   synchronize_box(domain, samerank);
 
   sync_box = 0;
@@ -158,10 +164,12 @@ void FixAlchemy::setup(int vflag)
 void FixAlchemy::post_integrate()
 {
   // synchronize atom positions
+
   const int nall = atom->nlocal + atom->nghost;
   MPI_Bcast(&atom->x[0][0], 3 * nall, MPI_DOUBLE, 0, samerank);
 
   // synchronize box dimensions, if needed
+
   if (sync_box) synchronize_box(domain, samerank);
 }
 
@@ -192,6 +200,7 @@ void FixAlchemy::post_force(int /*vflag*/)
   MPI_Allreduce(commbuf, f, nall, MPI_DOUBLE, MPI_SUM, samerank);
 
   // sum up potential energy
+
   const double scalefac = 1.0 / comm->nprocs;
   commbuf[0] = commbuf[1] = commbuf[2] = 0.0;
   commbuf[universe->iworld] = scalefac * pe->compute_scalar();
@@ -200,12 +209,14 @@ void FixAlchemy::post_force(int /*vflag*/)
   pe->addstep(update->ntimestep + 1);
 
   // sum up pressure
+
   press->compute_vector();
   for (int i = 0; i < 6; ++i) commbuf[i] = lambda * scalefac * press->vector[i];
   MPI_Allreduce(commbuf, pressure, 6, MPI_DOUBLE, MPI_SUM, universe->uworld);
   press->addstep(update->ntimestep + 1);
 
   // print progress info
+
   if (universe->me == 0) {
     int status = static_cast<int>(100.0 - lambda * 100.0);
     if ((status / 10) > (progress / 10)) {