diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index aaa899ef66..c2753b2afc 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -104,6 +104,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`pe/tally ` * :doc:`plasticity/atom ` * :doc:`pressure ` + * :doc:`pressure/alchemy ` * :doc:`pressure/uef ` * :doc:`property/atom ` * :doc:`property/chunk ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index d275b33eba..5b331af1aa 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -29,6 +29,7 @@ OPT. * :doc:`adapt/fep ` * :doc:`addforce ` * :doc:`addtorque ` + * :doc:`alchemy ` * :doc:`amoeba/bitorsion ` * :doc:`amoeba/pitorsion ` * :doc:`append/atoms ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 72ea31bbb4..3d74b3884e 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -258,6 +258,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`pe/tally ` - potential energy between two groups of atoms via the tally callback mechanism * :doc:`plasticity/atom ` - Peridynamic plasticity for each atom * :doc:`pressure ` - total pressure and pressure tensor +* :doc:`pressure/alchemy ` - mixed system total pressure and pressure tensor for :doc:`fix alchemy ` runs * :doc:`pressure/uef ` - pressure tensor in the reference frame of an applied flow field * :doc:`property/atom ` - convert atom attributes to per-atom vectors/arrays * :doc:`property/chunk ` - extract various per-chunk attributes diff --git a/doc/src/compute_pressure_alchemy.rst b/doc/src/compute_pressure_alchemy.rst new file mode 100644 index 0000000000..bdf9802e20 --- /dev/null +++ b/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 ` command +* pressure/alchemy = style name of this compute command +* fix-ID = ID of :doc:`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 ` command to +transform one topology to another. This can be used in combination with +either :doc:`thermo_modify press ` or :doc:`fix_modify +press ` to output and access a pressure consistent with the +simulated combined two topology system. + +The actual pressure is determined with :doc:`compute pressure +` commands that are internally used by :doc:`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 ` 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 `. + +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 +` page for more info. + + +Related commands +"""""""""""""""" + +:doc:`fix alchemy `, :doc:`compute pressure `, +:doc:`thermo_modify `, :doc:`fix_modify ` + +Default +""""""" + +none diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 2dfe97a3ec..71f8bf9e16 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -181,6 +181,7 @@ accelerated styles exist. * :doc:`adapt/fep ` - enhanced version of fix adapt * :doc:`addforce ` - add a force to each atom * :doc:`addtorque ` - add a torque to a group of atoms +* :doc:`alchemy ` - perform an "alchemical transformation" between two partitions * :doc:`amoeba/bitorsion ` - torsion/torsion terms in AMOEBA force field * :doc:`amoeba/pitorsion ` - 6-body terms in AMOEBA force field * :doc:`append/atoms ` - append atoms to a running simulation diff --git a/doc/src/fix_alchemy.rst b/doc/src/fix_alchemy.rst new file mode 100644 index 0000000000..1ee9acc629 --- /dev/null +++ b/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 ` 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 ` 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 ` or :doc:`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 ` +or :doc:`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 `. +None of the :doc:`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 ` 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 ` 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 ` 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 `. + +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 +` 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 ` +with *exactly* two partitions. + +This fix is *not* compatible with :doc:`load balancing `. + +Related commands +"""""""""""""""" + +:doc:`compute pressure/alchemy ` command, +:doc:`fix adapt ` command + +Default +""""""" + +none diff --git a/src/REPLICA/fix_alchemy.cpp b/src/REPLICA/fix_alchemy.cpp index 6a65c0d116..18e55603b9 100644 --- a/src/REPLICA/fix_alchemy.cpp +++ b/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(100.0 - lambda * 100.0); if ((status / 10) > (progress / 10)) {