diff --git a/examples/PACKAGES/alchemy/in.twowater b/examples/PACKAGES/alchemy/in.twowater index dd246666bf..63089d23c5 100644 --- a/examples/PACKAGES/alchemy/in.twowater +++ b/examples/PACKAGES/alchemy/in.twowater @@ -1,8 +1,8 @@ # Example for alchemical transformation of two water molecules into a hydronium and hydroxyl ion -# WARNING: This input is for demonstrating the method. -# Force field parameters are made up and not suitable for production simulations. +# WARNING: This input is intended for demonstrating the method, only +# The force field parameters are made up and NOT suitable for production simulations. -# set up names for two partitions with different topologies for output files. +# set up different names for two partitions variable name world twowater twoions units real @@ -36,19 +36,19 @@ create_atoms 0 single 2.0 0.0 0.0 mol water 767353 # ... and put them in a group group transform id 1:6 -# now fill the box with more water -create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.4 +# now fill the rest of the box with more water +create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.33 # change topology and settings for the two states # we cannot simply create a different topology directly or # load a different data file because the order and position # of all atoms must be maintained across both replica -# for that we first have to remove everything. +# we first have to remove all topology data in the transform group delete_bonds transform bond 1 delete_bonds transform angle 1 remove -# generate different topologies with the same number of atoms +# then generate different topologies for the two partitions. select by name. if "${name} == twowater" then & "create_bonds single/bond 2 1 2" & "create_bonds single/bond 2 1 3" & @@ -73,14 +73,22 @@ else & velocity all create 300.0 5463576 timestep 0.2 -fix 1 all nvt temp 300 300 1.0 -fix 2 all alchemy -fix 3 all shake 0.0001 50 1000 b 1 a 1 +# since the trajectory and forces are kept identical through fix alchemy +# we can even do fix npt simulations, but we need to use the "mixed" pressure +fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0 +fix transform all alchemy +compute pressure all pressure/alchemy transform +fix_modify integrate press pressure + +# only need to output a dump file form one partition if "${name} == twowater" then & "dump 1 all atom 100 ${name}.lammpstrj" & "dump_modify 1 sort id" -thermo_style custom step temp press etotal density pe ke f_2[*] +thermo_style custom step temp press etotal density pe ke f_transform[1] f_transform[4] +thermo_modify colname f_transform[1] lambda colname f_transform[4] EPot_mixed +thermo_modify press pressure + thermo 100 run 20000 diff --git a/src/.gitignore b/src/.gitignore index 9da19ca982..0d73a7911d 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1465,6 +1465,8 @@ /zstd_file_writer.cpp /zstd_file_writer.h +/compute_pressure_alchemy.cpp +/compute_pressure_alchemy.h /atom_vec_smd.cpp /atom_vec_smd.h /compute_saed.cpp diff --git a/src/REPLICA/compute_pressure_alchemy.cpp b/src/REPLICA/compute_pressure_alchemy.cpp new file mode 100644 index 0000000000..d19ea1b893 --- /dev/null +++ b/src/REPLICA/compute_pressure_alchemy.cpp @@ -0,0 +1,102 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "compute_pressure_alchemy.h" + +#include "domain.h" +#include "error.h" +#include "fix.h" +#include "modify.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputePressureAlchemy::ComputePressureAlchemy(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR, "Illegal compute pressure/alchemy command"); + if (igroup) error->all(FLERR, "Compute pressure/alchemy must use group all"); + + scalar_flag = vector_flag = 1; + size_vector = 6; + extscalar = 0; + extvector = 0; + pressflag = 1; + timeflag = 1; + + id_fix = arg[3]; + if (!modify->get_fix_by_id(id_fix)) + error->all(FLERR, "Could not find compute pressure/alchemy fix ID {} for fix alchemy", id_fix); + + vector = new double[size_vector]; +} + +/* ---------------------------------------------------------------------- */ + +ComputePressureAlchemy::~ComputePressureAlchemy() +{ + delete[] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputePressureAlchemy::init() +{ + + fix = modify->get_fix_by_id(id_fix); + if (!fix) + error->all(FLERR, "Could not find compute pressure/alchemy fix ID {} for fix alchemy", id_fix); + + int dim = 0; + void *ptr = fix->extract("pressure", dim); + if (!ptr || (dim != 1)) error->all(FLERR, "Could not extract pressure from fix alchemy"); +} + +/* ---------------------------------------------------------------------- + compute total pressure from tensor, averaged over Pxx, Pyy, Pzz +------------------------------------------------------------------------- */ + +double ComputePressureAlchemy::compute_scalar() +{ + invoked_scalar = update->ntimestep; + if (update->vflag_global != invoked_scalar) + error->all(FLERR, "Virial was not tallied on needed timestep"); + + compute_vector(); + + if (domain->dimension == 3) { + scalar = (vector[0] + vector[1] + vector[2]) / 3.0; + } else { + scalar = (vector[0] + vector[1]) / 2.0; + } + return scalar; +} + +/* ---------------------------------------------------------------------- + extract compute combined system pressure tensor from alchemy fix +------------------------------------------------------------------------- */ + +void ComputePressureAlchemy::compute_vector() +{ + invoked_vector = update->ntimestep; + if (update->vflag_global != invoked_vector) + error->all(FLERR, "Virial was not tallied on needed timestep"); + + int dim = 0; + double *pressure = (double *) fix->extract("pressure", dim); + if (!pressure || (dim != 1)) error->all(FLERR, "Could not extract pressure from fix alchemy"); + + for (int i = 0; i < 6; i++) vector[i] = pressure[i]; +} diff --git a/src/REPLICA/compute_pressure_alchemy.h b/src/REPLICA/compute_pressure_alchemy.h new file mode 100644 index 0000000000..0ce0039587 --- /dev/null +++ b/src/REPLICA/compute_pressure_alchemy.h @@ -0,0 +1,41 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(pressure/alchemy,ComputePressureAlchemy); +// clang-format on +#else + +#ifndef LMP_COMPUTE_PRESSURE_ALCHEMY_H +#define LMP_COMPUTE_PRESSURE_ALCHEMY_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputePressureAlchemy : public Compute { + public: + ComputePressureAlchemy(class LAMMPS *, int, char **); + ~ComputePressureAlchemy() override; + void init() override; + double compute_scalar() override; + void compute_vector() override; + + protected: + class Fix *fix; + std::string id_fix; +}; +} // namespace LAMMPS_NS +#endif +#endif diff --git a/src/REPLICA/fix_alchemy.cpp b/src/REPLICA/fix_alchemy.cpp index 104212cb9b..75b4c4cb73 100644 --- a/src/REPLICA/fix_alchemy.cpp +++ b/src/REPLICA/fix_alchemy.cpp @@ -23,6 +23,8 @@ #include "universe.h" #include "update.h" +#include + using namespace LAMMPS_NS; using namespace FixConst; @@ -38,26 +40,46 @@ static double get_lambda(const bigint &step, const bigint &begin, const bigint & /* ---------------------------------------------------------------------- */ -FixAlchemy::FixAlchemy(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +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"); lambda = epot[0] = epot[1] = 0.0; + for (int i = 0; i < 6; ++i) pressure[i] = 0.0; no_change_box = 1; vector_flag = 1; - size_vector = 4; + size_vector = 10; extvector = 0; ilevel_respa = 0; + nmax = 6; - // set up communicators + // set up rank-to-rank communicator if (universe->nworlds != 2) error->all(FLERR, "Must use exactly two partitions"); 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; + MPI_Allreduce(my_nlocal, all_nlocal, 2, MPI_INT, MPI_SUM, samerank); + int fail = (all_nlocal[0] == all_nlocal[1]) ? 0 : 1; + int allfail = 0; + MPI_Allreduce(&fail, &allfail, 1, MPI_INT, MPI_MAX, universe->uworld); + if (allfail) + error->all(FLERR, "Number of atoms and domain decomposition must match for both partitions"); + id_pe = std::string(id) + "_pe"; - modify->add_compute(id_pe + " all pe"); + pe = modify->add_compute(id_pe + " all pe"); + pe->addstep(update->ntimestep); + id_temp = std::string(id) + "_temp"; + temp = modify->add_compute(id_temp + " all temp"); + temp->addstep(update->ntimestep); + id_press = std::string(id) + "_press"; + press = modify->add_compute(id_press + " all pressure " + id_temp); + press->addstep(update->ntimestep); } /* ---------------------------------------------------------------------- */ @@ -66,6 +88,9 @@ FixAlchemy::~FixAlchemy() { MPI_Comm_free(&samerank); modify->delete_compute(id_pe); + modify->delete_compute(id_temp); + modify->delete_compute(id_press); + memory->destroy(commbuf); } /* ---------------------------------------------------------------------- */ @@ -80,7 +105,13 @@ int FixAlchemy::setmask() /* ---------------------------------------------------------------------- */ -void FixAlchemy::init() {} +void FixAlchemy::init() +{ + int onenmax = MAX(nmax, 3 * atom->nmax); + MPI_Allreduce(&onenmax, &nmax, 1, MPI_INT, MPI_MAX, universe->uworld); + memory->destroy(commbuf); + memory->create(commbuf, sizeof(double) * nmax, "alchemy:nmax"); +} /* ---------------------------------------------------------------------- */ @@ -108,26 +139,31 @@ void FixAlchemy::post_integrate() void FixAlchemy::post_force(int /*vflag*/) { - const int nall = atom->nlocal; - auto f = atom->f; + if (3 * atom->nmax > nmax) { + nmax = 3 * atom->nmax; + memory->grow(commbuf, sizeof(double) * atom->nmax, "alchemy:commbuf"); + } + + const int nall = 3 * atom->nlocal; + double *f = &atom->f[0][0]; lambda = get_lambda(update->ntimestep, update->beginstep, update->endstep, universe->iworld); - double *commbuf = new double[3 * nall]; - int m = 0; - for (int i = 0; i < nall; ++i) { - commbuf[m++] = f[i][0] * lambda; - commbuf[m++] = f[i][1] * lambda; - commbuf[m++] = f[i][2] * lambda; - } - MPI_Allreduce(commbuf, &f[0][0], 3 * nall, MPI_DOUBLE, MPI_SUM, samerank); + for (int i = 0; i < nall; ++i) commbuf[i] = f[i] * lambda; + MPI_Allreduce(commbuf, f, nall, MPI_DOUBLE, MPI_SUM, samerank); - auto pe = modify->get_compute_by_id(id_pe); + // sum up potential energy + const double scalefac = 1.0 / comm->nprocs; commbuf[0] = commbuf[1] = commbuf[2] = 0.0; - commbuf[universe->iworld] = pe->compute_scalar() / comm->nprocs; - commbuf[2] = lambda * pe->compute_scalar() / comm->nprocs; + commbuf[universe->iworld] = scalefac * pe->compute_scalar(); + commbuf[2] = lambda * scalefac * pe->compute_scalar(); MPI_Allreduce(commbuf, epot, 3, MPI_DOUBLE, MPI_SUM, universe->uworld); - delete[] commbuf; 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); } /* ---------------------------------------------------------------------- */ @@ -135,5 +171,18 @@ void FixAlchemy::post_force(int /*vflag*/) double FixAlchemy::compute_vector(int n) { if (n == 0) return lambda; - return epot[n - 1]; + if ((n > 0) && (n < 4)) return epot[n - 1]; + return pressure[n - 4]; +} + +/* ---------------------------------------------------------------------- */ + +void *FixAlchemy::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str, "lambda") == 0) { return λ } + if (strcmp(str, "pe") == 0) { return &epot[2]; } + dim = 1; + if (strcmp(str, "pressure") == 0) { return pressure; } + return nullptr; } diff --git a/src/REPLICA/fix_alchemy.h b/src/REPLICA/fix_alchemy.h index 0457d12735..9794629db3 100644 --- a/src/REPLICA/fix_alchemy.h +++ b/src/REPLICA/fix_alchemy.h @@ -35,13 +35,18 @@ class FixAlchemy : public Fix { void post_integrate() override; void post_force(int) override; double compute_vector(int) override; + void *extract(const char *, int &) override; protected: MPI_Comm samerank; - std::string id_pe; - double lambda; // changes from 0 to 1 during run - double epot[3]; // last (unscaled) potential energy from each replica and combined energy + double *commbuf; + class Compute *pe, *temp, *press; + std::string id_pe, id_temp, id_press; + double lambda; // changes from 0 to 1 during run + double epot[3]; // last (unscaled) potential energy from each replica and combined energy + double pressure[6]; // joined pressure int ilevel_respa; + int nmax; }; } // namespace LAMMPS_NS