Project

General

Profile

0003-Fixed-multiple-issues-with-gmx-spatial.patch

Fixes (3)-(6) - Shun Sakuraba, 11/25/2019 04:01 PM

View differences:

src/gromacs/gmxana/gmx_spatial.cpp
149 149
    int               nbin[3];
150 150
    FILE             *flp;
151 151
    int               x, y, z, minx, miny, minz, maxx, maxy, maxz;
152
    int               boutput[DIM], eoutput[DIM];
152 153
    int               numfr, numcu;
153
    int               tot, maxval, minval;
154
    int               maxval, minval;
155
    int64_t           tot;
154 156
    double            norm;
155 157
    gmx_output_env_t *oenv;
156 158
    gmx_rmpbc_t       gpbc = nullptr;
......
306 308
    if (!bCUTDOWN)
307 309
    {
308 310
        minx = miny = minz = 0;
309
        maxx = nbin[XX];
310
        maxy = nbin[YY];
311
        maxz = nbin[ZZ];
311
        maxx = nbin[XX] - 1;
312
        maxy = nbin[YY] - 1;
313
        maxz = nbin[ZZ] - 1;
312 314
    }
313 315

  
316
    boutput[0] = minx + std::max(iIGNOREOUTER, 0);
317
    boutput[1] = miny + std::max(iIGNOREOUTER, 0);
318
    boutput[2] = minz + std::max(iIGNOREOUTER, 0);
319
    eoutput[0] = maxx - std::max(iIGNOREOUTER, 0);
320
    eoutput[1] = maxy - std::max(iIGNOREOUTER, 0);
321
    eoutput[2] = maxz - std::max(iIGNOREOUTER, 0);
322

  
314 323
    /* OUTPUT */
315 324
    flp = gmx_ffopen("grid.cube", "w");
316 325
    fprintf(flp, "Spatial Distribution Function\n");
317 326
    fprintf(flp, "test\n");
318
    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);
319
    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxx-minx+1-(2*iIGNOREOUTER), rBINWIDTH*10./bohr, 0., 0.);
320
    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxy-miny+1-(2*iIGNOREOUTER), 0., rBINWIDTH*10./bohr, 0.);
321
    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", maxz-minz+1-(2*iIGNOREOUTER), 0., 0., rBINWIDTH*10./bohr);
327
    /*
328
       Values in .cube file represent the density at the grid point.
329
       Corresponding coordinates to the binned value is the center of the bin, i.e. + 0.5 to the bin index.
330
     */
331
    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", nidxp,
332
            (MINBIN[XX]+(boutput[XX] + 0.5)*rBINWIDTH)*10./bohr,
333
            (MINBIN[YY]+(boutput[YY] + 0.5)*rBINWIDTH)*10./bohr,
334
            (MINBIN[ZZ]+(boutput[ZZ] + 0.5)*rBINWIDTH)*10./bohr);
335
    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", eoutput[XX] - boutput[XX], rBINWIDTH*10./bohr, 0., 0.);
336
    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", eoutput[YY] - boutput[YY], 0., rBINWIDTH*10./bohr, 0.);
337
    fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", eoutput[ZZ] - boutput[ZZ], 0., 0., rBINWIDTH*10./bohr);
322 338
    for (i = 0; i < nidxp; i++)
323 339
    {
324 340
        v = 2;
......
376 392

  
377 393
    minval = 999;
378 394
    maxval = 0;
379
    for (k = 0; k < nbin[XX]; k++)
395
    for (k = boutput[XX]; k < eoutput[XX]; k++)
380 396
    {
381
        if (k < minx+iIGNOREOUTER || k > maxx-iIGNOREOUTER)
397
        for (j = boutput[YY]; j < eoutput[YY]; j++)
382 398
        {
383
            continue;
384
        }
385
        for (j = 0; j < nbin[YY]; j++)
386
        {
387
            if (j < miny+iIGNOREOUTER || j > maxy-iIGNOREOUTER)
399
            for (i = boutput[ZZ]; i < eoutput[ZZ]; i++)
388 400
            {
389
                continue;
390
            }
391
            for (i = 0; i < nbin[ZZ]; i++)
392
            {
393
                if (i < minz+iIGNOREOUTER || i > maxz-iIGNOREOUTER)
394
                {
395
                    continue;
396
                }
397 401
                tot += bin[k][j][i];
398 402
                if (bin[k][j][i] > maxval)
399 403
                {
......
407 411
        }
408 412
    }
409 413

  
410
    numcu = (maxx-minx+1-(2*iIGNOREOUTER))*(maxy-miny+1-(2*iIGNOREOUTER))*(maxz-minz+1-(2*iIGNOREOUTER));
414
    numcu = (eoutput[XX] - boutput[XX]) * (eoutput[YY] - boutput[YY]) * (eoutput[ZZ] - boutput[ZZ]);
411 415
    if (bCALCDIV)
412 416
    {
413
        norm = static_cast<double>(numcu*numfr)/tot;
417
        norm = static_cast<double>(static_cast<int64_t>(numcu)*numfr)/tot;
414 418
    }
415 419
    else
416 420
    {
417 421
        norm = 1.0;
418 422
    }
419 423

  
420
    for (k = 0; k < nbin[XX]; k++)
424
    for (k = boutput[XX]; k < eoutput[XX]; k++)
421 425
    {
422
        if (k < minx+iIGNOREOUTER || k > maxx-iIGNOREOUTER)
423
        {
424
            continue;
425
        }
426
        for (j = 0; j < nbin[YY]; j++)
426
        for (j = boutput[YY]; j < eoutput[YY]; j++)
427 427
        {
428
            if (j < miny+iIGNOREOUTER || j > maxy-iIGNOREOUTER)
429
            {
430
                continue;
431
            }
432
            for (i = 0; i < nbin[ZZ]; i++)
428
            for (i = boutput[ZZ]; i < eoutput[ZZ]; i++)
433 429
            {
434
                if (i < minz+iIGNOREOUTER || i > maxz-iIGNOREOUTER)
435
                {
436
                    continue;
437
                }
438 430
                fprintf(flp, "%12.6f ", static_cast<double>(norm*bin[k][j][i])/numfr);
439 431
            }
440 432
            fprintf(flp, "\n");
......
451 443
    else
452 444
    {
453 445
        printf("grid.cube contains counts per frame in all %d cubes\n", numcu);
454
        printf("Raw data: average %le, min %le, max %le\n", 1.0/norm, static_cast<double>(minval)/numfr, static_cast<double>(maxval)/numfr);
446
        printf("Raw data: average %le, min %le, max %le\n", static_cast<double>(tot)/numfr/numcu, static_cast<double>(minval)/numfr, static_cast<double>(maxval)/numfr);
455 447
    }
456 448

  
457 449
    return 0;
458
-