src/gromacs/gmxana/gmx_densorder.cpp

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