src/gromacs/gmxana/gmx_hydorder.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
  36 #include "gmxpre.h"
  37
  38 #include <cmath>
  39 #include <cstring>
  40
  41 #include "gromacs/commandline/pargs.h"
  42 #include "gromacs/fileio/confio.h"
  43 #include "gromacs/fileio/matio.h"
  44 #include "gromacs/fileio/trxio.h"
  45 #include "gromacs/fileio/xvgr.h"
  46 #include "gromacs/gmxana/binsearch.h"
  47 #include "gromacs/gmxana/gmx_ana.h"
  48 #include "gromacs/gmxana/gstat.h"
  49 #include "gromacs/gmxana/powerspect.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/pbcutil/pbc.h"
  52 #include "gromacs/pbcutil/rmpbc.h"
  53 #include "gromacs/topology/index.h"
  54 #include "gromacs/topology/topology.h"
  55 #include "gromacs/utility/arraysize.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/gmxassert.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 static void find_tetra_order_grid(t_topology top, int ePBC,
  63                                   int natoms, matrix box,
  64                                   rvec x[], int maxidx, int index[],
  65                                   real *sgmean, real *skmean,
  66                                   int nslicex, int nslicey, int nslicez,
  67                                   real ***sggrid, real ***skgrid)
  68 {
  69     int         ix, jx, i, j, k, l, n, *nn[4];
  70     rvec        dx, rj, rk, urk, urj;
  71     real        cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
  72     t_pbc       pbc;
  73     int         slindex_x, slindex_y, slindex_z;
  74     int      ***sl_count;
  75     real        onethird = 1.0/3.0;
  76     gmx_rmpbc_t gpbc;
  77
  78     /*  dmat = init_mat(maxidx, FALSE); */
  79
  80     box2 = box[XX][XX] * box[XX][XX];
  81
  82     /* Initialize expanded sl_count array */
  83     snew(sl_count, nslicex);
  84     for (i = 0; i < nslicex; i++)
  85     {
  86         snew(sl_count[i], nslicey);
  87         for (j = 0; j < nslicey; j++)
  88         {
  89             snew(sl_count[i][j], nslicez);
  90         }
  91     }
  92
  93
  94     for (i = 0; (i < 4); i++)
  95     {
  96         snew(r_nn[i], natoms);
  97         snew(nn[i], natoms);
  98
  99         for (j = 0; (j < natoms); j++)
 100         {
 101             r_nn[i][j] = box2;
 102         }
 103     }
 104
 105     snew(sgmol, maxidx);
 106     snew(skmol, maxidx);
 107
 108     /* Must init pbc every step because of pressure coupling */
 109     set_pbc(&pbc, ePBC, box);
 110     gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 111     gmx_rmpbc(gpbc, natoms, box, x);
 112
 113     *sgmean = 0.0;
 114     *skmean = 0.0;
 115     l       = 0;
 116     for (i = 0; (i < maxidx); i++)
 117     {   /* loop over index file */
 118         ix = index[i];
 119         for (j = 0; (j < maxidx); j++)
 120         {
 121
 122             if (i == j)
 123             {
 124                 continue;
 125             }
 126
 127             jx = index[j];
 128
 129             pbc_dx(&pbc, x[ix], x[jx], dx);
 130             r2 = iprod(dx, dx);
 131
 132             /* set_mat_entry(dmat,i,j,r2); */
 133
 134             /* determine the nearest neighbours */
 135             if (r2 < r_nn[0][i])
 136             {
 137                 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
 138                 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
 139                 r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
 140                 r_nn[0][i] = r2;         nn[0][i] = j;
 141             }
 142             else if (r2 < r_nn[1][i])
 143             {
 144                 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
 145                 r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
 146                 r_nn[1][i] = r2;         nn[1][i] = j;
 147             }
 148             else if (r2 < r_nn[2][i])
 149             {
 150                 r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
 151                 r_nn[2][i] = r2;         nn[2][i] = j;
 152             }
 153             else if (r2 < r_nn[3][i])
 154             {
 155                 r_nn[3][i] = r2;         nn[3][i] = j;
 156             }
 157         }
 158
 159
 160         /* calculate mean distance between nearest neighbours */
 161         rmean = 0;
 162         for (j = 0; (j < 4); j++)
 163         {
 164             r_nn[j][i] = std::sqrt(r_nn[j][i]);
 165             rmean     += r_nn[j][i];
 166         }
 167         rmean /= 4;
 168
 169         n        = 0;
 170         sgmol[i] = 0.0;
 171         skmol[i] = 0.0;
 172
 173         /* Chau1998a eqn 3 */
 174         /* angular part tetrahedrality order parameter per atom */
 175         for (j = 0; (j < 3); j++)
 176         {
 177             for (k = j+1; (k < 4); k++)
 178             {
 179                 pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
 180                 pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
 181
 182                 unitv(rk, urk);
 183                 unitv(rj, urj);
 184
 185                 cost  = iprod(urk, urj) + onethird;
 186                 cost2 = cost * cost;
 187
 188                 sgmol[i] += cost2;
 189                 l++;
 190                 n++;
 191             }
 192         }
 193         /* normalize sgmol between 0.0 and 1.0 */
 194         sgmol[i] = 3*sgmol[i]/32;
 195         *sgmean += sgmol[i];
 196
 197         /* distance part tetrahedrality order parameter per atom */
 198         rmean2 = 4 * 3 * rmean * rmean;
 199         for (j = 0; (j < 4); j++)
 200         {
 201             skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
 202             /*      printf("%d %f (%f %f %f %f) \n",
 203                     i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
 204              */
 205         }
 206
 207         *skmean += skmol[i];
 208
 209         /* Compute sliced stuff in x y z*/
 210         slindex_x = static_cast<int>(std::round((1+x[i][XX]/box[XX][XX])*nslicex)) % nslicex;
 211         slindex_y = static_cast<int>(std::round((1+x[i][YY]/box[YY][YY])*nslicey)) % nslicey;
 212         slindex_z = static_cast<int>(std::round((1+x[i][ZZ]/box[ZZ][ZZ])*nslicez)) % nslicez;
 213         sggrid[slindex_x][slindex_y][slindex_z] += sgmol[i];
 214         skgrid[slindex_x][slindex_y][slindex_z] += skmol[i];
 215         (sl_count[slindex_x][slindex_y][slindex_z])++;
 216     } /* loop over entries in index file */
 217
 218     *sgmean /= maxidx;
 219     *skmean /= maxidx;
 220
 221     for (i = 0; (i < nslicex); i++)
 222     {
 223         for (j = 0; j < nslicey; j++)
 224         {
 225             for (k = 0; k < nslicez; k++)
 226             {
 227                 if (sl_count[i][j][k] > 0)
 228                 {
 229                     sggrid[i][j][k] /= sl_count[i][j][k];
 230                     skgrid[i][j][k] /= sl_count[i][j][k];
 231                 }
 232             }
 233         }
 234     }
 235
 236     sfree(sl_count);
 237     sfree(sgmol);
 238     sfree(skmol);
 239     for (i = 0; (i < 4); i++)
 240     {
 241         sfree(r_nn[i]);
 242         sfree(nn[i]);
 243     }
 244 }
 245
 246 /*Determines interface from tetrahedral order parameter in box with specified binwidth.  */
 247 /*Outputs interface positions(bins), the number of timeframes, and the number of surface-mesh points in xy*/
 248
 249 static void calc_tetra_order_interface(const char *fnNDX, const char *fnTPS, const char *fnTRX, real binw, int tblock,
 250                                        int *nframes,  int *nslicex, int *nslicey,
 251                                        real sgang1, real sgang2, real ****intfpos,
 252                                        gmx_output_env_t *oenv)
 253 {
 254     FILE         *fpsg   = NULL, *fpsk = NULL;
 255     t_topology    top;
 256     int           ePBC;
 257     t_trxstatus  *status;
 258     int           natoms;
 259     real          t;
 260     rvec         *xtop, *x;
 261     matrix        box;
 262     real          sg, sk, sgintf;
 263     int         **index   = NULL;
 264     char        **grpname = NULL;
 265     int           i, j, k, n, *isize, ng, nslicez, framenr;
 266     real       ***sg_grid = NULL, ***sk_grid = NULL, ***sg_fravg = NULL, ***sk_fravg = NULL, ****sk_4d = NULL, ****sg_4d = NULL;
 267     int          *perm;
 268     int           ndx1, ndx2;
 269     int           bins;
 270     const real    onehalf = 1.0/2.0;
 271     /* real   ***intfpos[2]; pointers to arrays of two interface positions zcoord(framenr,xbin,ybin): intfpos[interface_index][t][nslicey*x+y]
 272      * i.e 1D Row-major order in (t,x,y) */
 273
 274
 275     read_tps_conf(fnTPS, &top, &ePBC, &xtop, NULL, box, FALSE);
 276
 277     *nslicex = static_cast<int>(box[XX][XX]/binw + onehalf); /*Calculate slicenr from binwidth*/
 278     *nslicey = static_cast<int>(box[YY][YY]/binw + onehalf);
 279     nslicez  = static_cast<int>(box[ZZ][ZZ]/binw +  onehalf);
 280
 281
 282
 283     ng = 1;
 284     /* get index groups */
 285     printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
 286     snew(grpname, ng);
 287     snew(index, ng);
 288     snew(isize, ng);
 289     get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
 290
 291     /* Analyze trajectory */
 292     natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
 293     if (natoms > top.atoms.nr)
 294     {
 295         gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
 296                   top.atoms.nr, natoms);
 297     }
 298     check_index(NULL, ng, index[0], NULL, natoms);
 299
 300
 301     /*Prepare structures for temporary storage of frame info*/
 302     snew(sg_grid, *nslicex);
 303     snew(sk_grid, *nslicex);
 304     for (i = 0; i < *nslicex; i++)
 305     {
 306         snew(sg_grid[i], *nslicey);
 307         snew(sk_grid[i], *nslicey);
 308         for (j = 0; j < *nslicey; j++)
 309         {
 310             snew(sg_grid[i][j], nslicez);
 311             snew(sk_grid[i][j], nslicez);
 312         }
 313     }
 314
 315     sg_4d    = NULL;
 316     sk_4d    = NULL;
 317     *nframes = 0;
 318     framenr  = 0;
 319
 320 /* Loop over frames*/
 321     do
 322     {
 323         /*Initialize box meshes (temporary storage for each tblock frame -reinitialise every tblock steps */
 324         if (framenr%tblock == 0)
 325         {
 326             srenew(sk_4d, *nframes+1);
 327             srenew(sg_4d, *nframes+1);
 328             snew(sg_fravg, *nslicex);
 329             snew(sk_fravg, *nslicex);
 330             for (i = 0; i < *nslicex; i++)
 331             {
 332                 snew(sg_fravg[i], *nslicey);
 333                 snew(sk_fravg[i], *nslicey);
 334                 for (j = 0; j < *nslicey; j++)
 335                 {
 336                     snew(sg_fravg[i][j], nslicez);
 337                     snew(sk_fravg[i][j], nslicez);
 338                 }
 339             }
 340         }
 341
 342         find_tetra_order_grid(top, ePBC, natoms, box, x, isize[0], index[0],
 343                               &sg, &sk, *nslicex, *nslicey, nslicez, sg_grid, sk_grid);
 344         GMX_RELEASE_ASSERT(sk_fravg != NULL, "Trying to dereference NULL sk_fravg pointer");
 345         for (i = 0; i < *nslicex; i++)
 346         {
 347             for (j = 0; j < *nslicey; j++)
 348             {
 349                 for (k = 0; k < nslicez; k++)
 350                 {
 351                     sk_fravg[i][j][k] += sk_grid[i][j][k]/tblock;
 352                     sg_fravg[i][j][k] += sg_grid[i][j][k]/tblock;
 353                 }
 354             }
 355         }
 356
 357         framenr++;
 358
 359         if (framenr%tblock == 0)
 360         {
 361             GMX_RELEASE_ASSERT(sk_4d != NULL, "Trying to dereference NULL sk_4d pointer");
 362             sk_4d[*nframes] = sk_fravg;
 363             sg_4d[*nframes] = sg_fravg;
 364             (*nframes)++;
 365         }
 366
 367     }
 368     while (read_next_x(oenv, status, &t, x, box));
 369     close_trj(status);
 370
 371     sfree(grpname);
 372     sfree(index);
 373     sfree(isize);
 374
 375     /*Debugging for printing out the entire order parameter meshes.*/
 376     if (debug)
 377     {
 378         fpsg = xvgropen("sg_ang_mesh", "S\\sg\\N Angle Order Parameter / Meshpoint", "(nm)", "S\\sg\\N", oenv);
 379         fpsk = xvgropen("sk_dist_mesh", "S\\sk\\N Distance Order Parameter / Meshpoint", "(nm)", "S\\sk\\N", oenv);
 380         for (n = 0; n < (*nframes); n++)
 381         {
 382             fprintf(fpsg, "%i\n", n);
 383             fprintf(fpsk, "%i\n", n);
 384             for (i = 0; (i < *nslicex); i++)
 385             {
 386                 for (j = 0; j < *nslicey; j++)
 387                 {
 388                     for (k = 0; k < nslicez; k++)
 389                     {
 390                         fprintf(fpsg, "%4f  %4f  %4f  %8f\n", (i+0.5)*box[XX][XX]/(*nslicex), (j+0.5)*box[YY][YY]/(*nslicey), (k+0.5)*box[ZZ][ZZ]/nslicez, sg_4d[n][i][j][k]);
 391                         fprintf(fpsk, "%4f  %4f  %4f  %8f\n", (i+0.5)*box[XX][XX]/(*nslicex), (j+0.5)*box[YY][YY]/(*nslicey), (k+0.5)*box[ZZ][ZZ]/nslicez, sk_4d[n][i][j][k]);
 392                     }
 393                 }
 394             }
 395         }
 396         xvgrclose(fpsg);
 397         xvgrclose(fpsk);
 398     }
 399
 400
 401     /* Find positions of interface z by scanning orderparam for each frame and for each xy-mesh cylinder along z*/
 402
 403     /*Simple trial: assume interface is in the middle of -sgang1 and sgang2*/
 404     sgintf = 0.5*(sgang1+sgang2);
 405
 406
 407     /*Allocate memory for interface arrays; */
 408     snew((*intfpos), 2);
 409     snew((*intfpos)[0], *nframes);
 410     snew((*intfpos)[1], *nframes);
 411
 412     bins = (*nslicex)*(*nslicey);
 413
 414
 415     snew(perm, nslicez);  /*permutation array for sorting along normal coordinate*/
 416
 417
 418     for (n = 0; n < *nframes; n++)
 419     {
 420         snew((*intfpos)[0][n], bins);
 421         snew((*intfpos)[1][n], bins);
 422         for (i = 0; i < *nslicex; i++)
 423         {
 424             for (j = 0; j < *nslicey; j++)
 425             {
 426                 rangeArray(perm, nslicez); /*reset permutation array to identity*/
 427                 /*Binsearch returns 2 bin-numbers where the order param is  <= setpoint sgintf*/
 428                 ndx1 = start_binsearch(sg_4d[n][i][j], perm, 0, nslicez/2-1, sgintf, 1);
 429                 ndx2 = start_binsearch(sg_4d[n][i][j], perm, nslicez/2, nslicez-1, sgintf, -1);
 430                 /*Use linear interpolation to smooth out the interface position*/
 431
 432                 /*left interface (0)*/
 433                 /*if((sg_4d[n][i][j][perm[ndx1+1]]-sg_4d[n][i][j][perm[ndx1]])/sg_4d[n][i][j][perm[ndx1]] > 0.01){
 434                    pos=( (sgintf-sg_4d[n][i][j][perm[ndx1]])*perm[ndx1+1]+(sg_4d[n][i][j][perm[ndx1+1]]-sgintf)*perm[ndx1 ])*/
 435                 (*intfpos)[0][n][j+*nslicey*i] = (perm[ndx1]+onehalf)*binw;
 436                 /*right interface (1)*/
 437                 /*alpha=(sgintf-sg_4d[n][i][j][perm[ndx2]])/(sg_4d[n][i][j][perm[ndx2]+1]-sg_4d[n][i][j][perm[ndx2]]);*/
 438                 /*(*intfpos)[1][n][j+*nslicey*i]=((1-alpha)*perm[ndx2]+alpha*(perm[ndx2]+1)+onehalf)*box[ZZ][ZZ]/nslicez;*/
 439                 (*intfpos)[1][n][j+*nslicey*i] = (perm[ndx2]+onehalf)*binw;
 440             }
 441         }
 442     }
 443
 444
 445     sfree(sk_4d);
 446     sfree(sg_4d);
 447 }
 448
 449 static void writesurftoxpms(real ***surf, int tblocks, int xbins, int ybins, real bw, char **outfiles, int maplevels )
 450 {
 451
 452     char   numbuf[8];
 453     int    n, i, j;
 454     real **profile1, **profile2;
 455     real   max1, max2, min1, min2, *xticks, *yticks;
 456     t_rgb  lo = {1, 1, 1};
 457     t_rgb  hi = {0, 0, 0};
 458     FILE  *xpmfile1, *xpmfile2;
 459
 460 /*Prepare xpm structures for output*/
 461
 462 /*Allocate memory to tick's and matrices*/
 463     snew (xticks, xbins+1);
 464     snew (yticks, ybins+1);
 465
 466     profile1 = mk_matrix(xbins, ybins, FALSE);
 467     profile2 = mk_matrix(xbins, ybins, FALSE);
 468
 469     for (i = 0; i < xbins+1; i++)
 470     {
 471         xticks[i] += bw;
 472     }
 473     for (j = 0; j < ybins+1; j++)
 474     {
 475         yticks[j] += bw;
 476     }
 477
 478     xpmfile1 = gmx_ffopen(outfiles[0], "w");
 479     xpmfile2 = gmx_ffopen(outfiles[1], "w");
 480
 481     max1 = max2 = 0.0;
 482     min1 = min2 = 1000.00;
 483
 484     for (n = 0; n < tblocks; n++)
 485     {
 486         sprintf(numbuf, "%5d", n);
 487         /*Filling matrices for inclusion in xpm-files*/
 488         for (i = 0; i < xbins; i++)
 489         {
 490             for (j = 0; j < ybins; j++)
 491             {
 492                 profile1[i][j] = (surf[0][n][j+ybins*i]);
 493                 profile2[i][j] = (surf[1][n][j+ybins*i]);
 494                 /*Finding max and min values*/
 495                 if (profile1[i][j] > max1)
 496                 {
 497                     max1 = profile1[i][j];
 498                 }
 499                 if (profile1[i][j] < min1)
 500                 {
 501                     min1 = profile1[i][j];
 502                 }
 503                 if (profile2[i][j] > max2)
 504                 {
 505                     max2 = profile2[i][j];
 506                 }
 507                 if (profile2[i][j] < min2)
 508                 {
 509                     min2 = profile2[i][j];
 510                 }
 511             }
 512         }
 513
 514         write_xpm(xpmfile1, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile1, min1, max1, lo, hi, &maplevels);
 515         write_xpm(xpmfile2, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile2, min2, max2, lo, hi, &maplevels);
 516     }
 517
 518     gmx_ffclose(xpmfile1);
 519     gmx_ffclose(xpmfile2);
 520
 521
 522
 523     sfree(profile1);
 524     sfree(profile2);
 525     sfree(xticks);
 526     sfree(yticks);
 527 }
 528
 529
 530 static void writeraw(real ***surf, int tblocks, int xbins, int ybins, char **fnms)
 531 {
 532     FILE *raw1, *raw2;
 533     int   i, j, n;
 534
 535     raw1 = gmx_ffopen(fnms[0], "w");
 536     raw2 = gmx_ffopen(fnms[1], "w");
 537     fprintf(raw1, "#Legend\n#TBlock\n#Xbin Ybin Z t\n");
 538     fprintf(raw2, "#Legend\n#TBlock\n#Xbin Ybin Z t\n");
 539     for (n = 0; n < tblocks; n++)
 540     {
 541         fprintf(raw1, "%5d\n", n);
 542         fprintf(raw2, "%5d\n", n);
 543         for (i = 0; i < xbins; i++)
 544         {
 545             for (j = 0; j < ybins; j++)
 546             {
 547                 fprintf(raw1, "%i  %i  %8.5f\n", i, j, (surf[0][n][j+ybins*i]));
 548                 fprintf(raw2, "%i  %i  %8.5f\n", i, j, (surf[1][n][j+ybins*i]));
 549             }
 550         }
 551     }
 552
 553     gmx_ffclose(raw1);
 554     gmx_ffclose(raw2);
 555 }
 556
 557
 558
 559 int gmx_hydorder(int argc, char *argv[])
 560 {
 561     static const char *desc[] = {
 562         "[THISMODULE] computes the tetrahedrality order parameters around a ",
 563         "given atom. Both angle an distance order parameters are calculated. See",
 564         "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
 565         "for more details.[PAR]"
 566         "[THISMODULE] calculates the order parameter in a 3d-mesh in the box, and",
 567         "with 2 phases in the box gives the user the option to define a 2D interface in time",
 568         "separating the faces by specifying parameters [TT]-sgang1[tt] and",
 569         "[TT]-sgang2[tt] (it is important to select these judiciously)."
 570     };
 571
 572     int                axis      = 0;
 573     static int         nsttblock = 1;
 574     static int         nlevels   = 100;
 575     static real        binwidth  = 1.0; /* binwidth in mesh           */
 576     static real        sg1       = 1;
 577     static real        sg2       = 1;   /* order parameters for bulk phases */
 578     static gmx_bool    bFourier  = FALSE;
 579     static gmx_bool    bRawOut   = FALSE;
 580     int                frames, xslices, yslices; /* Dimensions of interface arrays*/
 581     real            ***intfpos;                  /* Interface arrays (intfnr,t,xy) -potentially large */
 582     static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
 583
 584     t_pargs            pa[] = {
 585         { "-d",   FALSE, etENUM, {normal_axis},
 586           "Direction of the normal on the membrane" },
 587         { "-bw",  FALSE, etREAL, {&binwidth},
 588           "Binwidth of box mesh" },
 589         { "-sgang1", FALSE, etREAL, {&sg1},
 590           "tetrahedral angle parameter in Phase 1 (bulk)" },
 591         { "-sgang2", FALSE, etREAL, {&sg2},
 592           "tetrahedral angle parameter in Phase 2 (bulk)" },
 593         { "-tblock", FALSE, etINT, {&nsttblock},
 594           "Number of frames in one time-block average"},
 595         { "-nlevel", FALSE, etINT, {&nlevels},
 596           "Number of Height levels in 2D - XPixMaps"}
 597     };
 598
 599     t_filenm           fnm[] = {                      /* files for g_order    */
 600         { efTRX, "-f", NULL,  ffREAD },               /* trajectory file              */
 601         { efNDX, "-n", NULL,  ffREAD },               /* index file           */
 602         { efTPR, "-s", NULL,  ffREAD },               /* topology file                */
 603         { efXPM, "-o", "intf",  ffWRMULT},            /* XPM- surface maps      */
 604         { efOUT, "-or", "raw", ffOPTWRMULT },         /* xvgr output file           */
 605         { efOUT, "-Spect", "intfspect", ffOPTWRMULT}, /* Fourier spectrum interfaces */
 606     };
 607 #define NFILE asize(fnm)
 608
 609     /*Filenames*/
 610     const char       *ndxfnm, *tpsfnm, *trxfnm;
 611     char            **spectra, **intfn, **raw;
 612     int               nfspect, nfxpm, nfraw;
 613     gmx_output_env_t *oenv;
 614
 615     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 616                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 617     {
 618         return 0;
 619     }
 620     bFourier = opt2bSet("-Spect", NFILE, fnm);
 621     bRawOut  = opt2bSet("-or", NFILE, fnm);
 622
 623     if (binwidth < 0.0)
 624     {
 625         gmx_fatal(FARGS, "Can not have binwidth < 0");
 626     }
 627
 628     ndxfnm = ftp2fn(efNDX, NFILE, fnm);
 629     tpsfnm = ftp2fn(efTPR, NFILE, fnm);
 630     trxfnm = ftp2fn(efTRX, NFILE, fnm);
 631
 632     /* Calculate axis */
 633     GMX_RELEASE_ASSERT(normal_axis[0] != NULL, "Option setting inconsistency; normal_axis[0] is NULL");
 634     if (std::strcmp(normal_axis[0], "x") == 0)
 635     {
 636         axis = XX;
 637     }
 638     else if (std::strcmp(normal_axis[0], "y") == 0)
 639     {
 640         axis = YY;
 641     }
 642     else if (std::strcmp(normal_axis[0], "z") == 0)
 643     {
 644         axis = ZZ;
 645     }
 646     else
 647     {
 648         gmx_fatal(FARGS, "Invalid axis, use x, y or z");
 649     }
 650
 651     switch (axis)
 652     {
 653         case 0:
 654             fprintf(stderr, "Taking x axis as normal to the membrane\n");
 655             break;
 656         case 1:
 657             fprintf(stderr, "Taking y axis as normal to the membrane\n");
 658             break;
 659         case 2:
 660             fprintf(stderr, "Taking z axis as normal to the membrane\n");
 661             break;
 662     }
 663
 664     /* tetraheder order parameter */
 665     /* If either of the options is set we compute both */
 666     nfxpm = opt2fns(&intfn, "-o", NFILE, fnm);
 667     if (nfxpm != 2)
 668     {
 669         gmx_fatal(FARGS, "No or not correct number (2) of output-files: %d", nfxpm);
 670     }
 671     calc_tetra_order_interface(ndxfnm, tpsfnm, trxfnm, binwidth, nsttblock, &frames, &xslices, &yslices, sg1, sg2, &intfpos, oenv);
 672     writesurftoxpms(intfpos, frames, xslices, yslices, binwidth, intfn, nlevels);
 673
 674     if (bFourier)
 675     {
 676         nfspect = opt2fns(&spectra, "-Spect", NFILE, fnm);
 677         if (nfspect != 2)
 678         {
 679             gmx_fatal(FARGS, "No or not correct number (2) of output-files: %d", nfspect);
 680         }
 681         powerspectavg(intfpos, frames, xslices, yslices, spectra);
 682     }
 683
 684     if (bRawOut)
 685     {
 686         nfraw = opt2fns(&raw, "-or", NFILE, fnm);
 687         if (nfraw != 2)
 688         {
 689             gmx_fatal(FARGS, "No or not correct number (2) of output-files: %d", nfraw);
 690         }
 691         writeraw(intfpos, frames, xslices, yslices, raw);
 692     }
 693
 694     return 0;
 695 }