src/tools/gmx_hydorder.c

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