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