add support for recomputing normalization factors and finite size correction during

doc/src/compute_rdf.txt

@@ -180,9 +180,18 @@ will register an arbitrarily large spike at whatever distance they
 happen to be at, and zero everywhere else.  Coord(r) will show a step
 change from zero to one at the location of the spike in g(r).
 
+NOTE: compute rdf can handle dynamic groups and systems where atoms
+are added or removed, but this causes that certain normalization
+parameters need to be recomputed in every step and include collective
+communication operations. This will reduce performance and limit
+parallel efficiency and scaling. For systems, where only the type
+of atoms changes (e.g. when using "fix atom/swap"_fix_atom_swap.html),
+you need to explicitly request the dynamic normalization updates
+via "compute_modify dynamic yes"_compute_modify.html
+
 [Related commands:]
 
-"fix ave/time"_fix_ave_time.html
+"fix ave/time"_fix_ave_time.html, "compute_modify"_compute_modify.html
 
 [Default:]
 

src/compute_rdf.cpp

@@ -48,7 +48,6 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
 
   array_flag = 1;
   extarray = 0;
-  dynamic_group_allow = 0;
 
   nbin = force->inumeric(FLERR,arg[3]);
   if (nbin < 1) error->all(FLERR,"Illegal compute rdf command");
@@ -125,6 +124,9 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
   icount = new int[npairs];
   jcount = new int[npairs];
   duplicates = new int[npairs];
+
+  dynamic = 0;
+  natoms_old = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -150,7 +152,6 @@ ComputeRDF::~ComputeRDF()
 
 void ComputeRDF::init()
 {
-  int i,j,m;
 
   if (!force->pair && !cutflag)
     error->all(FLERR,"Compute rdf requires a pair style be defined "
@@ -184,12 +185,50 @@ void ComputeRDF::init()
   for (int i = 0; i < nbin; i++)
     array[i][0] = (i+0.5) * delr;
 
+  // initialize normalization, finite size correction, and changing atom counts
+
+  natoms_old = atom->natoms;
+  dynamic = group->dynamic[igroup];
+  if (dynamic_user) dynamic = 1;
+  init_norm();
+
+  // need an occasional half neighbor list
+  // if user specified, request a cutoff = cutoff_user + skin
+  // skin is included b/c Neighbor uses this value similar
+  //   to its cutneighmax = force cutoff + skin
+  // also, this NeighList may be used by this compute for multiple steps
+  //   (until next reneighbor), so it needs to contain atoms further
+  //   than cutoff_user apart, just like a normal neighbor list does
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->occasional = 1;
+  if (cutflag) {
+    neighbor->requests[irequest]->cut = 1;
+    neighbor->requests[irequest]->cutoff = mycutneigh;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRDF::init_list(int id, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeRDF::init_norm()
+{
+  int i,j,m;
+
   // count atoms of each type that are also in group
 
-  int *mask = atom->mask;
-  int *type = atom->type;
-  int nlocal = atom->nlocal;
-  int ntypes = atom->ntypes;
+  const int nlocal = atom->nlocal;
+  const int ntypes = atom->ntypes;
+  const int * const mask = atom->mask;
+  const int * const type = atom->type;
 
   for (i = 1; i <= ntypes; i++) typecount[i] = 0;
   for (i = 0; i < nlocal; i++)
@@ -218,30 +257,6 @@ void ComputeRDF::init()
   MPI_Allreduce(duplicates,scratch,npairs,MPI_INT,MPI_SUM,world);
   for (i = 0; i < npairs; i++) duplicates[i] = scratch[i];
   delete [] scratch;
-
-  // need an occasional half neighbor list
-  // if user specified, request a cutoff = cutoff_user + skin
-  // skin is included b/c Neighbor uses this value similar
-  //   to its cutneighmax = force cutoff + skin
-  // also, this NeighList may be used by this compute for multiple steps
-  //   (until next reneighbor), so it needs to contain atoms further
-  //   than cutoff_user apart, just like a normal neighbor list does
-
-  int irequest = neighbor->request(this,instance_me);
-  neighbor->requests[irequest]->pair = 0;
-  neighbor->requests[irequest]->compute = 1;
-  neighbor->requests[irequest]->occasional = 1;
-  if (cutflag) {
-    neighbor->requests[irequest]->cut = 1;
-    neighbor->requests[irequest]->cutoff = mycutneigh;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeRDF::init_list(int id, NeighList *ptr)
-{
-  list = ptr;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -253,6 +268,17 @@ void ComputeRDF::compute_array()
   int *ilist,*jlist,*numneigh,**firstneigh;
   double factor_lj,factor_coul;
 
+  if (natoms_old != atom->natoms) {
+    dynamic = 1;
+    natoms_old = atom->natoms;
+  }
+
+  // if the number of atoms has changed or we have a dynamic group
+  // or dynamic updates are requested (e.g. when changing atom types)
+  // we need to recompute some normalization parameters
+
+  if (dynamic) init_norm();
+
   invoked_array = update->ntimestep;
 
   // invoke half neighbor list (will copy or build if necessary)

src/compute_rdf.h

@@ -51,6 +51,8 @@ class ComputeRDF : public Compute {
   int *duplicates;
 
   class NeighList *list; // half neighbor list
+  void init_norm();
+  bigint natoms_old;
 };
 
 }