ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14082 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-09-28 12:37:59 +00:00
parent 4d4d56175a
commit a6d7932ef4
2 changed files with 15 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
src/compute_rdf.cpp
View File

@ -95,7 +95,6 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
typecount = new int[ntypes+1]; typecount = new int[ntypes+1];
icount = new int[npairs]; icount = new int[npairs];
jcount = new int[npairs]; jcount = new int[npairs];
duplicates = new int[npairs];
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -114,14 +113,13 @@ ComputeRDF::~ComputeRDF()
delete [] typecount; delete [] typecount;
delete [] icount; delete [] icount;
delete [] jcount; delete [] jcount;
delete [] duplicates;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputeRDF::init() void ComputeRDF::init()
{ {
int i,j,m; int i,m;
if (force->pair) delr = force->pair->cutforce / nbin; if (force->pair) delr = force->pair->cutforce / nbin;
else error->all(FLERR,"Compute rdf requires a pair style be defined"); 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 // icount = # of I atoms participating in I,J pairs for each histogram
// jcount = # of J 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++) { for (m = 0; m < npairs; m++) {
icount[m] = 0; icount[m] = 0;
for (i = ilo[m]; i <= ihi[m]; i++) icount[m] += typecount[i]; for (i = ilo[m]; i <= ihi[m]; i++) icount[m] += typecount[i];
jcount[m] = 0; jcount[m] = 0;
for (i = jlo[m]; i <= jhi[m]; i++) jcount[m] += typecount[i]; 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]; int *scratch = new int[npairs];
@ -163,11 +156,10 @@ void ComputeRDF::init()
for (i = 0; i < npairs; i++) icount[i] = scratch[i]; for (i = 0; i < npairs; i++) icount[i] = scratch[i];
MPI_Allreduce(jcount,scratch,npairs,MPI_INT,MPI_SUM,world); MPI_Allreduce(jcount,scratch,npairs,MPI_INT,MPI_SUM,world);
for (i = 0; i < npairs; i++) jcount[i] = scratch[i]; 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; delete [] scratch;
// need an occasional half neighbor list // need an occasional half neighbor list
int irequest = neighbor->request(this,instance_me); int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0; neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1; 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); 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 // convert counts to g(r) and coord(r) and copy into output array
// vfrac = fraction of volume in shell m // nideal = # of J atoms surrounding single I atom in a single bin
// npairs = number of pairs, corrected for duplicates // assuming J atoms are at uniform density
// duplicates = pairs in which both atoms are the same
double constant,vfrac,gr,ncoord,rlower,rupper; double constant,nideal,gr,ncoord,rlower,rupper;
int pairtot;
if (domain->dimension == 3) { if (domain->dimension == 3) {
constant = 4.0*MY_PI / (3.0*domain->xprd*domain->yprd*domain->zprd); 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++) { for (ibin = 0; ibin < nbin; ibin++) {
rlower = ibin*delr; rlower = ibin*delr;
rupper = (ibin+1)*delr; rupper = (ibin+1)*delr;
vfrac = constant * (rupper*rupper*rupper - rlower*rlower*rlower); nideal = constant *
pairtot = icount[m] * jcount[m] - duplicates[m]; (rupper*rupper*rupper - rlower*rlower*rlower) * jcount[m];
if (vfrac * pairtot != 0.0) if (icount[m]*nideal != 0.0)
gr = histall[m][ibin] / (vfrac * pairtot); gr = histall[m][ibin] / (icount[m]*nideal);
else gr = 0.0; else gr = 0.0;
if (icount[m] != 0) ncoord += gr*nideal;
ncoord += gr * vfrac * pairtot / icount[m];
array[ibin][1+2*m] = gr; array[ibin][1+2*m] = gr;
array[ibin][2+2*m] = ncoord; array[ibin][2+2*m] = ncoord;
} }
@ -308,13 +297,11 @@ void ComputeRDF::compute_array()
for (ibin = 0; ibin < nbin; ibin++) { for (ibin = 0; ibin < nbin; ibin++) {
rlower = ibin*delr; rlower = ibin*delr;
rupper = (ibin+1)*delr; rupper = (ibin+1)*delr;
vfrac = constant * (rupper*rupper - rlower*rlower); nideal = constant * (rupper*rupper - rlower*rlower) * jcount[m];
pairtot = icount[m] * jcount[m] - duplicates[m]; if (icount[m]*nideal != 0.0)
if (vfrac * pairtot != 0.0) gr = histall[m][ibin] / (icount[m]*nideal);
gr = histall[m][ibin] / (vfrac * pairtot);
else gr = 0.0; else gr = 0.0;
if (icount[m] != 0) ncoord += gr*nideal;
ncoord += gr * vfrac * pairtot / icount[m];
array[ibin][1+2*m] = gr; array[ibin][1+2*m] = gr;
array[ibin][2+2*m] = ncoord; array[ibin][2+2*m] = ncoord;
} }

2
src/compute_rdf.h
View File

@ -44,7 +44,7 @@ class ComputeRDF : public Compute {
double **histall; // summed histogram bins across all procs double **histall; // summed histogram bins across all procs
int *typecount; int *typecount;
int *icount,*jcount,*duplicates; int *icount,*jcount;
class NeighList *list; // half neighbor list class NeighList *list; // half neighbor list
}; };