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