From 59b1c3ca603fd0dd399bc763176f03cf849b175e Mon Sep 17 00:00:00 2001 From: Shun Sakuraba Date: Fri, 22 Nov 2019 20:05:32 +0900 Subject: [PATCH 3/3] Fixed multiple issues with gmx spatial gmx spatial had multiple problems which is now fixed by this commit: (1) Having negative iIGNOREOUTER (default -1) led to incorrect number of grid points (2) The coordinates of the grid points were incorrect irrespective to the iIGNOREOUTER (3) Norm calculation was incorrect (4) Norm may become negative due to overflow --- src/gromacs/gmxana/gmx_spatial.cpp | 74 +++++++++++++----------------- 1 file changed, 33 insertions(+), 41 deletions(-) diff --git a/src/gromacs/gmxana/gmx_spatial.cpp b/src/gromacs/gmxana/gmx_spatial.cpp index c28bcc714b..39221b6a84 100644 --- a/src/gromacs/gmxana/gmx_spatial.cpp +++ b/src/gromacs/gmxana/gmx_spatial.cpp @@ -149,8 +149,10 @@ int gmx_spatial(int argc, char *argv[]) int nbin[3]; FILE *flp; int x, y, z, minx, miny, minz, maxx, maxy, maxz; + int boutput[DIM], eoutput[DIM]; int numfr, numcu; - int tot, maxval, minval; + int maxval, minval; + int64_t tot; double norm; gmx_output_env_t *oenv; gmx_rmpbc_t gpbc = nullptr; @@ -306,19 +308,33 @@ int gmx_spatial(int argc, char *argv[]) if (!bCUTDOWN) { minx = miny = minz = 0; - maxx = nbin[XX]; - maxy = nbin[YY]; - maxz = nbin[ZZ]; + maxx = nbin[XX] - 1; + maxy = nbin[YY] - 1; + maxz = nbin[ZZ] - 1; } + boutput[0] = minx + std::max(iIGNOREOUTER, 0); + boutput[1] = miny + std::max(iIGNOREOUTER, 0); + boutput[2] = minz + std::max(iIGNOREOUTER, 0); + eoutput[0] = maxx - std::max(iIGNOREOUTER, 0); + eoutput[1] = maxy - std::max(iIGNOREOUTER, 0); + eoutput[2] = maxz - std::max(iIGNOREOUTER, 0); + /* OUTPUT */ flp = gmx_ffopen("grid.cube", "w"); fprintf(flp, "Spatial Distribution Function\n"); fprintf(flp, "test\n"); - fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", nidxp, (MINBIN[XX]+(minx+iIGNOREOUTER)*rBINWIDTH)*10./bohr, (MINBIN[YY]+(miny+iIGNOREOUTER)*rBINWIDTH)*10./bohr, (MINBIN[ZZ]+(minz+iIGNOREOUTER)*rBINWIDTH)*10./bohr); - fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxx-minx+1-(2*iIGNOREOUTER), rBINWIDTH*10./bohr, 0., 0.); - fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxy-miny+1-(2*iIGNOREOUTER), 0., rBINWIDTH*10./bohr, 0.); - fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxz-minz+1-(2*iIGNOREOUTER), 0., 0., rBINWIDTH*10./bohr); + /* + Values in .cube file represent the density at the grid point. + Corresponding coordinates to the binned value is the center of the bin, i.e. + 0.5 to the bin index. + */ + fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", nidxp, + (MINBIN[XX]+(boutput[XX] + 0.5)*rBINWIDTH)*10./bohr, + (MINBIN[YY]+(boutput[YY] + 0.5)*rBINWIDTH)*10./bohr, + (MINBIN[ZZ]+(boutput[ZZ] + 0.5)*rBINWIDTH)*10./bohr); + fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", eoutput[XX] - boutput[XX], rBINWIDTH*10./bohr, 0., 0.); + fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", eoutput[YY] - boutput[YY], 0., rBINWIDTH*10./bohr, 0.); + fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", eoutput[ZZ] - boutput[ZZ], 0., 0., rBINWIDTH*10./bohr); for (i = 0; i < nidxp; i++) { v = 2; @@ -376,24 +392,12 @@ int gmx_spatial(int argc, char *argv[]) minval = 999; maxval = 0; - for (k = 0; k < nbin[XX]; k++) + for (k = boutput[XX]; k < eoutput[XX]; k++) { - if (k < minx+iIGNOREOUTER || k > maxx-iIGNOREOUTER) + for (j = boutput[YY]; j < eoutput[YY]; j++) { - continue; - } - for (j = 0; j < nbin[YY]; j++) - { - if (j < miny+iIGNOREOUTER || j > maxy-iIGNOREOUTER) + for (i = boutput[ZZ]; i < eoutput[ZZ]; i++) { - continue; - } - for (i = 0; i < nbin[ZZ]; i++) - { - if (i < minz+iIGNOREOUTER || i > maxz-iIGNOREOUTER) - { - continue; - } tot += bin[k][j][i]; if (bin[k][j][i] > maxval) { @@ -407,34 +411,22 @@ int gmx_spatial(int argc, char *argv[]) } } - numcu = (maxx-minx+1-(2*iIGNOREOUTER))*(maxy-miny+1-(2*iIGNOREOUTER))*(maxz-minz+1-(2*iIGNOREOUTER)); + numcu = (eoutput[XX] - boutput[XX]) * (eoutput[YY] - boutput[YY]) * (eoutput[ZZ] - boutput[ZZ]); if (bCALCDIV) { - norm = static_cast(numcu*numfr)/tot; + norm = static_cast(static_cast(numcu)*numfr)/tot; } else { norm = 1.0; } - for (k = 0; k < nbin[XX]; k++) + for (k = boutput[XX]; k < eoutput[XX]; k++) { - if (k < minx+iIGNOREOUTER || k > maxx-iIGNOREOUTER) - { - continue; - } - for (j = 0; j < nbin[YY]; j++) + for (j = boutput[YY]; j < eoutput[YY]; j++) { - if (j < miny+iIGNOREOUTER || j > maxy-iIGNOREOUTER) - { - continue; - } - for (i = 0; i < nbin[ZZ]; i++) + for (i = boutput[ZZ]; i < eoutput[ZZ]; i++) { - if (i < minz+iIGNOREOUTER || i > maxz-iIGNOREOUTER) - { - continue; - } fprintf(flp, "%12.6f ", static_cast(norm*bin[k][j][i])/numfr); } fprintf(flp, "\n"); @@ -451,7 +443,7 @@ int gmx_spatial(int argc, char *argv[]) else { printf("grid.cube contains counts per frame in all %d cubes\n", numcu); - printf("Raw data: average %le, min %le, max %le\n", 1.0/norm, static_cast(minval)/numfr, static_cast(maxval)/numfr); + printf("Raw data: average %le, min %le, max %le\n", static_cast(tot)/numfr/numcu, static_cast(minval)/numfr, static_cast(maxval)/numfr); } return 0; -- 2.17.1