From ed822f8002d0438d751f668f1b7707dbb572750a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 24 Feb 2023 15:41:04 -0500 Subject: [PATCH] WIP --- examples/PACKAGES/alchemy/h2o.mol | 62 ++++++++++++ examples/PACKAGES/alchemy/h3o+.mol | 39 ++++++++ examples/PACKAGES/alchemy/in.twowater | 53 ++++++++++ examples/PACKAGES/alchemy/oh-.mol | 24 +++++ src/REPLICA/fix_alchemy.cpp | 134 ++++++++++++++++++++++++++ src/REPLICA/fix_alchemy.h | 49 ++++++++++ 6 files changed, 361 insertions(+) create mode 100644 examples/PACKAGES/alchemy/h2o.mol create mode 100644 examples/PACKAGES/alchemy/h3o+.mol create mode 100644 examples/PACKAGES/alchemy/in.twowater create mode 100644 examples/PACKAGES/alchemy/oh-.mol create mode 100644 src/REPLICA/fix_alchemy.cpp create mode 100644 src/REPLICA/fix_alchemy.h diff --git a/examples/PACKAGES/alchemy/h2o.mol b/examples/PACKAGES/alchemy/h2o.mol new file mode 100644 index 0000000000..e5a5e4fe93 --- /dev/null +++ b/examples/PACKAGES/alchemy/h2o.mol @@ -0,0 +1,62 @@ +# Water molecule. SPC/E model. + +3 atoms +2 bonds +1 angles + +Coords + +1 1.12456 0.09298 1.27452 +2 1.53683 0.75606 1.89928 +3 0.49482 0.56390 0.65678 + +Types + +1 1 +2 2 +3 2 + +Charges + +1 -0.8472 +2 0.4236 +3 0.4236 + +Bonds + +1 1 1 2 +2 1 1 3 + +Angles + +1 1 2 1 3 + +Shake Flags + +1 1 +2 1 +3 1 + +Shake Atoms + +1 1 2 3 +2 1 2 3 +3 1 2 3 + +Shake Bond Types + +1 1 1 1 +2 1 1 1 +3 1 1 1 + +Special Bond Counts + +1 2 0 0 +2 1 1 0 +3 1 1 0 + +Special Bonds + +1 2 3 +2 1 3 +3 1 2 diff --git a/examples/PACKAGES/alchemy/h3o+.mol b/examples/PACKAGES/alchemy/h3o+.mol new file mode 100644 index 0000000000..699438e7a4 --- /dev/null +++ b/examples/PACKAGES/alchemy/h3o+.mol @@ -0,0 +1,39 @@ +# Hydronium ion. SPC/E model. + +4 atoms +3 bonds +3 angles + +Coords + +1 0.00000 1.00000 0.00000 +2 0.00000 0.00000 0.00000 +3 0.00000 -0.50000 -0.8660 +4 0.00000 -0.50000 0.8660 + +Types + +1 2 +2 1 +3 2 +4 2 + +Charges + +1 0.56775 +2 -0.70305 +3 0.56775 +4 0.56775 + +Bonds + +1 1 2 1 +2 1 2 3 +3 1 2 4 + +Angles + +1 1 1 2 3 +1 1 3 2 4 +1 1 1 2 4 + diff --git a/examples/PACKAGES/alchemy/in.twowater b/examples/PACKAGES/alchemy/in.twowater new file mode 100644 index 0000000000..62e0dacd60 --- /dev/null +++ b/examples/PACKAGES/alchemy/in.twowater @@ -0,0 +1,53 @@ +# set up names for two partitions with different topologies +variable name world twowater twoions + +units real +atom_style full +atom_modify map array +region box block -10 10 -10 10 -10 10 +boundary m m m +create_box 2 box bond/types 1 angle/types 1 & + extra/bond/per/atom 3 extra/angle/per/atom 3 extra/special/per/atom 3 + +mass 1 15.9994 +mass 2 1.008 + +pair_style lj/cut/coul/cut 10.0 +pair_coeff 1 1 0.1553 3.166 +pair_coeff 1 2 0.0 1.0 +pair_coeff 2 2 0.0 1.0 + +bond_style harmonic +bond_coeff 1 1000.0 1.0 + +angle_style harmonic +angle_coeff 1 100.0 109.47 + +molecule water h2o.mol +molecule hydronium h3o+.mol +molecule hydroxyl oh-.mol + + +timestep 0.1 + +# generate different topologies with the same number of atoms +if "${name} == twowater" then & + "create_atoms 0 single -2.0 0.0 0.0 mol water 453624" & + "create_atoms 0 single 2.0 0.0 0.0 mol water 767353" & +else & + "create_atoms 0 single 2.0 0.0 0.0 mol hydroxyl 767353" & + "create_atoms 0 single -2.0 0.0 0.0 mol hydronium 453624" + +velocity all create 100 5463576 + +fix 1 all nve +fix 2 all alchemy + +dump 1 all atom 100 ${name}.lammpstrj +dump_modify 1 sort id + +thermo_style custom step temp press etotal pe ke f_2[*] +thermo 100 + +reset_timestep 100 +run 10000 start 0 stop 20000 diff --git a/examples/PACKAGES/alchemy/oh-.mol b/examples/PACKAGES/alchemy/oh-.mol new file mode 100644 index 0000000000..782d1fa88f --- /dev/null +++ b/examples/PACKAGES/alchemy/oh-.mol @@ -0,0 +1,24 @@ +# Hydroxyl ion. SPC/E model. + +2 atoms +1 bonds + +Coords + +1 0.5 0.0 0.0 +2 -0.5 0.0 0.0 + +Types + +1 1 +2 2 + +Charges + +1 -1.1354 +2 0.1354 + +Bonds + +1 1 1 2 + diff --git a/src/REPLICA/fix_alchemy.cpp b/src/REPLICA/fix_alchemy.cpp new file mode 100644 index 0000000000..da10f8261f --- /dev/null +++ b/src/REPLICA/fix_alchemy.cpp @@ -0,0 +1,134 @@ +/* ---------------------------------------------------------------------- + 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 "fix_alchemy.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "memory.h" +#include "universe.h" +#include "update.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +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; +} + +/* ---------------------------------------------------------------------- */ + +FixAlchemy::FixAlchemy(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), coordbuf(nullptr) +{ + lambda = epot[0] = epot[1] = 0.0; + + no_change_box = 1; + vector_flag = 1; + size_vector = 3; + extvector = 0; + + // set up communicators + 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); + + if (narg != 3) error->all(FLERR, "Incorrect number of arguments for fix alchemy"); +} + +/* ---------------------------------------------------------------------- */ + +FixAlchemy::~FixAlchemy() +{ + MPI_Comm_free(&samerank); +} + +/* ---------------------------------------------------------------------- */ + +int FixAlchemy::setmask() +{ + int mask = 0; + mask |= POST_INTEGRATE; + mask |= POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::init() {} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::setup(int /*vflag*/) +{ + const int nall = atom->nlocal; + const auto tag = atom->tag; + const auto x = atom->x; + + memory->destroy(coordbuf); + memory->create(coordbuf, 4 * sizeof(double) * nall, "alchemy:coordbuf"); + + if (universe->iworld == 0) { + int m = 0; + for (int i = 0; i < nall; ++i) { + coordbuf[m++] = ubuf(tag[i]).d; + coordbuf[m++] = x[i][0]; + coordbuf[m++] = x[i][1]; + coordbuf[m++] = x[i][2]; + } + } + + MPI_Bcast(coordbuf, 4 * nall, MPI_DOUBLE, 0, samerank); + + if (universe->iworld == 1) { + int m = 0; + for (int i = 0; i < nall; ++i) { + tagint mytag = (tagint) ubuf(coordbuf[m++]).i; + x[atom->map(mytag)][0] = coordbuf[m++]; + x[atom->map(mytag)][1] = coordbuf[m++]; + x[atom->map(mytag)][2] = coordbuf[m++]; + } + } + + lambda = get_lambda(update->ntimestep, update->beginstep, update->endstep, universe->iworld); +} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::post_integrate() +{ + int nall = 3 * (atom->nlocal + atom->nghost); + // MPI_Bcast(&atom->x[0][0], nall, MPI_DOUBLE, 0, samerank); +} + +/* ---------------------------------------------------------------------- */ + +void FixAlchemy::post_force(int /*vflag*/) +{ + lambda = get_lambda(update->ntimestep, update->beginstep, update->endstep, universe->iworld); +} + +/* ---------------------------------------------------------------------- */ + +double FixAlchemy::compute_vector(int n) +{ + if (n == 0) return lambda; + return epot[n - 1]; +} diff --git a/src/REPLICA/fix_alchemy.h b/src/REPLICA/fix_alchemy.h new file mode 100644 index 0000000000..2c48b5b4e7 --- /dev/null +++ b/src/REPLICA/fix_alchemy.h @@ -0,0 +1,49 @@ +/* -*- 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 FIX_CLASS +// clang-format off +FixStyle(alchemy,FixAlchemy); +// clang-format on +#else + +#ifndef LMP_FIX_ALCHEMY_H +#define LMP_FIX_ALCHEMY_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAlchemy : public Fix { + public: + FixAlchemy(class LAMMPS *, int, char **); + ~FixAlchemy() override; + + int setmask() override; + void init() override; + void setup(int) override; + void post_integrate() override; + void post_force(int) override; + double compute_vector(int) override; + + protected: + MPI_Comm samerank; + MPI_Comm rankzero; + double *coordbuf; + double lambda; // changes from 0 to 1 during run + double epot[2]; // last (unscaled) potential energy from each replica +}; +} // namespace LAMMPS_NS + +#endif +#endif