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, 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) * AMU
 232                  / (box[ax1][ax1] * box[ax2][ax2] * box[axis][axis] * nsttblock * (NANO * NANO * NANO));
 233
 234         if (bCenter)
 235         {
 236             center_coords(&top->atoms, box, x0, axis);
 237         }
 238
 239
 240         for (j = 0; j < gnx[0]; j++)
 241         { /*Loop over all atoms in selected index*/
 242             x = x0[index[0][j]][ax1];
 243             y = x0[index[0][j]][ax2];
 244             z = x0[index[0][j]][axis];
 245             while (x < 0)
 246             {
 247                 x += box[ax1][ax1];
 248             }
 249             while (x > box[ax1][ax1])
 250             {
 251                 x -= box[ax1][ax1];
 252             }
 253
 254             while (y < 0)
 255             {
 256                 y += box[ax2][ax2];
 257             }
 258             while (y > box[ax2][ax2])
 259             {
 260                 y -= box[ax2][ax2];
 261             }
 262
 263             while (z < 0)
 264             {
 265                 z += box[axis][axis];
 266             }
 267             while (z > box[axis][axis])
 268             {
 269                 z -= box[axis][axis];
 270             }
 271
 272             slicex = static_cast<int>(x / bbww[XX]) % *xslices;
 273             slicey = static_cast<int>(y / bbww[YY]) % *yslices;
 274             slicez = static_cast<int>(z / bbww[ZZ]) % *zslices;
 275             Densslice[slicex][slicey][slicez] += (top->atoms.atom[index[0][j]].m * dscale);
 276         }
 277
 278         framenr++;
 279
 280         if (framenr % nsttblock == 0)
 281         {
 282             /*Implicit incrementation of Densdevel via renewal of Densslice*/
 283             /*only every nsttblock steps*/
 284             (*tblock)++;
 285         }
 286
 287     } while (read_next_x(oenv, status, &t, x0, box));
 288
 289
 290     /*Free memory we no longer need and exit.*/
 291     gmx_rmpbc_done(gpbc);
 292     close_trx(status);
 293
 294     if (/* DISABLES CODE */ (false))
 295     {
 296         FILE* fp;
 297         fp = fopen("koko.xvg", "w");
 298         for (j = 0; (j < *zslices); j++)
 299         {
 300             fprintf(fp, "%5d", j);
 301             for (i = 0; (i < *tblock); i++)
 302             {
 303                 fprintf(fp, "  %10g", (*Densdevel)[i][9][1][j]);
 304             }
 305             fprintf(fp, "\n");
 306         }
 307         fclose(fp);
 308     }
 309 }
 310
 311 static void outputfield(const char* fldfn, real**** Densmap, int xslices, int yslices, int zslices, int tdim)
 312 {
 313     /*Debug-filename and filehandle*/
 314     FILE* fldH;
 315     int   n, i, j, k;
 316     int   dim[4];
 317     real  totdens = 0;
 318
 319     dim[0] = tdim;
 320     dim[1] = xslices;
 321     dim[2] = yslices;
 322     dim[3] = zslices;
 323
 324     fldH = gmx_ffopen(fldfn, "w");
 325     fwrite(dim, sizeof(int), 4, fldH);
 326     for (n = 0; n < tdim; n++)
 327     {
 328         for (i = 0; i < xslices; i++)
 329         {
 330             for (j = 0; j < yslices; j++)
 331             {
 332                 for (k = 0; k < zslices; k++)
 333                 {
 334                     fwrite(&(Densmap[n][i][j][k]), sizeof(real), 1, fldH);
 335                     totdens += (Densmap[n][i][j][k]);
 336                 }
 337             }
 338         }
 339     }
 340     totdens /= (xslices * yslices * zslices * tdim);
 341     fprintf(stderr, "Total density [kg/m^3]  %8f", totdens);
 342     gmx_ffclose(fldH);
 343 }
 344
 345 static void filterdensmap(real**** Densmap, int xslices, int yslices, int zslices, int tblocks, int ftsize)
 346 {
 347     real* kernel;
 348     real  std, var;
 349     int   i, j, n, order;
 350     order = ftsize / 2;
 351     std   = order / 2.0;
 352     var   = std * std;
 353     snew(kernel, ftsize);
 354     gausskernel(kernel, ftsize, var);
 355     for (n = 0; n < tblocks; n++)
 356     {
 357         for (i = 0; i < xslices; i++)
 358         {
 359             for (j = 0; j < yslices; j++)
 360             {
 361                 periodic_convolution(zslices, Densmap[n][i][j], ftsize, kernel);
 362             }
 363         }
 364     }
 365 }
 366
 367
 368 static void interfaces_txy(real****                Densmap,
 369                            int                     xslices,
 370                            int                     yslices,
 371                            int                     zslices,
 372                            int                     tblocks,
 373                            real                    binwidth,
 374                            int                     method,
 375                            real                    dens1,
 376                            real                    dens2,
 377                            t_interf****            intf1,
 378                            t_interf****            intf2,
 379                            const gmx_output_env_t* oenv)
 380 {
 381     /*Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)*/
 382     FILE*         xvg;
 383     real*         zDensavg; /* zDensavg[z]*/
 384     int           i, j, k, n;
 385     int           xysize;
 386     int           ndx1, ndx2, *zperm;
 387     real          densmid;
 388     real          splitpoint, startpoint, endpoint;
 389     real *        sigma1, *sigma2;
 390     double        beginfit1[4];
 391     double        beginfit2[4];
 392     double *      fit1 = nullptr, *fit2 = nullptr;
 393     const double* avgfit1;
 394     const double* avgfit2;
 395     const real    onehalf = 1.00 / 2.00;
 396     t_interf ***  int1    = nullptr,
 397              ***int2      = nullptr; /*Interface matrices [t][x,y] - last index in row-major order*/
 398     /*Create int1(t,xy) and int2(t,xy) arrays with correct number of interf_t elements*/
 399     xysize = xslices * yslices;
 400     snew(int1, tblocks);
 401     snew(int2, tblocks);
 402     for (i = 0; i < tblocks; i++)
 403     {
 404         snew(int1[i], xysize);
 405         snew(int2[i], xysize);
 406         for (j = 0; j < xysize; j++)
 407         {
 408             snew(int1[i][j], 1);
 409             snew(int2[i][j], 1);
 410             init_interf(int1[i][j]);
 411             init_interf(int2[i][j]);
 412         }
 413     }
 414
 415     if (method == methBISECT)
 416     {
 417         densmid = onehalf * (dens1 + dens2);
 418         snew(zperm, zslices);
 419         for (n = 0; n < tblocks; n++)
 420         {
 421             for (i = 0; i < xslices; i++)
 422             {
 423                 for (j = 0; j < yslices; j++)
 424                 {
 425                     rangeArray(zperm, zslices); /*reset permutation array to identity*/
 426                     /*Binsearch returns slice-nr where the order param is  <= setpoint sgmid*/
 427                     ndx1 = start_binsearch(Densmap[n][i][j], zperm, 0, zslices / 2 - 1, densmid, 1);
 428                     ndx2 = start_binsearch(Densmap[n][i][j], zperm, zslices / 2, zslices - 1, densmid, -1);
 429
 430                     /* Linear interpolation (for use later if time allows)
 431                      * rho_1s= Densmap[n][i][j][zperm[ndx1]]
 432                      * rho_1e =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
 433                      * rho_2s =Densmap[n][i][j][zperm[ndx2+1]]
 434                      * rho_2e =Densmap[n][i][j][zperm[ndx2]]
 435                      * For 1st interface we have:
 436                        densl= Densmap[n][i][j][zperm[ndx1]];
 437                        densr= Densmap[n][i][j][zperm[ndx1+1]];
 438                        alpha=(densmid-densl)/(densr-densl);
 439                        deltandx=zperm[ndx1+1]-zperm[ndx1];
 440
 441                        if(debug){
 442                        printf("Alpha, Deltandx  %f %i\n", alpha,deltandx);
 443                        }
 444                        if(abs(alpha)>1.0 || abs(deltandx)>3){
 445                        pos=zperm[ndx1];
 446                        spread=-1;
 447                        }
 448                        else {
 449                        pos=zperm[ndx1]+alpha*deltandx;
 450                        spread=binwidth*deltandx;
 451                        }
 452                      * For the 2nd interface  can use the same formulation, since alpha should become negative ie:
 453                      * alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/(Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
 454                      * deltandx=zperm[ndx2+1]-zperm[ndx2];
 455                      * pos=zperm[ndx2]+alpha*deltandx;   */
 456
 457                     /*After filtering we use the direct approach        */
 458                     int1[n][j + (i * yslices)]->Z = (zperm[ndx1] + onehalf) * binwidth;
 459                     int1[n][j + (i * yslices)]->t = binwidth;
 460                     int2[n][j + (i * yslices)]->Z = (zperm[ndx2] + onehalf) * binwidth;
 461                     int2[n][j + (i * yslices)]->t = binwidth;
 462                 }
 463             }
 464         }
 465     }
 466
 467     if (method == methFUNCFIT)
 468     {
 469         /*Assume a box divided in 2 along midpoint of z for starters*/
 470         startpoint = 0.0;
 471         endpoint   = binwidth * zslices;
 472         splitpoint = (startpoint + endpoint) / 2.0;
 473         /*Initial fit proposals*/
 474         beginfit1[0] = dens1;
 475         beginfit1[1] = dens2;
 476         beginfit1[2] = (splitpoint / 2);
 477         beginfit1[3] = 0.5;
 478
 479         beginfit2[0] = dens2;
 480         beginfit2[1] = dens1;
 481         beginfit2[2] = (3 * splitpoint / 2);
 482         beginfit2[3] = 0.5;
 483
 484         snew(zDensavg, zslices);
 485         snew(sigma1, zslices);
 486         snew(sigma2, zslices);
 487
 488         for (k = 0; k < zslices; k++)
 489         {
 490             sigma1[k] = sigma2[k] = 1;
 491         }
 492         /*Calculate average density along z - avoid smoothing by using coarse-grained-mesh*/
 493         for (k = 0; k < zslices; k++)
 494         {
 495             for (n = 0; n < tblocks; n++)
 496             {
 497                 for (i = 0; i < xslices; i++)
 498                 {
 499                     for (j = 0; j < yslices; j++)
 500                     {
 501                         zDensavg[k] += (Densmap[n][i][j][k] / (xslices * yslices * tblocks));
 502                     }
 503                 }
 504             }
 505         }
 506
 507         if (debug)
 508         {
 509             xvg = xvgropen(
 510                     "DensprofileonZ.xvg", "Averaged Densityprofile on Z", "z[nm]", "Density[kg/m^3]", oenv);
 511             for (k = 0; k < zslices; k++)
 512             {
 513                 fprintf(xvg, "%4f.3   %8f.4\n", k * binwidth, zDensavg[k]);
 514             }
 515             xvgrclose(xvg);
 516         }
 517
 518         /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
 519
 520         /*Fit 1st half of box*/
 521         do_lmfit(zslices, zDensavg, sigma1, binwidth, nullptr, startpoint, splitpoint, oenv, FALSE, effnERF, beginfit1, 8, nullptr);
 522         /*Fit 2nd half of box*/
 523         do_lmfit(zslices, zDensavg, sigma2, binwidth, nullptr, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 8, nullptr);
 524
 525         /*Initialise the const arrays for storing the average fit parameters*/
 526         avgfit1 = beginfit1;
 527         avgfit2 = beginfit2;
 528
 529
 530         /*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*/
 531         for (n = 0; n < tblocks; n++)
 532         {
 533             for (i = 0; i < xslices; i++)
 534             {
 535                 for (j = 0; j < yslices; j++)
 536                 {
 537                     /*Reinitialise fit for each mesh-point*/
 538                     srenew(fit1, 4);
 539                     srenew(fit2, 4);
 540                     for (k = 0; k < 4; k++)
 541                     {
 542                         fit1[k] = avgfit1[k];
 543                         fit2[k] = avgfit2[k];
 544                     }
 545                     /*Now fit and store in structures in row-major order int[n][i][j]*/
 546                     do_lmfit(zslices,
 547                              Densmap[n][i][j],
 548                              sigma1,
 549                              binwidth,
 550                              nullptr,
 551                              startpoint,
 552                              splitpoint,
 553                              oenv,
 554                              FALSE,
 555                              effnERF,
 556                              fit1,
 557                              0,
 558                              nullptr);
 559                     int1[n][j + (yslices * i)]->Z = fit1[2];
 560                     int1[n][j + (yslices * i)]->t = fit1[3];
 561                     do_lmfit(zslices,
 562                              Densmap[n][i][j],
 563                              sigma2,
 564                              binwidth,
 565                              nullptr,
 566                              splitpoint,
 567                              endpoint,
 568                              oenv,
 569                              FALSE,
 570                              effnERF,
 571                              fit2,
 572                              0,
 573                              nullptr);
 574                     int2[n][j + (yslices * i)]->Z = fit2[2];
 575                     int2[n][j + (yslices * i)]->t = fit2[3];
 576                 }
 577             }
 578         }
 579     }
 580
 581
 582     *intf1 = int1;
 583     *intf2 = int2;
 584 }
 585
 586 static void writesurftoxpms(t_interf***                      surf1,
 587                             t_interf***                      surf2,
 588                             int                              tblocks,
 589                             int                              xbins,
 590                             int                              ybins,
 591                             int                              zbins,
 592                             real                             bw,
 593                             real                             bwz,
 594                             gmx::ArrayRef<const std::string> outfiles,
 595                             int                              maplevels)
 596 {
 597     char   numbuf[STRLEN];
 598     int    n, i, j;
 599     real **profile1, **profile2;
 600     real   max1, max2, min1, min2, *xticks, *yticks;
 601     t_rgb  lo = { 0, 0, 0 };
 602     t_rgb  hi = { 1, 1, 1 };
 603     FILE * xpmfile1, *xpmfile2;
 604
 605     /*Prepare xpm structures for output*/
 606
 607     /*Allocate memory to tick's and matrices*/
 608     snew(xticks, xbins + 1);
 609     snew(yticks, ybins + 1);
 610
 611     profile1 = mk_matrix(xbins, ybins, FALSE);
 612     profile2 = mk_matrix(xbins, ybins, FALSE);
 613
 614     for (i = 0; i < xbins + 1; i++)
 615     {
 616         xticks[i] += bw;
 617     }
 618     for (j = 0; j < ybins + 1; j++)
 619     {
 620         yticks[j] += bw;
 621     }
 622
 623     xpmfile1 = gmx_ffopen(outfiles[0], "w");
 624     xpmfile2 = gmx_ffopen(outfiles[1], "w");
 625
 626     max1 = max2 = 0.0;
 627     min1 = min2 = zbins * bwz;
 628
 629     for (n = 0; n < tblocks; n++)
 630     {
 631         sprintf(numbuf, "tblock: %4i", n);
 632         /*Filling matrices for inclusion in xpm-files*/
 633         for (i = 0; i < xbins; i++)
 634         {
 635             for (j = 0; j < ybins; j++)
 636             {
 637                 profile1[i][j] = (surf1[n][j + ybins * i])->Z;
 638                 profile2[i][j] = (surf2[n][j + ybins * i])->Z;
 639                 /*Finding max and min values*/
 640                 if (profile1[i][j] > max1)
 641                 {
 642                     max1 = profile1[i][j];
 643                 }
 644                 if (profile1[i][j] < min1)
 645                 {
 646                     min1 = profile1[i][j];
 647                 }
 648                 if (profile2[i][j] > max2)
 649                 {
 650                     max2 = profile2[i][j];
 651                 }
 652                 if (profile2[i][j] < min2)
 653                 {
 654                     min2 = profile2[i][j];
 655                 }
 656             }
 657         }
 658
 659         write_xpm(xpmfile1, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile1, min1, max1, lo, hi, &maplevels);
 660         write_xpm(xpmfile2, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile2, min2, max2, lo, hi, &maplevels);
 661     }
 662
 663     gmx_ffclose(xpmfile1);
 664     gmx_ffclose(xpmfile2);
 665
 666
 667     sfree(profile1);
 668     sfree(profile2);
 669     sfree(xticks);
 670     sfree(yticks);
 671 }
 672
 673 static void writeraw(t_interf***                      int1,
 674                      t_interf***                      int2,
 675                      int                              tblocks,
 676                      int                              xbins,
 677                      int                              ybins,
 678                      gmx::ArrayRef<const std::string> fnms,
 679                      const gmx_output_env_t*          oenv)
 680 {
 681     FILE *raw1, *raw2;
 682     int   i, j, n;
 683
 684     raw1 = gmx_ffopen(fnms[0], "w");
 685     raw2 = gmx_ffopen(fnms[1], "w");
 686     try
 687     {
 688         gmx::BinaryInformationSettings settings;
 689         settings.generatedByHeader(true);
 690         settings.linePrefix("# ");
 691         gmx::printBinaryInformation(raw1, output_env_get_program_context(oenv), settings);
 692         gmx::printBinaryInformation(raw2, output_env_get_program_context(oenv), settings);
 693     }
 694     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 695     fprintf(raw1, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
 696     fprintf(raw2, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
 697     fprintf(raw1, "%i %i %i\n", tblocks, xbins, ybins);
 698     fprintf(raw2, "%i %i %i\n", tblocks, xbins, ybins);
 699     for (n = 0; n < tblocks; n++)
 700     {
 701         for (i = 0; i < xbins; i++)
 702         {
 703             for (j = 0; j < ybins; j++)
 704             {
 705                 fprintf(raw1,
 706                         "%i  %i  %8.5f  %6.4f\n",
 707                         i,
 708                         j,
 709                         (int1[n][j + ybins * i])->Z,
 710                         (int1[n][j + ybins * i])->t);
 711                 fprintf(raw2,
 712                         "%i  %i  %8.5f  %6.4f\n",
 713                         i,
 714                         j,
 715                         (int2[n][j + ybins * i])->Z,
 716                         (int2[n][j + ybins * i])->t);
 717             }
 718         }
 719     }
 720
 721     gmx_ffclose(raw1);
 722     gmx_ffclose(raw2);
 723 }
 724
 725
 726 int gmx_densorder(int argc, char* argv[])
 727 {
 728     static const char* desc[] = { "[THISMODULE] reduces a two-phase density distribution",
 729                                   "along an axis, computed over a MD trajectory,",
 730                                   "to 2D surfaces fluctuating in time, by a fit to",
 731                                   "a functional profile for interfacial densities.",
 732                                   "A time-averaged spatial representation of the",
 733                                   "interfaces can be output with the option [TT]-tavg[tt]." };
 734
 735     /* Extra arguments - but note how you always get the begin/end
 736      * options when running the program, without mentioning them here!
 737      */
 738
 739     gmx_output_env_t*  oenv;
 740     t_topology*        top;
 741     char**             grpname;
 742     PbcType            pbcType;
 743     int*               ngx;
 744     static real        binw      = 0.2;
 745     static real        binwz     = 0.05;
 746     static real        dens1     = 0.00;
 747     static real        dens2     = 1000.00;
 748     static int         ftorder   = 0;
 749     static int         nsttblock = 100;
 750     static int         axis      = 2;
 751     static const char* axtitle   = "Z";
 752     int**              index; /* Index list for single group*/
 753     int                xslices, yslices, zslices, tblock;
 754     static gmx_bool    bGraph   = FALSE;
 755     static gmx_bool    bCenter  = FALSE;
 756     static gmx_bool    bFourier = FALSE;
 757     static gmx_bool    bRawOut  = FALSE;
 758     static gmx_bool    bOut     = FALSE;
 759     static gmx_bool    b1d      = FALSE;
 760     static int         nlevels  = 100;
 761     /*Densitymap - Densmap[t][x][y][z]*/
 762     real**** Densmap = nullptr;
 763     /* Surfaces surf[t][surf_x,surf_y]*/
 764     t_interf ***surf1, ***surf2;
 765
 766     static const char* meth[] = { nullptr, "bisect", "functional", nullptr };
 767     int                eMeth;
 768
 769     t_pargs pa[] = {
 770         { "-1d", FALSE, etBOOL, { &b1d }, "Pseudo-1d interface geometry" },
 771         { "-bw",
 772           FALSE,
 773           etREAL,
 774           { &binw },
 775           "Binwidth of density distribution tangential to interface" },
 776         { "-bwn",
 777           FALSE,
 778           etREAL,
 779           { &binwz },
 780           "Binwidth of density distribution normal to interface" },
 781         { "-order",
 782           FALSE,
 783           etINT,
 784           { &ftorder },
 785           "Order of Gaussian filter, order 0 equates to NO filtering" },
 786         { "-axis", FALSE, etSTR, { &axtitle }, "Axis Direction - X, Y or Z" },
 787         { "-method", FALSE, etENUM, { meth }, "Interface location method" },
 788         { "-d1", FALSE, etREAL, { &dens1 }, "Bulk density phase 1 (at small z)" },
 789         { "-d2", FALSE, etREAL, { &dens2 }, "Bulk density phase 2 (at large z)" },
 790         { "-tblock", FALSE, etINT, { &nsttblock }, "Number of frames in one time-block average" },
 791         { "-nlevel", FALSE, etINT, { &nlevels }, "Number of Height levels in 2D - XPixMaps" }
 792     };
 793
 794
 795     t_filenm fnm[] = {
 796         { efTPR, "-s", nullptr, ffREAD }, /* this is for the topology */
 797         { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
 798         { efNDX, "-n", nullptr, ffREAD }, /* this is to select groups */
 799         { efDAT, "-o", "Density4D", ffOPTWR }, /* This is for outputting the entire 4D densityfield in binary format */
 800         { efOUT, "-or", nullptr, ffOPTWRMULT }, /* This is for writing out the entire information in the t_interf arrays */
 801         { efXPM, "-og", "interface", ffOPTWRMULT }, /* This is for writing out the interface meshes - one xpm-file per tblock*/
 802         { efOUT, "-Spect", "intfspect", ffOPTWRMULT }, /* This is for the trajectory averaged Fourier-spectra*/
 803     };
 804
 805 #define NFILE asize(fnm)
 806
 807     /* This is the routine responsible for adding default options,
 808      * calling the X/motif interface, etc. */
 809     if (!parse_common_args(
 810                 &argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 811     {
 812         return 0;
 813     }
 814
 815
 816     eMeth    = nenum(meth);
 817     bFourier = opt2bSet("-Spect", NFILE, fnm);
 818     bRawOut  = opt2bSet("-or", NFILE, fnm);
 819     bGraph   = opt2bSet("-og", NFILE, fnm);
 820     bOut     = opt2bSet("-o", NFILE, fnm);
 821     top      = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType);
 822     snew(grpname, 1);
 823     snew(index, 1);
 824     snew(ngx, 1);
 825
 826     /* Calculate axis */
 827     axis = toupper(axtitle[0]) - 'X';
 828
 829     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, ngx, index, grpname);
 830
 831     density_in_time(ftp2fn(efTRX, NFILE, fnm),
 832                     index,
 833                     ngx,
 834                     binw,
 835                     binwz,
 836                     nsttblock,
 837                     &Densmap,
 838                     &xslices,
 839                     &yslices,
 840                     &zslices,
 841                     &tblock,
 842                     top,
 843                     pbcType,
 844                     axis,
 845                     bCenter,
 846                     b1d,
 847                     oenv);
 848
 849     if (ftorder > 0)
 850     {
 851         filterdensmap(Densmap, xslices, yslices, zslices, tblock, 2 * ftorder + 1);
 852     }
 853
 854     if (bOut)
 855     {
 856         outputfield(opt2fn("-o", NFILE, fnm), Densmap, xslices, yslices, zslices, tblock);
 857     }
 858
 859     interfaces_txy(
 860             Densmap, xslices, yslices, zslices, tblock, binwz, eMeth, dens1, dens2, &surf1, &surf2, oenv);
 861
 862     if (bGraph)
 863     {
 864
 865         /*Output surface-xpms*/
 866         gmx::ArrayRef<const std::string> graphFiles = opt2fns("-og", NFILE, fnm);
 867         if (graphFiles.size() != 2)
 868         {
 869             gmx_fatal(FARGS, "No or not correct number (2) of output-files: %td", graphFiles.ssize());
 870         }
 871         writesurftoxpms(surf1, surf2, tblock, xslices, yslices, zslices, binw, binwz, graphFiles, zslices);
 872     }
 873
 874
 875     /*Output raw-data*/
 876     if (bRawOut)
 877     {
 878         gmx::ArrayRef<const std::string> rawFiles = opt2fns("-or", NFILE, fnm);
 879         if (rawFiles.size() != 2)
 880         {
 881             gmx_fatal(FARGS, "No or not correct number (2) of output-files: %td", rawFiles.ssize());
 882         }
 883         writeraw(surf1, surf2, tblock, xslices, yslices, rawFiles, oenv);
 884     }
 885
 886
 887     if (bFourier)
 888     {
 889         gmx::ArrayRef<const std::string> spectra = opt2fns("-Spect", NFILE, fnm);
 890         if (spectra.size() != 2)
 891         {
 892             gmx_fatal(FARGS, "No or not correct number (2) of output-file-series: %td", spectra.ssize());
 893         }
 894         powerspectavg_intf(surf1, surf2, tblock, xslices, yslices, spectra);
 895     }
 896
 897     sfree(Densmap);
 898     if (bGraph || bFourier || bRawOut)
 899     {
 900         sfree(surf1);
 901         sfree(surf2);
 902     }
 903
 904     return 0;
 905 }