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