diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp index fe4be0cd3b..7c5f433bb5 100644 --- a/src/compute_rdf.cpp +++ b/src/compute_rdf.cpp @@ -95,7 +95,6 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) : typecount = new int[ntypes+1]; icount = new int[npairs]; jcount = new int[npairs]; - duplicates = new int[npairs]; } /* ---------------------------------------------------------------------- */ @@ -114,14 +113,13 @@ ComputeRDF::~ComputeRDF() delete [] typecount; delete [] icount; delete [] jcount; - delete [] duplicates; } /* ---------------------------------------------------------------------- */ void ComputeRDF::init() { - int i,j,m; + int i,m; if (force->pair) delr = force->pair->cutforce / nbin; else error->all(FLERR,"Compute rdf requires a pair style be defined"); @@ -145,17 +143,12 @@ void ComputeRDF::init() // 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]; @@ -163,11 +156,10 @@ void ComputeRDF::init() 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; // need an occasional half neighbor list + int irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->pair = 0; neighbor->requests[irequest]->compute = 1; @@ -273,12 +265,10 @@ void ComputeRDF::compute_array() MPI_Allreduce(hist[0],histall[0],npairs*nbin,MPI_DOUBLE,MPI_SUM,world); // convert counts to g(r) and coord(r) and copy into output array - // vfrac = fraction of volume in shell m - // npairs = number of pairs, corrected for duplicates - // duplicates = pairs in which both atoms are the same + // nideal = # of J atoms surrounding single I atom in a single bin + // assuming J atoms are at uniform density - double constant,vfrac,gr,ncoord,rlower,rupper; - int pairtot; + double constant,nideal,gr,ncoord,rlower,rupper; if (domain->dimension == 3) { constant = 4.0*MY_PI / (3.0*domain->xprd*domain->yprd*domain->zprd); @@ -288,13 +278,12 @@ void ComputeRDF::compute_array() for (ibin = 0; ibin < nbin; ibin++) { rlower = ibin*delr; rupper = (ibin+1)*delr; - vfrac = constant * (rupper*rupper*rupper - rlower*rlower*rlower); - pairtot = icount[m] * jcount[m] - duplicates[m]; - if (vfrac * pairtot != 0.0) - gr = histall[m][ibin] / (vfrac * pairtot); + nideal = constant * + (rupper*rupper*rupper - rlower*rlower*rlower) * jcount[m]; + if (icount[m]*nideal != 0.0) + gr = histall[m][ibin] / (icount[m]*nideal); else gr = 0.0; - if (icount[m] != 0) - ncoord += gr * vfrac * pairtot / icount[m]; + ncoord += gr*nideal; array[ibin][1+2*m] = gr; array[ibin][2+2*m] = ncoord; } @@ -308,13 +297,11 @@ void ComputeRDF::compute_array() for (ibin = 0; ibin < nbin; ibin++) { rlower = ibin*delr; rupper = (ibin+1)*delr; - vfrac = constant * (rupper*rupper - rlower*rlower); - pairtot = icount[m] * jcount[m] - duplicates[m]; - if (vfrac * pairtot != 0.0) - gr = histall[m][ibin] / (vfrac * pairtot); + nideal = constant * (rupper*rupper - rlower*rlower) * jcount[m]; + if (icount[m]*nideal != 0.0) + gr = histall[m][ibin] / (icount[m]*nideal); else gr = 0.0; - if (icount[m] != 0) - ncoord += gr * vfrac * pairtot / icount[m]; + ncoord += gr*nideal; array[ibin][1+2*m] = gr; array[ibin][2+2*m] = ncoord; } diff --git a/src/compute_rdf.h b/src/compute_rdf.h index 0da56609e4..794fcc9ee6 100644 --- a/src/compute_rdf.h +++ b/src/compute_rdf.h @@ -44,7 +44,7 @@ class ComputeRDF : public Compute { double **histall; // summed histogram bins across all procs int *typecount; - int *icount,*jcount,*duplicates; + int *icount,*jcount; class NeighList *list; // half neighbor list };