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