Merge pull request #566 from akohlmey/compute_rdf_dynamic

Dynamic group and normalization support for compute rdf

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,40 +185,12 @@ void ComputeRDF::init()
   for (int i = 0; i < nbin; i++)
     array[i][0] = (i+0.5) * delr;
 
-  // count atoms of each type that are also in group
+  // initialize normalization, finite size correction, and changing atom counts
 
-  int *mask = atom->mask;
-  int *type = atom->type;
-  int nlocal = atom->nlocal;
-  int ntypes = atom->ntypes;
-
-  for (i = 1; i <= ntypes; i++) typecount[i] = 0;
-  for (i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) typecount[type[i]]++;
-
-  // icount = # of I atoms participating in I,J pairs for each histogram
-  // jcount = # of J atoms participating in I,J pairs for each histogram
-  // duplicates = # of atoms in both groups I and J for each histogram
-
-  for (m = 0; m < npairs; m++) {
-    icount[m] = 0;
-    for (i = ilo[m]; i <= ihi[m]; i++) icount[m] += typecount[i];
-    jcount[m] = 0;
-    for (i = jlo[m]; i <= jhi[m]; i++) jcount[m] += typecount[i];
-    duplicates[m] = 0;
-    for (i = ilo[m]; i <= ihi[m]; i++)
-      for (j = jlo[m]; j <= jhi[m]; j++)
-        if (i == j) duplicates[m] += typecount[i];
-  }
-
-  int *scratch = new int[npairs];
-  MPI_Allreduce(icount,scratch,npairs,MPI_INT,MPI_SUM,world);
-  for (i = 0; i < npairs; i++) icount[i] = scratch[i];
-  MPI_Allreduce(jcount,scratch,npairs,MPI_INT,MPI_SUM,world);
-  for (i = 0; i < npairs; i++) jcount[i] = scratch[i];
-  MPI_Allreduce(duplicates,scratch,npairs,MPI_INT,MPI_SUM,world);
-  for (i = 0; i < npairs; i++) duplicates[i] = scratch[i];
-  delete [] scratch;
+  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
@@ -246,6 +219,48 @@ void ComputeRDF::init_list(int id, NeighList *ptr)
 
 /* ---------------------------------------------------------------------- */
 
+void ComputeRDF::init_norm()
+{
+  int i,j,m;
+
+  // count atoms of each type that are also in group
+
+  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++)
+    if (mask[i] & groupbit) typecount[type[i]]++;
+
+  // icount = # of I atoms participating in I,J pairs for each histogram
+  // jcount = # of J atoms participating in I,J pairs for each histogram
+  // duplicates = # of atoms in both groups I and J for each histogram
+
+  for (m = 0; m < npairs; m++) {
+    icount[m] = 0;
+    for (i = ilo[m]; i <= ihi[m]; i++) icount[m] += typecount[i];
+    jcount[m] = 0;
+    for (i = jlo[m]; i <= jhi[m]; i++) jcount[m] += typecount[i];
+    duplicates[m] = 0;
+    for (i = ilo[m]; i <= ihi[m]; i++)
+      for (j = jlo[m]; j <= jhi[m]; j++)
+        if (i == j) duplicates[m] += typecount[i];
+  }
+
+  int *scratch = new int[npairs];
+  MPI_Allreduce(icount,scratch,npairs,MPI_INT,MPI_SUM,world);
+  for (i = 0; i < npairs; i++) icount[i] = scratch[i];
+  MPI_Allreduce(jcount,scratch,npairs,MPI_INT,MPI_SUM,world);
+  for (i = 0; i < npairs; i++) jcount[i] = scratch[i];
+  MPI_Allreduce(duplicates,scratch,npairs,MPI_INT,MPI_SUM,world);
+  for (i = 0; i < npairs; i++) duplicates[i] = scratch[i];
+  delete [] scratch;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ComputeRDF::compute_array()
 {
   int i,j,m,ii,jj,inum,jnum,itype,jtype,ipair,jpair,ibin,ihisto;
@@ -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;
 };
 
 }