From acd0a62de72bbb6d5431be3b65e8e5b8ad4f36e8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Aug 2020 14:44:18 -0400 Subject: [PATCH] add new fix style momentum/chunk --- doc/src/Commands_fix.rst | 1 + doc/src/bond_style.rst | 1 + doc/src/fix.rst | 1 + doc/src/fix_momentum.rst | 51 +++-- src/USER-MISC/fix_momentum_chunk.cpp | 287 +++++++++++++++++++++++++++ src/USER-MISC/fix_momentum_chunk.h | 63 ++++++ 6 files changed, 388 insertions(+), 16 deletions(-) create mode 100644 src/USER-MISC/fix_momentum_chunk.cpp create mode 100644 src/USER-MISC/fix_momentum_chunk.h diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 863d09e4ac..4b5894cd2e 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -98,6 +98,7 @@ OPT. * :doc:`manifoldforce ` * :doc:`meso/move ` * :doc:`momentum (k) ` + * :doc:`momentum/chunk ` * :doc:`move ` * :doc:`mscg ` * :doc:`msst ` diff --git a/doc/src/bond_style.rst b/doc/src/bond_style.rst index 52cc761da1..64fdeb95c9 100644 --- a/doc/src/bond_style.rst +++ b/doc/src/bond_style.rst @@ -98,6 +98,7 @@ accelerated styles exist. * :doc:`oxdna2/fene ` - same as oxdna but used with different pair styles * :doc:`oxrna2/fene ` - modified FENE bond suitable for RNA modeling * :doc:`quartic ` - breakable quartic bond +* :doc:`special ` - enable special bond exclusions for 1-5 pairs and beyond * :doc:`table ` - tabulated by bond length ---------- diff --git a/doc/src/fix.rst b/doc/src/fix.rst index d0fe0574da..3528645317 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -241,6 +241,7 @@ accelerated styles exist. * :doc:`manifoldforce ` - restrain atoms to a manifold during minimization * :doc:`meso/move ` - move mesoscopic SPH/SDPD particles in a prescribed fashion * :doc:`momentum ` - zero the linear and/or angular momentum of a group of atoms +* :doc:`momentum/chunk ` - zero the linear and/or angular momentum of a chunk of atoms * :doc:`move ` - move atoms in a prescribed fashion * :doc:`mscg ` - apply MSCG method for force-matching to generate coarse grain models * :doc:`msst ` - multi-scale shock technique (MSST) integration diff --git a/doc/src/fix_momentum.rst b/doc/src/fix_momentum.rst index 4e11859446..dbf599f29c 100644 --- a/doc/src/fix_momentum.rst +++ b/doc/src/fix_momentum.rst @@ -6,6 +6,9 @@ fix momentum command fix momentum/kk command ======================= +fix momentum/chunk command +========================== + Syntax """""" @@ -16,6 +19,16 @@ Syntax * ID, group-ID are documented in :doc:`fix ` command * momentum = style name of this fix command * N = adjust the momentum every this many timesteps + one or more keyword/value pairs may be appended +* keyword = *linear* or *angular* or *rescale* + + fix ID group-ID momentum/chunk N chunkID keyword values ... + +* ID, group-ID are documented in :doc:`fix ` command +* momentum/chunk = style name of this fix command +* N = adjust the momentum per chunk every this many timesteps +* chunkID = ID of :doc:`compute chunk/atom ` command + one or more keyword/value pairs may be appended * keyword = *linear* or *angular* or *rescale* @@ -24,9 +37,6 @@ Syntax *linear* values = xflag yflag zflag xflag,yflag,zflag = 0/1 to exclude/include each dimension *angular* values = none - - .. parsed-literal:: - *rescale* values = none Examples @@ -37,19 +47,22 @@ Examples fix 1 all momentum 1 linear 1 1 0 fix 1 all momentum 1 linear 1 1 1 rescale fix 1 all momentum 100 linear 1 1 1 angular + fix 1 all momentum/chunk 100 molchunk linear 1 1 1 angular Description """"""""""" -Zero the linear and/or angular momentum of the group of atoms every N -timesteps by adjusting the velocities of the atoms. One (or both) of -the *linear* or *angular* keywords must be specified. +Fix momentum zeroes the linear and/or angular momentum of the group of +atoms every N timesteps by adjusting the velocities of the atoms. +Fix momentum/chunk works equivalently, but operates on a per-chunk basis. + +One (or both) of the *linear* or *angular* keywords **must** be specified. If the *linear* keyword is used, the linear momentum is zeroed by -subtracting the center-of-mass velocity of the group from each atom. -This does not change the relative velocity of any pair of atoms. One -or more dimensions can be excluded from this operation by setting the -corresponding flag to 0. +subtracting the center-of-mass velocity of the group or chunk from each +atom. This does not change the relative velocity of any pair of atoms. +One or more dimensions can be excluded from this operation by setting +the corresponding flag to 0. If the *angular* keyword is used, the angular momentum is zeroed by subtracting a rotational component from each atom. @@ -60,7 +73,7 @@ to random perturbations (e.g. :doc:`fix langevin ` thermostatting). The *rescale* keyword enables conserving the kinetic energy of the group -of atoms by rescaling the velocities after the momentum was removed. +or chunk of atoms by rescaling the velocities after the momentum was removed. Note that the :doc:`velocity ` command can be used to create initial velocities with zero aggregate linear and/or angular momentum. @@ -71,15 +84,21 @@ initial velocities with zero aggregate linear and/or angular momentum. **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. No global or per-atom quantities are stored -by this fix for access by various :doc:`output commands `. +No information about this fix is written to :doc:`binary restart files +`. None of the :doc:`fix_modify ` options are +relevant to this fix. No global or per-atom quantities are stored by +this fix for access by various :doc:`output commands `. No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. Restrictions """""""""""" - none + +Fix momentum/chunk is part of the USER-MISC package. It is only enabled +if LAMMPS was built with that package. See the :doc:`Build package +` doc page for more info. + Related commands """""""""""""""" diff --git a/src/USER-MISC/fix_momentum_chunk.cpp b/src/USER-MISC/fix_momentum_chunk.cpp new file mode 100644 index 0000000000..67183ace52 --- /dev/null +++ b/src/USER-MISC/fix_momentum_chunk.cpp @@ -0,0 +1,287 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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_momentum_chunk.h" +#include +#include +#include +#include "atom.h" +#include "compute.h" +#include "compute_chunk_atom.h" +#include "compute_com_chunk.h" +#include "domain.h" +#include "group.h" +#include "error.h" +#include "force.h" +#include "modify.h" +#include "fmt/format.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- + Contributing author: Jiang Xiao (Hong Kong Polytechnic University) +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- */ + +FixMomentumChunk::FixMomentumChunk(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), + cchunk(nullptr), ccom(nullptr), cvcm(nullptr), comega(nullptr) +{ + if (narg < 5) error->all(FLERR,"Illegal fix momentum/chunk command"); + + nevery = force->inumeric(FLERR,arg[3]); + if (nevery <= 0) error->all(FLERR,"Illegal fix momentum/chunk command"); + + id_chunk = arg[4]; + int icompute = modify->find_compute(id_chunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for fix momentum/chunk"); + + id_com.clear(); + id_vcm.clear(); + id_omega.clear(); + + linear = angular = rescale = 0; + + int iarg = 5; + while (iarg < narg) { + if (strcmp(arg[iarg],"linear") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix momentum command"); + linear = 1; + xflag = force->inumeric(FLERR,arg[iarg+1]); + yflag = force->inumeric(FLERR,arg[iarg+2]); + zflag = force->inumeric(FLERR,arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg],"angular") == 0) { + angular = 1; + iarg += 1; + } else if (strcmp(arg[iarg],"rescale") == 0) { + rescale = 1; + iarg += 1; + } else error->all(FLERR,"Illegal fix momentum/chunk command"); + } + + if (linear == 0 && angular == 0) + error->all(FLERR,"Illegal fix momentum/chunk command"); + + if (linear) + if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 || + zflag < 0 || zflag > 1) + error->all(FLERR,"Illegal fix momentum/chunk command"); + + dynamic_group_allow = 0; +} + +FixMomentumChunk::~FixMomentumChunk() +{ + if (!id_com.empty()) modify->delete_compute(id_com); + if (!id_vcm.empty()) modify->delete_compute(id_vcm); + if (!id_omega.empty()) modify->delete_compute(id_omega); +} + +/* ---------------------------------------------------------------------- */ + +int FixMomentumChunk::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + mask |= POST_RUN; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixMomentumChunk::init() +{ + // current indices for idchunk and idcom + + int icompute = modify->find_compute(id_chunk); + if (icompute < 0) + error->all(FLERR,"Chunk/atom compute does not exist for fix momentum/chunk"); + cchunk = (ComputeChunkAtom *) modify->compute[icompute]; + if (strcmp(cchunk->style,"chunk/atom") != 0) + error->all(FLERR,"Fix momentum/chunk does not use chunk/atom compute"); + + // create computes dependent on chunks + + id_com = id + id_chunk + "_com"; + auto cmd = fmt::format("{} {} com/chunk {}",id_com,group->names[igroup],id_chunk); + modify->add_compute(cmd); + icompute = modify->find_compute(id_com); + ccom = (ComputeCOMChunk *) modify->compute[icompute]; + + id_vcm = id + id_chunk + "_vcm"; + cmd = fmt::format("{} {} vcm/chunk {}",id_vcm,group->names[igroup],id_chunk); + modify->add_compute(cmd); + icompute = modify->find_compute(id_vcm); + cvcm = modify->compute[icompute]; + + id_omega = id + id_chunk + "_omega"; + cmd = fmt::format("{} {} omega/chunk {}",id_omega,group->names[igroup],id_chunk); + modify->add_compute(cmd); + icompute = modify->find_compute(id_omega); + comega = modify->compute[icompute]; +} + +/* ---------------------------------------------------------------------- */ + +void FixMomentumChunk::end_of_step() +{ + // calculate per-chunk properties. + // this will also trigger a compute/update of the chunks if needed. + + ccom->compute_array(); + cvcm->compute_array(); + comega->compute_array(); + + nchunk = cchunk->nchunk; + int *ichunk = cchunk->ichunk; + double **com = ccom->array; + double **vcm = cvcm->array; + double **omega = comega->array; + + + // apply removing translational and rotational velocity from atoms in each chunk + + double **v = atom->v; + int *mask = atom->mask; + const int nlocal = atom->nlocal; + + // compute per-chunk kinetic energy before momentum removal, if needed + + double *ke_chunk_old,*ke_chunk_new,*ke_chunk_local,*factor; + if (rescale) { + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + ke_chunk_old = new double[nchunk]; + ke_chunk_local = new double[nchunk]; + memset(ke_chunk_local,0,nchunk*sizeof(double)); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + int m = ichunk[i]-1; + if (m < 0) continue; + + if (rmass) + ke_chunk_local[m] += rmass[i] * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + else + ke_chunk_local[m] += mass[type[i]] * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + } + } + MPI_Allreduce(ke_chunk_local,ke_chunk_old,nchunk,MPI_DOUBLE,MPI_SUM,world); + } + + if (linear) { + + // adjust velocities by vcm to zero linear momentum + // only adjust a component if flag is set + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + int m = ichunk[i]-1; + if (m < 0) continue; + if (xflag) v[i][0] -= vcm[m][0]; + if (yflag) v[i][1] -= vcm[m][1]; + if (zflag) v[i][2] -= vcm[m][2]; + } + } + } + + if (angular) { + + // adjust velocities to zero omega + // vnew_i = v_i - w x r_i + // must use unwrapped coords to compute r_i correctly + + double **x = atom->x; + imageint *image = atom->image; + double dx,dy,dz; + double unwrap[3]; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + int m = ichunk[i]-1; + if (m < 0) continue; + domain->unmap(x[i],image[i],unwrap); + dx = unwrap[0] - com[m][0]; + dy = unwrap[1] - com[m][1]; + dz = unwrap[2] - com[m][2]; + v[i][0] -= omega[m][1]*dz - omega[m][2]*dy; + v[i][1] -= omega[m][2]*dx - omega[m][0]*dz; + v[i][2] -= omega[m][0]*dy - omega[m][1]*dx; + } + } + } + + // compute kinetic energy after momentum removal, if needed + + if (rescale) { + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + ke_chunk_new = new double[nchunk]; + factor = new double[nchunk]; + memset(ke_chunk_local,0,nchunk*sizeof(double)); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + int m = ichunk[i]-1; + if (m < 0) continue; + + if (rmass) + ke_chunk_local[m] += rmass[i] * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + else + ke_chunk_local[m] += mass[type[i]] * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + } + } + MPI_Allreduce(ke_chunk_local,ke_chunk_new,nchunk,MPI_DOUBLE,MPI_SUM,world); + + // get scaling factors + + for (int m = 0; m < nchunk; ++m) + factor[m] = (ke_chunk_new[0] > 0.0) ? sqrt(ke_chunk_old[m]/ke_chunk_new[m]) : 1.0; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + int m = ichunk[i]-1; + if (m < 0) continue; + v[i][0] *= factor[m]; + v[i][1] *= factor[m]; + v[i][2] *= factor[m]; + } + } + delete[] factor; + delete[] ke_chunk_local; + delete[] ke_chunk_old; + delete[] ke_chunk_new; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixMomentumChunk::post_run() +{ + modify->delete_compute(id_com); + modify->delete_compute(id_vcm); + modify->delete_compute(id_omega); + id_com.clear(); + id_vcm.clear(); + id_omega.clear(); +} diff --git a/src/USER-MISC/fix_momentum_chunk.h b/src/USER-MISC/fix_momentum_chunk.h new file mode 100644 index 0000000000..5e056f2494 --- /dev/null +++ b/src/USER-MISC/fix_momentum_chunk.h @@ -0,0 +1,63 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + 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 + +FixStyle(momentum/chunk,FixMomentumChunk) + +#else + +#ifndef LMP_FIX_MOMENTUM_CHUNK_H +#define LMP_FIX_MOMENTUM_CHUNK_H + +#include "fix.h" +#include + +namespace LAMMPS_NS { + +class FixMomentumChunk : public Fix { + public: + FixMomentumChunk(class LAMMPS *, int, char **); + virtual ~FixMomentumChunk(); + int setmask(); + void init(); + void end_of_step(); + void post_run(); + + protected: + std::string id_chunk,id_com,id_vcm,id_omega; + int nchunk,linear,angular,rescale; + int xflag,yflag,zflag; + + class ComputeChunkAtom *cchunk; + class Compute *ccom,*cvcm,*comega; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Fix momentum group has no atoms + +Self-explanatory. + +*/