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

reverse normalization between type pairs if the types were swapped on input

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-03-06 13:15:09 -05:00
parent 09c8dc07d8
commit 916ab55a31
2 changed files with 24 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

33
src/compute_rdf.cpp
View File

@ -43,8 +43,8 @@ using namespace MathConst;
ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
rdfpair(nullptr), nrdfpair(nullptr), ilo(nullptr), ihi(nullptr), jlo(nullptr), jhi(nullptr),
hist(nullptr), histall(nullptr), typecount(nullptr), icount(nullptr), jcount(nullptr),
duplicates(nullptr)
rev(nullptr), hist(nullptr), histall(nullptr), typecount(nullptr),
icount(nullptr), jcount(nullptr), duplicates(nullptr)
{
if (narg < 4) utils::missing_cmd_args(FLERR,"compute rdf", error);
@ -95,10 +95,12 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
ihi = new int[npairs];
jlo = new int[npairs];
jhi = new int[npairs];
rev = new int[npairs];
if (!nargpair) {
ilo[0] = 1; ihi[0] = ntypes;
jlo[0] = 1; jhi[0] = ntypes;
rev[0] = 0;
} else {
iarg = 4;
for (int ipair = 0; ipair < npairs; ipair++) {
@ -110,11 +112,12 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) :
if ( (ilo[ipair] == ihi[ipair]) &&
(jlo[ipair] == jhi[ipair]) &&
(ilo[ipair] > jlo[ipair]) ) {
rev[ipair] = 1;
jlo[ipair] = ihi[ipair];
ilo[ipair] = jhi[ipair];
ihi[ipair] = ilo[ipair];
jhi[ipair] = jlo[ipair];
}
} else rev[ipair] = 0;
iarg += 2;
}
@ -381,18 +384,28 @@ void ComputeRDF::compute_array()
constant = 4.0*MY_PI / (3.0*domain->xprd*domain->yprd*domain->zprd);
for (m = 0; m < npairs; m++) {
normfac = (icount[m] > 0) ? static_cast<double>(jcount[m])
- static_cast<double>(duplicates[m])/icount[m] : 0.0;
if (rev[m]) { // swap i and j because they were entered in different order
normfac = (jcount[m] > 0) ? static_cast<double>(icount[m])
- static_cast<double>(duplicates[m])/jcount[m] : 0.0;
} else {
normfac = (icount[m] > 0) ? static_cast<double>(jcount[m])
- static_cast<double>(duplicates[m])/icount[m] : 0.0;
}
ncoord = 0.0;
for (ibin = 0; ibin < nbin; ibin++) {
rlower = ibin*delr;
rupper = (ibin+1)*delr;
vfrac = constant * (rupper*rupper*rupper - rlower*rlower*rlower);
if (vfrac * normfac != 0.0)
gr = histall[m][ibin] / (vfrac * normfac * icount[m]);
else gr = 0.0;
if (icount[m] != 0)
ncoord += gr * vfrac * normfac;
gr = 0.0;
if (vfrac * normfac != 0.0) {
if (rev[m]) { // swap i and j because they were entered in different order
gr = histall[m][ibin] / (vfrac * normfac * jcount[m]);
if (jcount[m] != 0) ncoord += gr * vfrac * normfac;
} else {
gr = histall[m][ibin] / (vfrac * normfac * icount[m]);
if (icount[m] != 0) ncoord += gr * vfrac * normfac;
}
}
array[ibin][1+2*m] = gr;
array[ibin][2+2*m] = ncoord;
}

2
src/compute_rdf.h
View File

@ -41,7 +41,7 @@ class ComputeRDF : public Compute {
double mycutneigh; // user-specified cutoff + neighbor skin
int ***rdfpair; // map 2 type pair to rdf pair for each histo
int **nrdfpair; // # of histograms for each type pair
int *ilo, *ihi, *jlo, *jhi;
int *ilo, *ihi, *jlo, *jhi, *rev;
double **hist; // histogram bins
double **histall; // summed histogram bins across all procs