src/tools/gmx_densorder.c

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