diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index cfc896aa0e..aaf706b5df 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -42,6 +42,7 @@ OPT. * :doc:`gaussian ` * :doc:`gromos (o) ` * :doc:`harmonic (iko) ` + * :doc:`harmonic/restrain ` * :doc:`harmonic/shift (o) ` * :doc:`harmonic/shift/cut (o) ` * :doc:`lepton (o) ` diff --git a/doc/src/bond_harmonic_restrain.rst b/doc/src/bond_harmonic_restrain.rst new file mode 100644 index 0000000000..a9918d8b5d --- /dev/null +++ b/doc/src/bond_harmonic_restrain.rst @@ -0,0 +1,79 @@ +.. index:: bond_style harmonic/restrain + +bond_style harmonic/restrain command +==================================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + bond_style harmonic/restrain + +Examples +"""""""" + +.. code-block:: LAMMPS + + bond_style harmonic + bond_coeff 5 80.0 + +Description +""""""""""" + +The *harmonic/restrain* bond style uses the potential + +.. math:: + + E = K (r - r_{t=0})^2 + +where :math:`r_{t=0}` is the distance between the bond atoms at the +beginning of the first :doc:`run ` or :doc:`minimize ` +command after the bond style has been defined (*t=0*). Note that the +usual 1/2 factor is included in :math:`K`. This will effectively +restrain bonds to their initial length, whatever that is. + +The following coefficient must be defined for each bond type via the +:doc:`bond_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`read_restart ` commands + +* :math:`K` (energy/distance\^2) + +Restart info +"""""""""""" + +This bond style supports the :doc:`write_restart ` and +:doc:`read_restart ` commands. The state of the initial +bond lengths is stored with the restart files and read back. + +Restrictions +"""""""""""" + +This bond style can only be used if LAMMPS was built with the +EXTRA-MOLECULE package. See the :doc:`Build package ` +page for more info. + +This bond style does not write its coefficients to a data file since the +:math:`r_{t=0}` values are for individual bonds and not bond types. +Since it uses :doc:`fix property/atom ` internally, +you should use the "nofix" argument to the :doc:`write_data command +` to avoid writing a section to the data file can cannot be +read back in by the bond style. + +This bond style cannot be used with :doc:`fix shake or fix rattle +`, with :doc:`fix filter/corotate `, or +any :doc:`tip4p pair style ` since there is no specific +equilibrium distance for a given bond type. + + +Related commands +"""""""""""""""" + +:doc:`bond_coeff `, :doc:`delete_bonds `, +:doc:`bond_harmonic ` + +Default +""""""" + +none diff --git a/doc/src/bond_style.rst b/doc/src/bond_style.rst index 23b89d00a2..b33d0a9e9a 100644 --- a/doc/src/bond_style.rst +++ b/doc/src/bond_style.rst @@ -10,7 +10,7 @@ Syntax bond_style style args -* style = *none* or *zero* or *hybrid* or *bpm/rotational* or *bpm/spring* or *class2* or *fene* or *fene/expand* or *fene/nm* or *gaussian* or *gromos* or *harmonic* or *harmonic/shift* or *harmonic/shift/cut* or *lepton* or *morse* or *nonlinear* or *oxdna/fene* or *oxdena2/fene* or *oxrna2/fene* or *quartic* or *special* or *table* +* style = *none* or *zero* or *hybrid* or *bpm/rotational* or *bpm/spring* or *class2* or *fene* or *fene/expand* or *fene/nm* or *gaussian* or *gromos* or *harmonic* or *harmonic/restrain* *harmonic/shift* or *harmonic/shift/cut* or *lepton* or *morse* or *nonlinear* or *oxdna/fene* or *oxdena2/fene* or *oxrna2/fene* or *quartic* or *special* or *table* * args = none for any style except *hybrid* @@ -93,6 +93,7 @@ accelerated styles exist. * :doc:`gaussian ` - multicentered Gaussian-based bond potential * :doc:`gromos ` - GROMOS force field bond * :doc:`harmonic ` - harmonic bond +* :doc:`harmonic/restrain ` - harmonic bond to restrain to original bond distance * :doc:`harmonic/shift ` - shifted harmonic bond * :doc:`harmonic/shift/cut ` - shifted harmonic bond with a cutoff * :doc:`lepton ` - bond potential from evaluating a string diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index d03300439f..ff8f0d455b 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2463,6 +2463,7 @@ nodeless nodeset nodesets Noehring +nofix Noffset noforce noguess diff --git a/src/.gitignore b/src/.gitignore index df02ab3397..aa57927669 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -492,6 +492,8 @@ /bond_gromos.h /bond_harmonic.cpp /bond_harmonic.h +/bond_harmonic_restrain.cpp +/bond_harmonic_restrain.h /bond_harmonic_shift.cpp /bond_harmonic_shift.h /bond_harmonic_shift_cut.cpp diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp new file mode 100644 index 0000000000..413d00226d --- /dev/null +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.cpp @@ -0,0 +1,244 @@ +/* ---------------------------------------------------------------------- + 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 "bond_harmonic_restrain.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "fix_property_atom.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" + +#include +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondHarmonicRestrain::BondHarmonicRestrain(LAMMPS *_lmp) : Bond(_lmp), initial(nullptr) +{ + writedata = 0; + natoms = -1; +} + +/* ---------------------------------------------------------------------- */ + +BondHarmonicRestrain::~BondHarmonicRestrain() +{ + modify->delete_fix("BOND_RESTRAIN"); + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondHarmonicRestrain::compute(int eflag, int vflag) +{ + int i1, i2, n, type; + double delx, dely, delz, ebond, fbond; + double rsq, r, r0, dr, rk; + + ebond = 0.0; + ev_init(eflag, vflag); + + double **x = atom->x; + double **x0 = (double **) atom->extract("d2_x0"); + double **f = atom->f; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nbondlist; n++) { + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + + delx = x0[i1][0] - x0[i2][0]; + dely = x0[i1][1] - x0[i2][1]; + delz = x0[i1][2] - x0[i2][2]; + rsq = delx * delx + dely * dely + delz * delz; + r0 = sqrt(rsq); + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + + rsq = delx * delx + dely * dely + delz * delz; + r = sqrt(rsq); + dr = r - r0; + rk = k[type] * dr; + + // force & energy + + if (r > 0.0) + fbond = -2.0 * rk / r; + else + fbond = 0.0; + + if (eflag) ebond = rk * dr; + + // apply force to each of 2 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += delx * fbond; + f[i1][1] += dely * fbond; + f[i1][2] += delz * fbond; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= delx * fbond; + f[i2][1] -= dely * fbond; + f[i2][2] -= delz * fbond; + } + + if (evflag) ev_tally(i1, i2, nlocal, newton_bond, ebond, fbond, delx, dely, delz); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondHarmonicRestrain::allocate() +{ + allocated = 1; + const int np1 = atom->nbondtypes + 1; + + memory->create(k, np1, "bond:k"); + + memory->create(setflag, np1, "bond:setflag"); + for (int i = 1; i < np1; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::coeff(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR, "Incorrect args for bond coefficients"); + if (!allocated) allocate(); + + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error); + + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + initialize custom data storage +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::init_style() +{ + // create internal fix to store initial positions + initial = dynamic_cast( + modify->add_fix("BOND_RESTRAIN all property/atom d2_x0 3 ghost yes")); + if (!initial) error->all(FLERR, "Failure to create internal per-atom storage"); + + // store initial positions + if (natoms < 0) { + natoms = atom->natoms; + double *const *const x0 = (double **) atom->extract("d2_x0"); + const double *const *const x = atom->x; + for (int i = 0; i < atom->nlocal; ++i) + for (int j = 0; j < 3; ++j) x0[i][j] = x[i][j]; + } + + // must not add atoms + if (natoms < atom->natoms) + error->all(FLERR, "Bond style harmonic/restrain does not support adding atoms"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::write_restart(FILE *fp) +{ + fwrite(&natoms, sizeof(bigint), 1, fp); + fwrite(&k[1], sizeof(double), atom->nbondtypes, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR, &natoms, sizeof(bigint), 1, fp, nullptr, error); + utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + } + MPI_Bcast(&natoms, 1, MPI_LMP_BIGINT, 0, world); + MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondHarmonicRestrain::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) fprintf(fp, "%d %g\n", i, k[i]); +} + +/* ---------------------------------------------------------------------- */ + +double BondHarmonicRestrain::single(int type, double rsq, int i, int j, double &fforce) +{ + double **x0 = (double **) atom->extract("d2_x0"); + + double delx = x0[i][0] - x0[j][0]; + double dely = x0[i][1] - x0[j][1]; + double delz = x0[i][2] - x0[j][2]; + double r0 = sqrt(delx * delx + dely * dely + delz * delz); + + double r = sqrt(rsq); + double dr = r - r0; + double rk = k[type] * dr; + fforce = 0; + if (r > 0.0) fforce = -2.0 * rk / r; + return rk * dr; +} + +/* ---------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + return ptr to internal members upon request +------------------------------------------------------------------------ */ + +void *BondHarmonicRestrain::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "k") == 0) return (void *) k; + return nullptr; +} diff --git a/src/EXTRA-MOLECULE/bond_harmonic_restrain.h b/src/EXTRA-MOLECULE/bond_harmonic_restrain.h new file mode 100644 index 0000000000..e17c21e125 --- /dev/null +++ b/src/EXTRA-MOLECULE/bond_harmonic_restrain.h @@ -0,0 +1,52 @@ +/* -*- 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 BOND_CLASS +// clang-format off +BondStyle(harmonic/restrain,BondHarmonicRestrain); +// clang-format on +#else + +#ifndef LMP_BOND_HARMONIC_RESTRAIN_H +#define LMP_BOND_HARMONIC_RESTRAIN_H + +#include "bond.h" + +namespace LAMMPS_NS { + +class BondHarmonicRestrain : public Bond { + public: + BondHarmonicRestrain(class LAMMPS *); + ~BondHarmonicRestrain() override; + void compute(int, int) override; + void coeff(int, char **) override; + void init_style() override; + double equilibrium_distance(int) override { return -1.0; }; // return invalid value since not applicable + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_data(FILE *) override; + double single(int, double, int, int, double &) override; + void *extract(const char *, int &) override; + + protected: + double *k; + bigint natoms; + class FixPropertyAtom *initial; + + virtual void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif