src/gromacs/gmxana/gmx_densorder.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010-2018, The GROMACS development team.
   5  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include <cctype>
  39 #include <cmath>
  40 #include <cstring>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/correlationfunctions/autocorr.h"
  44 #include "gromacs/correlationfunctions/expfit.h"
  45 #include "gromacs/fileio/matio.h"
  46 #include "gromacs/fileio/tpxio.h"
  47 #include "gromacs/fileio/trxio.h"
  48 #include "gromacs/fileio/xvgr.h"
  49 #include "gromacs/gmxana/binsearch.h"
  50 #include "gromacs/gmxana/dens_filter.h"
  51 #include "gromacs/gmxana/gmx_ana.h"
  52 #include "gromacs/gmxana/gstat.h"
  53 #include "gromacs/gmxana/powerspect.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/pbcutil/rmpbc.h"
  58 #include "gromacs/topology/index.h"
  59 #include "gromacs/topology/topology.h"
  60 #include "gromacs/utility/arraysize.h"
  61 #include "gromacs/utility/binaryinformation.h"
  62 #include "gromacs/utility/cstringutil.h"
  63 #include "gromacs/utility/exceptions.h"
  64 #include "gromacs/utility/fatalerror.h"
  65 #include "gromacs/utility/futil.h"
  66 #include "gromacs/utility/smalloc.h"
  67
  68
  69 enum
  70 {
  71     methSEL,
  72     methBISECT,
  73     methFUNCFIT,
  74     methNR
  75 };
  76
  77 static void center_coords(const t_atoms* atoms, matrix box, rvec x0[], int axis)
  78 {
  79     int  i, m;
  80     real tmass, mm;
  81     rvec com, shift, box_center;
  82
  83     tmass = 0;
  84     clear_rvec(com);
  85     for (i = 0; (i < atoms->nr); i++)
  86     {
  87         mm = atoms->atom[i].m;
  88         tmass += mm;
  89         for (m = 0; (m < DIM); m++)
  90         {
  91             com[m] += mm * x0[i][m];
  92         }
  93     }
  94     for (m = 0; (m < DIM); m++)
  95     {
  96         com[m] /= tmass;
  97     }
  98     calc_box_center(ecenterDEF, box, box_center);
  99     rvec_sub(box_center, com, shift);
 100     shift[axis] -= box_center[axis];
 101
 102     for (i = 0; (i < atoms->nr); i++)
 103     {
 104         rvec_dec(x0[i], shift);
 105     }
 106 }
 107
 108
 109 static void density_in_time(const char*             fn,
 110                             int**                   index,
 111                             const int               gnx[],
 112                             real                    bw,
 113                             real                    bwz,
 114                             int                     nsttblock,
 115                             real*****               Densdevel,
 116                             int*                    xslices,
 117                             int*                    yslices,
 118                             int*                    zslices,
 119                             int*                    tblock,
 120                             const t_topology*       top,
 121                             PbcType                 pbcType,
 122                             int                     axis,
 123                             gmx_bool                bCenter,
 124                             gmx_bool                bps1d,
 125                             const gmx_output_env_t* oenv)
 126
 127 {
 128     /*
 129      * *****Densdevel pointer to array of density values in slices and frame-blocks
 130      * Densdevel[*nsttblock][*xslices][*yslices][*zslices] Densslice[x][y][z] nsttblock - nr of
 131      * frames in each time-block bw  widths of normal slices
 132      *
 133      * axis      - axis direction (normal to slices)
 134      * nndx - number ot atoms in **index
 135      * grpn      - group number in index
 136      */
 137     t_trxstatus* status;
 138     gmx_rmpbc_t  gpbc = nullptr;
 139     matrix       box;               /* Box - 3x3 -each step*/
 140     rvec*        x0;                /* List of Coord without PBC*/
 141     int          i, j,              /* loop indices, checks etc*/
 142             ax1 = 0, ax2 = 0,       /* tangent directions */
 143             framenr = 0,            /* frame number in trajectory*/
 144             slicex, slicey, slicez; /*slice # of x y z position */
 145     real*** Densslice = nullptr;    /* Density-slice in one frame*/
 146     real    dscale;                 /*physical scaling factor*/
 147     real    t, x, y, z;             /* time and coordinates*/
 148     rvec    bbww;
 149
 150     *tblock = 0; /* blocknr in block average - initialise to 0*/
 151     /* Axis: X=0, Y=1,Z=2 */
 152     switch (axis)
 153     {
 154         case 0:
 155             ax1 = YY;
 156             ax2 = ZZ; /*Surface: YZ*/
 157             break;
 158         case 1:
 159             ax1 = ZZ;
 160             ax2 = XX; /* Surface : XZ*/
 161             break;
 162         case 2:
 163             ax1 = XX;
 164             ax2 = YY; /* Surface XY*/
 165             break;
 166         default: gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 167     }
 168
 169     if (read_first_x(oenv, &status, fn, &t, &x0, box) == 0)
 170     {
 171         gmx_fatal(FARGS, "Could not read coordinates from file"); /* Open trajectory for read*/
 172     }
 173     *zslices = 1 + static_cast<int>(std::floor(box[axis][axis] / bwz));
 174     *yslices = 1 + static_cast<int>(std::floor(box[ax2][ax2] / bw));
 175     *xslices = 1 + static_cast<int>(std::floor(box[ax1][ax1] / bw));
 176     if (bps1d)
 177     {
 178         if (*xslices < *yslices)
 179         {
 180             *xslices = 1;
 181         }
 182         else
 183         {
 184             *yslices = 1;
 185         }
 186     }
 187     fprintf(stderr,
 188             "\nDividing the box in %5d x %5d x %5d slices with binw %f along axis %d\n",
 189             *xslices,
 190             *yslices,
 191             *zslices,
 192             bw,
 193             axis);
 194
 195
 196     /****Start trajectory processing***/
 197
 198     /*Initialize Densdevel and PBC-remove*/
 199     gpbc = gmx_rmpbc_init(&top->idef, pbcType, top->atoms.nr);
 200
 201     *Densdevel = nullptr;
 202
 203     do
 204     {
 205         bbww[XX] = box[ax1][ax1] / *xslices;
 206         bbww[YY] = box[ax2][ax2] / *yslices;
 207         bbww[ZZ] = box[axis][axis] / *zslices;
 208         gmx_rmpbc(gpbc, top->atoms.nr, box, x0);
 209         /*Reset Densslice every nsttblock steps*/
 210         /* The first conditional is for clang to understand that this branch is
 211          * always taken the first time. */
 212         if (Densslice == nullptr || framenr % nsttblock == 0)
 213         {
 214             snew(Densslice, *xslices);
 215             for (i = 0; i < *xslices; i++)
 216             {
 217                 snew(Densslice[i], *yslices);
 218                 for (j = 0; j < *yslices; j++)
 219                 {
 220                     snew(Densslice[i][j], *zslices);
 221                 }
 222             }
 223
 224             /* Allocate Memory to  extra frame in Densdevel -  rather stupid approach:
 225              * A single frame each time, although only every nsttblock steps.
 226              */
 227             srenew(*Densdevel, *tblock + 1);
 228             (*Densdevel)[*tblock] = Densslice;
 229         }
 230
 231         dscale = (*xslices) * (*yslices) * (*zslices) * gmx::c_amu
 232                  / (box[ax1][ax1] * box[ax2][ax2] * box[axis][axis] * nsttblock
 233                     * (gmx::c_nano * gmx::c_nano * gmx::c_nano));
 234
 235         if (bCenter)
 236         {
 237             center_coords(&top->atoms, box, x0, axis);
 238         }
 239
 240
 241         for (j = 0; j < gnx[0]; j++)
 242         { /*Loop over all atoms in selected index*/
 243             x = x0[index[0][j]][ax1];
 244             y = x0[index[0][j]][ax2];
 245             z = x0[index[0][j]][axis];
 246             while (x < 0)
 247             {
 248                 x += box[ax1][ax1];
 249             }
 250             while (x > box[ax1][ax1])
 251             {
 252                 x -= box[ax1][ax1];
 253             }
 254
 255             while (y < 0)
 256             {
 257                 y += box[ax2][ax2];
 258             }
 259             while (y > box[ax2][ax2])
 260             {
 261                 y -= box[ax2][ax2];
 262             }
 263
 264             while (z < 0)
 265             {
 266                 z += box[axis][axis];
 267             }
 268             while (z > box[axis][axis])
 269             {
 270                 z -= box[axis][axis];
 271             }
 272
 273             slicex = static_cast<int>(x / bbww[XX]) % *xslices;
 274             slicey = static_cast<int>(y / bbww[YY]) % *yslices;
 275             slicez = static_cast<int>(z / bbww[ZZ]) % *zslices;
 276             Densslice[slicex][slicey][slicez] += (top->atoms.atom[index[0][j]].m * dscale);
 277         }
 278
 279         framenr++;
 280
 281         if (framenr % nsttblock == 0)
 282         {
 283             /*Implicit incrementation of Densdevel via renewal of Densslice*/
 284             /*only every nsttblock steps*/
 285             (*tblock)++;
 286         }
 287
 288     } while (read_next_x(oenv, status, &t, x0, box));
 289
 290
 291     /*Free memory we no longer need and exit.*/
 292     gmx_rmpbc_done(gpbc);
 293     close_trx(status);
 294
 295     if (/* DISABLES CODE */ (false))
 296     {
 297         FILE* fp;
 298         fp = fopen("koko.xvg", "w");
 299         for (j = 0; (j < *zslices); j++)
 300         {
 301             fprintf(fp, "%5d", j);
 302             for (i = 0; (i < *tblock); i++)
 303             {
 304                 fprintf(fp, "  %10g", (*Densdevel)[i][9][1][j]);
 305             }
 306             fprintf(fp, "\n");
 307         }
 308         fclose(fp);
 309     }
 310 }
 311
 312 static void outputfield(const char* fldfn, real**** Densmap, int xslices, int yslices, int zslices, int tdim)
 313 {
 314     /*Debug-filename and filehandle*/
 315     FILE* fldH;
 316     int   n, i, j, k;
 317     int   dim[4];
 318     real  totdens = 0;
 319
 320     dim[0] = tdim;
 321     dim[1] = xslices;
 322     dim[2] = yslices;
 323     dim[3] = zslices;
 324
 325     fldH = gmx_ffopen(fldfn, "w");
 326     fwrite(dim, sizeof(int), 4, fldH);
 327     for (n = 0; n < tdim; n++)
 328     {
 329         for (i = 0; i < xslices; i++)
 330         {
 331             for (j = 0; j < yslices; j++)
 332             {
 333                 for (k = 0; k < zslices; k++)
 334                 {
 335                     fwrite(&(Densmap[n][i][j][k]), sizeof(real), 1, fldH);
 336                     totdens += (Densmap[n][i][j][k]);
 337                 }
 338             }
 339         }
 340     }
 341     totdens /= (xslices * yslices * zslices * tdim);
 342     fprintf(stderr, "Total density [kg/m^3]  %8f", totdens);
 343     gmx_ffclose(fldH);
 344 }
 345
 346 static void filterdensmap(real**** Densmap, int xslices, int yslices, int zslices, int tblocks, int ftsize)
 347 {
 348     real* kernel;
 349     real  std, var;
 350     int   i, j, n, order;
 351     order = ftsize / 2;
 352     std   = order / 2.0;
 353     var   = std * std;
 354     snew(kernel, ftsize);
 355     gausskernel(kernel, ftsize, var);
 356     for (n = 0; n < tblocks; n++)
 357     {
 358         for (i = 0; i < xslices; i++)
 359         {
 360             for (j = 0; j < yslices; j++)
 361             {
 362                 periodic_convolution(zslices, Densmap[n][i][j], ftsize, kernel);
 363             }
 364         }
 365     }
 366 }
 367
 368
 369 static void interfaces_txy(real****                Densmap,
 370                            int                     xslices,
 371                            int                     yslices,
 372                            int                     zslices,
 373                            int                     tblocks,
 374                            real                    binwidth,
 375                            int                     method,
 376                            real                    dens1,
 377                            real                    dens2,
 378                            t_interf****            intf1,
 379                            t_interf****            intf2,
 380                            const gmx_output_env_t* oenv)
 381 {
 382     /*Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)*/
 383     FILE*         xvg;
 384     real*         zDensavg; /* zDensavg[z]*/
 385     int           i, j, k, n;
 386     int           xysize;
 387     int           ndx1, ndx2, *zperm;
 388     real          densmid;
 389     real          splitpoint, startpoint, endpoint;
 390     real *        sigma1, *sigma2;
 391     double        beginfit1[4];
 392     double        beginfit2[4];
 393     double *      fit1 = nullptr, *fit2 = nullptr;
 394     const double* avgfit1;
 395     const double* avgfit2;
 396     const real    onehalf = 1.00 / 2.00;
 397     t_interf ***  int1    = nullptr,
 398              ***int2      = nullptr; /*Interface matrices [t][x,y] - last index in row-major order*/
 399     /*Create int1(t,xy) and int2(t,xy) arrays with correct number of interf_t elements*/
 400     xysize = xslices * yslices;
 401     snew(int1, tblocks);
 402     snew(int2, tblocks);
 403     for (i = 0; i < tblocks; i++)
 404     {
 405         snew(int1[i], xysize);
 406         snew(int2[i], xysize);
 407         for (j = 0; j < xysize; j++)
 408         {
 409             snew(int1[i][j], 1);
 410             snew(int2[i][j], 1);
 411             init_interf(int1[i][j]);
 412             init_interf(int2[i][j]);
 413         }
 414     }
 415
 416     if (method == methBISECT)
 417     {
 418         densmid = onehalf * (dens1 + dens2);
 419         snew(zperm, zslices);
 420         for (n = 0; n < tblocks; n++)
 421         {
 422             for (i = 0; i < xslices; i++)
 423             {
 424                 for (j = 0; j < yslices; j++)
 425                 {
 426                     rangeArray(zperm, zslices); /*reset permutation array to identity*/
 427                     /*Binsearch returns slice-nr where the order param is  <= setpoint sgmid*/
 428                     ndx1 = start_binsearch(Densmap[n][i][j], zperm, 0, zslices / 2 - 1, densmid, 1);
 429                     ndx2 = start_binsearch(Densmap[n][i][j], zperm, zslices / 2, zslices - 1, densmid, -1);
 430
 431                     /* Linear interpolation (for use later if time allows)
 432                      * rho_1s= Densmap[n][i][j][zperm[ndx1]]
 433                      * rho_1e =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
 434                      * rho_2s =Densmap[n][i][j][zperm[ndx2+1]]
 435                      * rho_2e =Densmap[n][i][j][zperm[ndx2]]
 436                      * For 1st interface we have:
 437                        densl= Densmap[n][i][j][zperm[ndx1]];
 438                        densr= Densmap[n][i][j][zperm[ndx1+1]];
 439                        alpha=(densmid-densl)/(densr-densl);
 440                        deltandx=zperm[ndx1+1]-zperm[ndx1];
 441
 442                        if(debug){
 443                        printf("Alpha, Deltandx  %f %i\n", alpha,deltandx);
 444                        }
 445                        if(abs(alpha)>1.0 || abs(deltandx)>3){
 446                        pos=zperm[ndx1];
 447                        spread=-1;
 448                        }
 449                        else {
 450                        pos=zperm[ndx1]+alpha*deltandx;
 451                        spread=binwidth*deltandx;
 452                        }
 453                      * For the 2nd interface  can use the same formulation, since alpha should become negative ie:
 454                      * alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/(Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
 455                      * deltandx=zperm[ndx2+1]-zperm[ndx2];
 456                      * pos=zperm[ndx2]+alpha*deltandx;   */
 457
 458                     /*After filtering we use the direct approach        */
 459                     int1[n][j + (i * yslices)]->Z = (zperm[ndx1] + onehalf) * binwidth;
 460                     int1[n][j + (i * yslices)]->t = binwidth;
 461                     int2[n][j + (i * yslices)]->Z = (zperm[ndx2] + onehalf) * binwidth;
 462                     int2[n][j + (i * yslices)]->t = binwidth;
 463                 }
 464             }
 465         }
 466     }
 467
 468     if (method == methFUNCFIT)
 469     {
 470         /*Assume a box divided in 2 along midpoint of z for starters*/
 471         startpoint = 0.0;
 472         endpoint   = binwidth * zslices;
 473         splitpoint = (startpoint + endpoint) / 2.0;
 474         /*Initial fit proposals*/
 475         beginfit1[0] = dens1;
 476         beginfit1[1] = dens2;
 477         beginfit1[2] = (splitpoint / 2);
 478         beginfit1[3] = 0.5;
 479
 480         beginfit2[0] = dens2;
 481         beginfit2[1] = dens1;
 482         beginfit2[2] = (3 * splitpoint / 2);
 483         beginfit2[3] = 0.5;
 484
 485         snew(zDensavg, zslices);
 486         snew(sigma1, zslices);
 487         snew(sigma2, zslices);
 488
 489         for (k = 0; k < zslices; k++)
 490         {
 491             sigma1[k] = sigma2[k] = 1;
 492         }
 493         /*Calculate average density along z - avoid smoothing by using coarse-grained-mesh*/
 494         for (k = 0; k < zslices; k++)
 495         {
 496             for (n = 0; n < tblocks; n++)
 497             {
 498                 for (i = 0; i < xslices; i++)
 499                 {
 500                     for (j = 0; j < yslices; j++)
 501                     {
 502                         zDensavg[k] += (Densmap[n][i][j][k] / (xslices * yslices * tblocks));
 503                     }
 504                 }
 505             }
 506         }
 507
 508         if (debug)
 509         {
 510             xvg = xvgropen(
 511                     "DensprofileonZ.xvg", "Averaged Densityprofile on Z", "z[nm]", "Density[kg/m^3]", oenv);
 512             for (k = 0; k < zslices; k++)
 513             {
 514                 fprintf(xvg, "%4f.3   %8f.4\n", k * binwidth, zDensavg[k]);
 515             }
 516             xvgrclose(xvg);
 517         }
 518
 519         /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
 520
 521         /*Fit 1st half of box*/
 522         do_lmfit(zslices, zDensavg, sigma1, binwidth, nullptr, startpoint, splitpoint, oenv, FALSE, effnERF, beginfit1, 8, nullptr);
 523         /*Fit 2nd half of box*/
 524         do_lmfit(zslices, zDensavg, sigma2, binwidth, nullptr, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 8, nullptr);
 525
 526         /*Initialise the const arrays for storing the average fit parameters*/
 527         avgfit1 = beginfit1;
 528         avgfit2 = beginfit2;
 529
 530
 531         /*Now do fit over each x  y and t slice to get Zint(x,y,t) - loop is very large, we potentially should average over time directly*/
 532         for (n = 0; n < tblocks; n++)
 533         {
 534             for (i = 0; i < xslices; i++)
 535             {
 536                 for (j = 0; j < yslices; j++)
 537                 {
 538                     /*Reinitialise fit for each mesh-point*/
 539                     srenew(fit1, 4);
 540                     srenew(fit2, 4);
 541                     for (k = 0; k < 4; k++)
 542                     {
 543                         fit1[k] = avgfit1[k];
 544                         fit2[k] = avgfit2[k];
 545                     }
 546                     /*Now fit and store in structures in row-major order int[n][i][j]*/
 547                     do_lmfit(zslices,
 548                              Densmap[n][i][j],
 549                              sigma1,
 550                              binwidth,
 551                              nullptr,
 552                              startpoint,
 553                              splitpoint,
 554                              oenv,
 555                              FALSE,
 556                              effnERF,
 557                              fit1,
 558                              0,
 559                              nullptr);
 560                     int1[n][j + (yslices * i)]->Z = fit1[2];
 561                     int1[n][j + (yslices * i)]->t = fit1[3];
 562                     do_lmfit(zslices,
 563                              Densmap[n][i][j],
 564                              sigma2,
 565                              binwidth,
 566                              nullptr,
 567                              splitpoint,
 568                              endpoint,
 569                              oenv,
 570                              FALSE,
 571                              effnERF,
 572                              fit2,
 573                              0,
 574                              nullptr);
 575                     int2[n][j + (yslices * i)]->Z = fit2[2];
 576                     int2[n][j + (yslices * i)]->t = fit2[3];
 577                 }
 578             }
 579         }
 580     }
 581
 582
 583     *intf1 = int1;
 584     *intf2 = int2;
 585 }
 586
 587 static void writesurftoxpms(t_interf***                      surf1,
 588                             t_interf***                      surf2,
 589                             int                              tblocks,
 590                             int                              xbins,
 591                             int                              ybins,
 592                             int                              zbins,
 593                             real                             bw,
 594                             real                             bwz,
 595                             gmx::ArrayRef<const std::string> outfiles,
 596                             int                              maplevels)
 597 {
 598     char   numbuf[STRLEN];
 599     int    n, i, j;
 600     real **profile1, **profile2;
 601     real   max1, max2, min1, min2, *xticks, *yticks;
 602     t_rgb  lo = { 0, 0, 0 };
 603     t_rgb  hi = { 1, 1, 1 };
 604     FILE * xpmfile1, *xpmfile2;
 605
 606     /*Prepare xpm structures for output*/
 607
 608     /*Allocate memory to tick's and matrices*/
 609     snew(xticks, xbins + 1);
 610     snew(yticks, ybins + 1);
 611
 612     profile1 = mk_matrix(xbins, ybins, FALSE);
 613     profile2 = mk_matrix(xbins, ybins, FALSE);
 614
 615     for (i = 0; i < xbins + 1; i++)
 616     {
 617         xticks[i] += bw;
 618     }
 619     for (j = 0; j < ybins + 1; j++)
 620     {
 621         yticks[j] += bw;
 622     }
 623
 624     xpmfile1 = gmx_ffopen(outfiles[0], "w");
 625     xpmfile2 = gmx_ffopen(outfiles[1], "w");
 626
 627     max1 = max2 = 0.0;
 628     min1 = min2 = zbins * bwz;
 629
 630     for (n = 0; n < tblocks; n++)
 631     {
 632         sprintf(numbuf, "tblock: %4i", n);
 633         /*Filling matrices for inclusion in xpm-files*/
 634         for (i = 0; i < xbins; i++)
 635         {
 636             for (j = 0; j < ybins; j++)
 637             {
 638                 profile1[i][j] = (surf1[n][j + ybins * i])->Z;
 639                 profile2[i][j] = (surf2[n][j + ybins * i])->Z;
 640                 /*Finding max and min values*/
 641                 if (profile1[i][j] > max1)
 642                 {
 643                     max1 = profile1[i][j];
 644                 }
 645                 if (profile1[i][j] < min1)
 646                 {
 647                     min1 = profile1[i][j];
 648                 }
 649                 if (profile2[i][j] > max2)
 650                 {
 651                     max2 = profile2[i][j];
 652                 }
 653                 if (profile2[i][j] < min2)
 654                 {
 655                     min2 = profile2[i][j];
 656                 }
 657             }
 658         }
 659
 660         write_xpm(xpmfile1, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile1, min1, max1, lo, hi, &maplevels);
 661         write_xpm(xpmfile2, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile2, min2, max2, lo, hi, &maplevels);
 662     }
 663
 664     gmx_ffclose(xpmfile1);
 665     gmx_ffclose(xpmfile2);
 666
 667
 668     sfree(profile1);
 669     sfree(profile2);
 670     sfree(xticks);
 671     sfree(yticks);
 672 }
 673
 674 static void writeraw(t_interf***                      int1,
 675                      t_interf***                      int2,
 676                      int                              tblocks,
 677                      int                              xbins,
 678                      int                              ybins,
 679                      gmx::ArrayRef<const std::string> fnms,
 680                      const gmx_output_env_t*          oenv)
 681 {
 682     FILE *raw1, *raw2;
 683     int   i, j, n;
 684
 685     raw1 = gmx_ffopen(fnms[0], "w");
 686     raw2 = gmx_ffopen(fnms[1], "w");
 687     try
 688     {
 689         gmx::BinaryInformationSettings settings;
 690         settings.generatedByHeader(true);
 691         settings.linePrefix("# ");
 692         gmx::printBinaryInformation(raw1, output_env_get_program_context(oenv), settings);
 693         gmx::printBinaryInformation(raw2, output_env_get_program_context(oenv), settings);
 694     }
 695     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 696     fprintf(raw1, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
 697     fprintf(raw2, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
 698     fprintf(raw1, "%i %i %i\n", tblocks, xbins, ybins);
 699     fprintf(raw2, "%i %i %i\n", tblocks, xbins, ybins);
 700     for (n = 0; n < tblocks; n++)
 701     {
 702         for (i = 0; i < xbins; i++)
 703         {
 704             for (j = 0; j < ybins; j++)
 705             {
 706                 fprintf(raw1,
 707                         "%i  %i  %8.5f  %6.4f\n",
 708                         i,
 709                         j,
 710                         (int1[n][j + ybins * i])->Z,
 711                         (int1[n][j + ybins * i])->t);
 712                 fprintf(raw2,
 713                         "%i  %i  %8.5f  %6.4f\n",
 714                         i,
 715                         j,
 716                         (int2[n][j + ybins * i])->Z,
 717                         (int2[n][j + ybins * i])->t);
 718             }
 719         }
 720     }
 721
 722     gmx_ffclose(raw1);
 723     gmx_ffclose(raw2);
 724 }
 725
 726
 727 int gmx_densorder(int argc, char* argv[])
 728 {
 729     static const char* desc[] = { "[THISMODULE] reduces a two-phase density distribution",
 730                                   "along an axis, computed over a MD trajectory,",
 731                                   "to 2D surfaces fluctuating in time, by a fit to",
 732                                   "a functional profile for interfacial densities.",
 733                                   "A time-averaged spatial representation of the",
 734                                   "interfaces can be output with the option [TT]-tavg[tt]." };
 735
 736     /* Extra arguments - but note how you always get the begin/end
 737      * options when running the program, without mentioning them here!
 738      */
 739
 740     gmx_output_env_t*  oenv;
 741     t_topology*        top;
 742     char**             grpname;
 743     PbcType            pbcType;
 744     int*               ngx;
 745     static real        binw      = 0.2;
 746     static real        binwz     = 0.05;
 747     static real        dens1     = 0.00;
 748     static real        dens2     = 1000.00;
 749     static int         ftorder   = 0;
 750     static int         nsttblock = 100;
 751     static int         axis      = 2;
 752     static const char* axtitle   = "Z";
 753     int**              index; /* Index list for single group*/
 754     int                xslices, yslices, zslices, tblock;
 755     static gmx_bool    bGraph   = FALSE;
 756     static gmx_bool    bCenter  = FALSE;
 757     static gmx_bool    bFourier = FALSE;
 758     static gmx_bool    bRawOut  = FALSE;
 759     static gmx_bool    bOut     = FALSE;
 760     static gmx_bool    b1d      = FALSE;
 761     static int         nlevels  = 100;
 762     /*Densitymap - Densmap[t][x][y][z]*/
 763     real**** Densmap = nullptr;
 764     /* Surfaces surf[t][surf_x,surf_y]*/
 765     t_interf ***surf1, ***surf2;
 766
 767     static const char* meth[] = { nullptr, "bisect", "functional", nullptr };
 768     int                eMeth;
 769
 770     t_pargs pa[] = {
 771         { "-1d", FALSE, etBOOL, { &b1d }, "Pseudo-1d interface geometry" },
 772         { "-bw",
 773           FALSE,
 774           etREAL,
 775           { &binw },
 776           "Binwidth of density distribution tangential to interface" },
 777         { "-bwn",
 778           FALSE,
 779           etREAL,
 780           { &binwz },
 781           "Binwidth of density distribution normal to interface" },
 782         { "-order",
 783           FALSE,
 784           etINT,
 785           { &ftorder },
 786           "Order of Gaussian filter, order 0 equates to NO filtering" },
 787         { "-axis", FALSE, etSTR, { &axtitle }, "Axis Direction - X, Y or Z" },
 788         { "-method", FALSE, etENUM, { meth }, "Interface location method" },
 789         { "-d1", FALSE, etREAL, { &dens1 }, "Bulk density phase 1 (at small z)" },
 790         { "-d2", FALSE, etREAL, { &dens2 }, "Bulk density phase 2 (at large z)" },
 791         { "-tblock", FALSE, etINT, { &nsttblock }, "Number of frames in one time-block average" },
 792         { "-nlevel", FALSE, etINT, { &nlevels }, "Number of Height levels in 2D - XPixMaps" }
 793     };
 794
 795
 796     t_filenm fnm[] = {
 797         { efTPR, "-s", nullptr, ffREAD }, /* this is for the topology */
 798         { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
 799         { efNDX, "-n", nullptr, ffREAD }, /* this is to select groups */
 800         { efDAT, "-o", "Density4D", ffOPTWR }, /* This is for outputting the entire 4D densityfield in binary format */
 801         { efOUT, "-or", nullptr, ffOPTWRMULT }, /* This is for writing out the entire information in the t_interf arrays */
 802         { efXPM, "-og", "interface", ffOPTWRMULT }, /* This is for writing out the interface meshes - one xpm-file per tblock*/
 803         { efOUT, "-Spect", "intfspect", ffOPTWRMULT }, /* This is for the trajectory averaged Fourier-spectra*/
 804     };
 805
 806 #define NFILE asize(fnm)
 807
 808     /* This is the routine responsible for adding default options,
 809      * calling the X/motif interface, etc. */
 810     if (!parse_common_args(
 811                 &argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 812     {
 813         return 0;
 814     }
 815
 816
 817     eMeth    = nenum(meth);
 818     bFourier = opt2bSet("-Spect", NFILE, fnm);
 819     bRawOut  = opt2bSet("-or", NFILE, fnm);
 820     bGraph   = opt2bSet("-og", NFILE, fnm);
 821     bOut     = opt2bSet("-o", NFILE, fnm);
 822     top      = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType);
 823     snew(grpname, 1);
 824     snew(index, 1);
 825     snew(ngx, 1);
 826
 827     /* Calculate axis */
 828     axis = toupper(axtitle[0]) - 'X';
 829
 830     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, ngx, index, grpname);
 831
 832     density_in_time(ftp2fn(efTRX, NFILE, fnm),
 833                     index,
 834                     ngx,
 835                     binw,
 836                     binwz,
 837                     nsttblock,
 838                     &Densmap,
 839                     &xslices,
 840                     &yslices,
 841                     &zslices,
 842                     &tblock,
 843                     top,
 844                     pbcType,
 845                     axis,
 846                     bCenter,
 847                     b1d,
 848                     oenv);
 849
 850     if (ftorder > 0)
 851     {
 852         filterdensmap(Densmap, xslices, yslices, zslices, tblock, 2 * ftorder + 1);
 853     }
 854
 855     if (bOut)
 856     {
 857         outputfield(opt2fn("-o", NFILE, fnm), Densmap, xslices, yslices, zslices, tblock);
 858     }
 859
 860     interfaces_txy(
 861             Densmap, xslices, yslices, zslices, tblock, binwz, eMeth, dens1, dens2, &surf1, &surf2, oenv);
 862
 863     if (bGraph)
 864     {
 865
 866         /*Output surface-xpms*/
 867         gmx::ArrayRef<const std::string> graphFiles = opt2fns("-og", NFILE, fnm);
 868         if (graphFiles.size() != 2)
 869         {
 870             gmx_fatal(FARGS, "No or not correct number (2) of output-files: %td", graphFiles.ssize());
 871         }
 872         writesurftoxpms(surf1, surf2, tblock, xslices, yslices, zslices, binw, binwz, graphFiles, zslices);
 873     }
 874
 875
 876     /*Output raw-data*/
 877     if (bRawOut)
 878     {
 879         gmx::ArrayRef<const std::string> rawFiles = opt2fns("-or", NFILE, fnm);
 880         if (rawFiles.size() != 2)
 881         {
 882             gmx_fatal(FARGS, "No or not correct number (2) of output-files: %td", rawFiles.ssize());
 883         }
 884         writeraw(surf1, surf2, tblock, xslices, yslices, rawFiles, oenv);
 885     }
 886
 887
 888     if (bFourier)
 889     {
 890         gmx::ArrayRef<const std::string> spectra = opt2fns("-Spect", NFILE, fnm);
 891         if (spectra.size() != 2)
 892         {
 893             gmx_fatal(FARGS, "No or not correct number (2) of output-file-series: %td", spectra.ssize());
 894         }
 895         powerspectavg_intf(surf1, surf2, tblock, xslices, yslices, spectra);
 896     }
 897
 898     sfree(Densmap);
 899     if (bGraph || bFourier || bRawOut)
 900     {
 901         sfree(surf1);
 902         sfree(surf2);
 903     }
 904
 905     return 0;
 906 }