Merge pull request #1529 from rupertnash/compute_momentum

add compute momentum command

doc/src/Commands_compute.txt

@@ -81,6 +81,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "meso/e/atom"_compute_meso_e_atom.html,
 "meso/rho/atom"_compute_meso_rho_atom.html,
 "meso/t/atom"_compute_meso_t_atom.html,
+"momentum"_compute_momentum.html,
 "msd"_compute_msd.html,
 "msd/chunk"_compute_msd_chunk.html,
 "msd/nongauss"_compute_msd_nongauss.html,

doc/src/compute.txt

@@ -228,6 +228,7 @@ compute"_Commands_compute.html doc page are followed by one or more of
 "meso/e/atom"_compute_meso_e_atom.html - per-atom internal energy of Smooth-Particle Hydrodynamics atoms
 "meso/rho/atom"_compute_meso_rho_atom.html - per-atom mesoscopic density of Smooth-Particle Hydrodynamics atoms
 "meso/t/atom"_compute_meso_t_atom.html - per-atom internal temperature of Smooth-Particle Hydrodynamics atoms
+"momentum"_compute_momentum.html - translational momentum
 "msd"_compute_msd.html - mean-squared displacement of group of atoms
 "msd/chunk"_compute_msd_chunk.html - mean-squared displacement for each chunk
 "msd/nongauss"_compute_msd_nongauss.html - MSD and non-Gaussian parameter of group of atoms

doc/src/compute_momentum.txt

@@ -0,0 +1,49 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+compute momentum command :h3
+
+[Syntax:]
+
+compute ID group-ID momentum :pre
+
+ID, group-ID are documented in "compute"_compute.html command
+momentum = style name of this compute command :ul
+
+[Examples:]
+
+compute 1 all momentum :pre
+
+[Description:]
+
+Define a computation that calculates the translational momentum
+of a group of particles.
+
+The momentum of each particles is computed as m v, where m and v are
+the mass and velocity of the particle.
+
+[Output info:]
+
+This compute calculates a global vector (the summed momentum) of
+length 3. This value can be used by any command that uses a global
+vector value from a compute as input. See the "Howto
+output"_Howto_output.html doc page for an overview of LAMMPS output
+options.
+
+The vector value calculated by this compute is "extensive". The vector
+value will be in mass*velocity "units"_units.html.
+
+[Restrictions:]
+
+This compute is part of the USER-MISC package.  It is only enabled if
+LAMMPS was built with that package.  See the "Build
+package"_Build_package.html doc page for more info.
+
+[Related commands:]
+
+[Default:] none

doc/src/computes.txt

@@ -58,6 +58,7 @@ Computes :h1
    compute_meso_e_atom
    compute_meso_rho_atom
    compute_meso_t_atom
+   compute_momentum
    compute_msd
    compute_msd_chunk
    compute_msd_nongauss

doc/src/lammps.book

@@ -469,6 +469,7 @@ compute_ke_rigid.html
 compute_meso_e_atom.html
 compute_meso_rho_atom.html
 compute_meso_t_atom.html
+compute_momentum.html
 compute_msd.html
 compute_msd_chunk.html
 compute_msd_nongauss.html

examples/USER/misc/momentum/in.momentum

@@ -0,0 +1,65 @@
+# Test compute momentum with by comparing it's output with compute reduce
+
+# Script will output the two computations of the total momentum as
+# thermo vars - they should be identical
+
+# 3D LJ Poiseuille flow simulation
+dimension	3
+boundary	p p p
+
+atom_style	atomic
+neighbor	0.3 bin
+neigh_modify	delay 5
+
+# create geometry
+lattice		hcp 0.5
+region		box block 0 20 0 10 0 10
+create_box	3 box
+create_atoms	1 box
+
+mass		1 1.0
+mass		2 1.0
+mass		3 1.0
+
+# LJ potentials
+pair_style	lj/cut 1.0
+pair_coeff	* * 1.0 1.0 1.0
+
+# define groups
+region	     1 block INF INF INF 1.0 INF INF
+group	     lower region 1
+region	     2 block INF INF 9.0 INF INF INF
+group	     upper region 2
+group	     boundary union lower upper
+group	     liquid subtract all boundary
+
+set	     group lower type 2
+set	     group upper type 3
+
+# temperature settings
+compute      mobile liquid temp
+velocity     liquid create 1.0 100 temp mobile
+fix          1 all nve
+fix          2 liquid temp/rescale 200 1.0 2.0 0.02 1.0
+fix_modify   2 temp mobile
+
+velocity     boundary set 0.0 0.0 0.0
+fix	     4 lower setforce 0.0 0.0 0.0
+fix	     5 upper setforce 0.0 0.0 0.0
+fix	     6 liquid addforce 0.1 0.0 0.0
+
+# Compute total momentum
+compute mom liquid momentum
+
+# Comparison (momentum variables + compute reduce)
+variable mx atom mass*vx
+variable my atom mass*vy
+variable mz atom mass*vz
+compute mom_red liquid reduce sum v_mx v_my v_mz
+
+# Run
+timestep	0.003
+thermo		1000
+thermo_style custom step c_mom_red[*] c_mom[*]
+
+run		100000

src/USER-MISC/README

@@ -32,6 +32,7 @@ compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017
 compute entropy/atom, Pablo Piaggi (EPFL), pablo.piaggi at epfl.ch, 15 June 2018
 compute gyration/shape, Evangelos Voyiatzis, evoyiatzis at gmail.com, 25 July 2019
 compute pressure/cylinder, Cody K. Addington (NCSU), , 2 Oct 2018
+compute momentum, Rupert Nash (University of Edinburgh), r.nash at epcc.ed.ac.uk, 28 June 2019
 compute stress/mop, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
 compute stress/mop/profile, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
 compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11

src/USER-MISC/compute_momentum.cpp

@@ -0,0 +1,56 @@
+
+#include "compute_momentum.h"
+#include <mpi.h>
+#include "atom.h"
+#include "error.h"
+#include "update.h"
+
+using namespace LAMMPS_NS;
+
+ComputeMomentum::ComputeMomentum(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all(FLERR,"Illegal compute momentum command");
+
+  vector_flag = 1;
+  size_vector = 3;
+  extvector = 1;
+  vector = new double[size_vector];
+}
+
+ComputeMomentum::~ComputeMomentum() {
+  delete[] vector;
+}
+
+void ComputeMomentum::init()
+{
+}
+
+void ComputeMomentum::compute_vector()
+{
+  invoked_vector = update->ntimestep;
+
+  double **v = atom->v;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *mask = atom->mask;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  double mom[3] = {0.0, 0.0, 0.0};
+
+  if (rmass) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        for(int j = 0; j < 3; ++j)
+          mom[j] += rmass[i] * v[i][j];
+      }
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+        for(int j = 0; j < 3; ++j)
+          mom[j] += mass[type[i]] * v[i][j];
+  }
+
+  MPI_Allreduce(&mom, vector, 3, MPI_DOUBLE, MPI_SUM, world);
+}

src/USER-MISC/compute_momentum.h

@@ -0,0 +1,36 @@
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(momentum,ComputeMomentum)
+
+#else
+
+#ifndef LMP_COMPUTE_MOMENTUM_H
+#define LMP_COMPUTE_MOMENTUM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeMomentum : public Compute {
+ public:
+  ComputeMomentum(class LAMMPS *, int, char **);
+  virtual ~ComputeMomentum();
+
+  virtual void init();
+  virtual void compute_vector();
+};
+
+}
+
+#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.
+
+*/