Added compute_rdf.cpp

src/compute_rdf.cpp

@@ -45,7 +45,7 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
   hist(NULL), histall(NULL), typecount(NULL), icount(NULL), jcount(NULL),
   duplicates(NULL)
 {
-  if (narg < 4 || (narg-4) % 2) error->all(FLERR,"Illegal compute rdf command");
+  if (narg < 4) error->all(FLERR,"Illegal compute rdf command");
 
   array_flag = 1;
   extarray = 0;
@@ -77,7 +77,10 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
   // pairwise args
 
   if (nargpair == 0) npairs = 1;
-  else npairs = nargpair/2;
+  else { 
+    if (nargpair % 2) error->all(FLERR,"Illegal compute rdf command");
+    npairs = nargpair/2;
+  }
 
   size_array_rows = nbin;
   size_array_cols = 1 + 2*npairs;
@@ -93,17 +96,13 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
   if (nargpair == 0) {
     ilo[0] = 1; ihi[0] = ntypes;
     jlo[0] = 1; jhi[0] = ntypes;
-    npairs = 1;
-
   } else {
-    npairs = 0;
     iarg = 4;
-    while (iarg < 4+nargpair) {
-      force->bounds(FLERR,arg[iarg],atom->ntypes,ilo[npairs],ihi[npairs]);
-      force->bounds(FLERR,arg[iarg+1],atom->ntypes,jlo[npairs],jhi[npairs]);
-      if (ilo[npairs] > ihi[npairs] || jlo[npairs] > jhi[npairs])
+    for (int ipair = 0; ipair < npairs; ipair++) {
+      force->bounds(FLERR,arg[iarg],atom->ntypes,ilo[ipair],ihi[ipair]);
+      force->bounds(FLERR,arg[iarg+1],atom->ntypes,jlo[ipair],jhi[ipair]);
+      if (ilo[ipair] > ihi[ipair] || jlo[ipair] > jhi[ipair])
         error->all(FLERR,"Illegal compute rdf command");
-      npairs++;
       iarg += 2;
     }
   }
@@ -113,10 +112,13 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
     for (j = 1; j <= ntypes; j++)
       nrdfpair[i][j] = 0;
 
+  int ihisto;
   for (int m = 0; m < npairs; m++)
     for (i = ilo[m]; i <= ihi[m]; i++)
-      for (j = jlo[m]; j <= jhi[m]; j++)
-        rdfpair[nrdfpair[i][j]++][i][j] = m;
+      for (j = jlo[m]; j <= jhi[m]; j++) {
+        ihisto = nrdfpair[i][j]++;
+        rdfpair[ihisto][i][j] = m;
+      }
 
   memory->create(hist,npairs,nbin,"rdf:hist");
   memory->create(histall,npairs,nbin,"rdf:histall");
@@ -347,13 +349,15 @@ void ComputeRDF::compute_array()
       ibin = static_cast<int> (r*delrinv);
       if (ibin >= nbin) continue;
 
-      if (ipair)
-        for (ihisto = 0; ihisto < ipair; ihisto++)
-          hist[rdfpair[ihisto][itype][jtype]][ibin] += 1.0;
+      for (ihisto = 0; ihisto < ipair; ihisto++) {
+        m = rdfpair[ihisto][itype][jtype];
+        hist[m][ibin] += 1.0;
+      }
       if (newton_pair || j < nlocal) {
-        if (jpair)
-          for (ihisto = 0; ihisto < jpair; ihisto++)
-            hist[rdfpair[ihisto][jtype][itype]][ibin] += 1.0;
+        for (ihisto = 0; ihisto < jpair; ihisto++) {
+          m = rdfpair[ihisto][jtype][itype];
+          hist[m][ibin] += 1.0;
+        }
       }
     }
   }