add new fix style momentum/chunk

   * :doc:`manifoldforce <fix_manifoldforce>`

   * :doc:`meso/move <fix_meso_move>`

   * :doc:`momentum (k) <fix_momentum>`

   * :doc:`momentum/chunk <fix_momentum>`

   * :doc:`move <fix_move>`

   * :doc:`mscg <fix_mscg>`

   * :doc:`msst <fix_msst>`

* :doc:`oxdna2/fene <bond_oxdna>` - same as oxdna but used with different pair styles

* :doc:`oxrna2/fene <bond_oxdna>` - modified FENE bond suitable for RNA modeling

* :doc:`quartic <bond_quartic>` - breakable quartic bond

* :doc:`special <bond_special>` - enable special bond exclusions for 1-5 pairs and beyond

* :doc:`table <bond_table>` - tabulated by bond length

----------

* :doc:`manifoldforce <fix_manifoldforce>` - restrain atoms to a manifold during minimization

* :doc:`meso/move <fix_meso_move>` - move mesoscopic SPH/SDPD particles in a prescribed fashion

* :doc:`momentum <fix_momentum>` - zero the linear and/or angular momentum of a group of atoms

* :doc:`momentum/chunk <fix_momentum>` - zero the linear and/or angular momentum of a chunk of atoms

* :doc:`move <fix_move>` - move atoms in a prescribed fashion

* :doc:`mscg <fix_mscg>` - apply MSCG method for force-matching to generate coarse grain models

* :doc:`msst <fix_msst>` - multi-scale shock technique (MSST) integration

fix momentum/kk command

=======================

fix momentum/chunk command

==========================

Syntax

""""""

* ID, group-ID are documented in :doc:`fix <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 <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 <compute_chunk_atom>` command

  one or more keyword/value pairs may be appended

* keyword = *linear* or *angular* or *rescale*

       *linear* values = xflag yflag zflag

         xflag,yflag,zflag = 0/1 to exclude/include each dimension

       *angular* values = none

  .. parsed-literal::

       *rescale* values = none

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.

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 <velocity>` command can be used to create

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 <restart>`.  None of the :doc:`fix_modify <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 <Howto_output>`.

No information about this fix is written to :doc:`binary restart files

<restart>`.  None of the :doc:`fix_modify <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 <Howto_output>`.

No parameter of this fix can be used with the *start/stop* keywords of

the :doc:`run <run>` command.  This fix is not invoked during :doc:`energy minimization <minimize>`.

the :doc:`run <run>` command.  This fix is not invoked during

:doc:`energy minimization <minimize>`.

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

<Build_package>` doc page for more info.

Related commands

""""""""""""""""

/* ----------------------------------------------------------------------

   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 <mpi.h>

#include <cmath>

#include <cstring>

#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();

}