--- /dev/null
- #include <stdlib.h>
- #include <stdio.h>
- #include <string.h>
- #include "statutil.h"
- #include "typedefs.h"
- #include "smalloc.h"
- #include "vec.h"
- #include "copyrite.h"
- #include "tpxio.h"
+/*
+ * $Id: densorder.c,v 0.9
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.0
+ *
+ * Copyright (c) 1991-2001
+ * BIOSON Research Institute, Dept. of Biophysical Chemistry
+ * University of Groningen, The Netherlands
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
+ *
+ * And Hey:
+ * Gyas ROwers Mature At Cryogenic Speed
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <math.h>
- #include "pdbio.h"
+#include "sysstuff.h"
++#include "string.h"
+#include "string2.h"
++#include "typedefs.h"
++#include "smalloc.h"
+#include "macros.h"
+#include "gstat.h"
++#include "vec.h"
+#include "xvgr.h"
+#include "pbc.h"
++#include "copyrite.h"
+#include "futil.h"
++#include "statutil.h"
+#include "index.h"
- #include "gstat.h"
++#include "tpxio.h"
+#include "physics.h"
- void center_coords(t_atoms *atoms,matrix box,rvec x0[],int axis)
+#include "matio.h"
++#include "dens_filter.h"
+#include "binsearch.h"
+#include "powerspect.h"
++#include "gmx_ana.h"
++
++#ifdef GMX_DOUBLE
++#define FLOOR(x) ((int) floor(x))
++#else
++#define FLOOR(x) ((int) floorf(x))
++#endif
+
+enum {methSEL, methBISECT, methFUNCFIT, methNR};
+
-
- /*Code needs bounds checking and array-pointer debugging . 29/4*/
-
- static void density_in_time (const char *fn, const char *field,
- atom_id **index ,int nndx, int grpn, real binwidth,
- int nsttblock, real *****Densdevel, int *xslices, int *yslices, int *zslices,
- int *tblock, t_topology *top, int ePBC, int axis, bool bCenter,output_env_t oenv)
++static void center_coords(t_atoms *atoms,matrix box,rvec x0[],int axis)
+{
+ int i,m;
+ real tmass,mm;
+ rvec com,shift,box_center;
+
+ tmass = 0;
+ clear_rvec(com);
+ for(i=0; (i<atoms->nr); i++) {
+ mm = atoms->atom[i].m;
+ tmass += mm;
+ for(m=0; (m<DIM); m++)
+ com[m] += mm*x0[i][m];
+ }
+ for(m=0; (m<DIM); m++)
+ com[m] /= tmass;
+ calc_box_center(ecenterDEF,box,box_center);
+ rvec_sub(box_center,com,shift);
+ shift[axis] -= box_center[axis];
+
+ for(i=0; (i<atoms->nr); i++)
+ rvec_dec(x0[i],shift);
+}
+
+
- * binwidth widths of normal slices
++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, const output_env_t oenv)
+
+{
+/*
+ * *****Densdevel pointer to array of density values in slices and frame-blocks Densdevel[*nsttblock][*xslices][*yslices][*zslices]
+ * Densslice[x][y][z]
+ * nsttblock - nr of frames in each time-block
- int natoms, status,i,j,k,n, /* loop indices, checks etc*/
++ * bw widths of normal slices
+ *
+ * axis - axis direction (normal to slices)
+ * nndx - number ot atoms in **index
+ * grpn - group number in index
+ */
++ t_trxstatus *status;
++ gmx_rmpbc_t gpbc=NULL;
+ matrix box; /* Box - 3x3 -each step*/
+ rvec *x0; /* List of Coord without PBC*/
- slicex, slicey, slicez; /*slice # of x y z position */
- real ***Densslice; /* Density-slice in one frame*/
- real dscale; /*physical scaling factor*/
- real t,x,y,z; /* time and coordinates*/
++ int natoms,i,j,k,n, /* loop indices, checks etc*/
+ ax1=0,ax2=0, /* tangent directions */
+ framenr=0, /* frame number in trajectory*/
-
++ slicex, slicey, slicez; //slice # of x y z position
++ real ***Densslice; // Density-slice in one frame
++ real dscale; //physical scaling factor
++ real t,x,y,z; // time and coordinates
+ *tblock=0;/* blocknr in block average - initialise to 0*/
- gmx_fatal(FARGS, "Could not read init. coordinates!"); /* Open trajectory for read*/
++
+ /* Axis: X=0, Y=1,Z=2 */
+ switch(axis)
+ {
+ case 0:
+ ax1=YY; ax2=ZZ; /*Surface: YZ*/
+ break;
+ case 1:
+ ax1=ZZ; ax2=XX; /* Surface : XZ*/
+ break;
+ case 2:
+ ax1 = XX; ax2 = YY; /* Surface XY*/
+ break;
+ default:
+ gmx_fatal(FARGS,"Invalid axes. Terminating\n");
+ }
+
+ if( (natoms= read_first_x(oenv,&status,fn,&t,&x0,box)==0))
- *zslices=(int)(box[axis][axis]/binwidth+0.5);
- *yslices=(int)(box[ax2][ax2]/binwidth+0.5);
- *xslices=(int)(box[ax1][ax1]/binwidth+0.5);
++ gmx_fatal(FARGS, "Could not read coordinates from file"); // Open trajectory for read
+
+
- "\nDividing the box in %5d x %5d x %5d slices with binw %f along axis %d\n",*xslices,*yslices,*zslices,binwidth,axis );
-
-
++ *zslices=1+FLOOR(box[axis][axis]/bwz);
++ *yslices=1+FLOOR(box[ax2][ax2]/bw);
++ *xslices=1+FLOOR(box[ax1][ax1]/bw);
+ fprintf(stderr,
- /*Initialize Densdevel*/
-
++ "\nDividing the box in %5d x %5d x %5d slices with binw %f along axis %d\n",*xslices,*yslices,*zslices,bw,axis );
+
++
+ /****Start trajectory processing***/
+
-
- /*Reset Densslice every nsttblock steps*/
++ //Initialize Densdevel and PBC-remove
++ //gpbc=gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
++
+ *Densdevel=NULL;
+
+ do {
- srenew(*Densdevel,*tblock+1);
++ //gmx_rmpbc(gpbc,natoms,box,x0);
++ //Reset Densslice every nsttblock steps
+ if ( framenr % nsttblock==0 ){
+ snew(Densslice,*xslices);
+ for (i=0;i<*xslices;i++) {
+ snew(Densslice[i],*yslices);
+ for (j=0;j<*yslices;j++){
+ snew(Densslice[i][j],*zslices);
+ }
+ }
+
+ /*Allocate Memory to extra frame in Densdevel - rather stupid approach: *A single frame each time, although only every nsttblock steps.*/
-
-
++ srenew(*Densdevel,*tblock+1);
+
+ }
- rm_pbc(&(top->idef),ePBC,top->atoms.nr,box,x0,x0);
+
- for (j=0;j<nndx;j++) { /*Loop over all atoms in selected index*/
- x=x0[index[grpn][j]][ax1];
- y=x0[index[grpn][j]][ax2];
- z=x0[index[grpn][j]][axis];
+
+ dscale=(*xslices)*(*yslices)*(*zslices)*AMU/ (box[ax1][ax1]*box[ax2][ax2]*box[axis][axis]*nsttblock*(NANO*NANO*NANO));
+
+ if (bCenter)
+ center_coords(&top->atoms,box,x0,axis);
+
+
- slicex=((int) (x/box[ax1][ax1]* *xslices )) % *xslices;
- slicey=((int) (y/box[ax2][ax2]* *yslices)) % *yslices;
- slicez=((int) (z/box[axis][axis]* *zslices)) % *zslices;
- Densslice[slicex][slicey][slicez]+=(top->atoms.atom[index[grpn][j]].m*dscale);
++ for (j=0;j<gnx[0];j++) { /*Loop over all atoms in selected index*/
++ x=x0[index[0][j]][ax1];
++ y=x0[index[0][j]][ax2];
++ z=x0[index[0][j]][axis];
+ while (x<0)
+ x+=box[ax1][ax1];
+ while(x>box[ax1][ax1])
+ x-=box[ax1][ax1];
+
+ while (y<0)
+ y+=box[ax2][ax2];
+ while(y>box[ax2][ax2])
+ y-=box[ax2][ax2];
+
+ while (z<0)
+ z+=box[axis][axis];
+ while(z>box[axis][axis])
+ z-=box[axis][axis];
+
- /*Implicit incrementation of Densdevel via renewal of Densslice*/
- /*only every nsttblock steps*/
++ slicex=(int) (x/bw);
++ slicey=(int) (y/bw);
++ slicez=(int) (z/bwz);
++ Densslice[slicex][slicey][slicez]+=(top->atoms.atom[index[0][j]].m*dscale);
+
+
+ }
+
+ framenr++;
+
+ if(framenr % nsttblock == 0){
- } while(read_next_x(oenv,status,&t,natoms,x0,box));
++ //Implicit incrementation of Densdevel via renewal of Densslice
++ //only every nsttblock steps
+ (*Densdevel)[*tblock]=Densslice;
+ (*tblock)++;
+ }
+
-
-
++ } while(read_next_x(oenv,status,&t,natoms,x0,box));
+
- /*Need to sum and divide Densdevel array with constant factor blocknr for each element to average in time.*/
- if(NULL != field){
- /*Debug-filename and filehandle*/
- FILE *fldH;
- char fff[STRLEN];
- int nnn=1;
- real ddd;
-
- sprintf(fff,"%s%%6.2f%%6.2f\n",pdbformat);
- fldH=ffopen(field,"w");
- for(i=0;i<*xslices;i++){
- for (j=0;j<*yslices;j++){
- for (k=0;k<*zslices;k++){
- ddd = 0;
- for(n=0;n<*tblock;n++){
- ddd += (*Densdevel)[n][i][j][k];
- }
- ddd = 0.01*ddd/(*tblock);
- fprintf(fldH,fff,"ATOM",nnn++,"HE","HE",' ',1,'A',4.0*i,4.0*j,4.0*k,1.0,ddd);
- }
- }
- }
-
- ffclose(fldH);
- }
- /*Free memory we no longer need.*/
- sfree(x0);
-
+
- static void printdensavg( real ****Densmap, real binwidth, int xslices, int yslices, int zslices, int tdim)
++ //Free memory we no longer need and exit.
++ //gmx_rmpbc_done(gpbc);
++ close_trj(status);
+}
+
+
+
- int n,i,j,k;
- real totdens=0;
++static void printdensavg(char *fldfn, real ****Densmap, int xslices, int yslices, int zslices, int tdim)
+{
- /*Allocate memory to array*/
-
- for(n=0;n<tdim;n++){
- for(i=0;i<xslices;i++){
- for (j=0;j<yslices;j++){
- for (k=0;k<zslices;k++){
- totdens+=((Densmap[n][i][j][k])/(xslices*yslices*zslices*tdim));
- }
- }
++/*Debug-filename and filehandle*/
++ FILE *fldH;
++ int n,i,j,k;
++ real totdens=0;
++ fldH=ffopen(fldfn,"w");
++ fprintf(fldH,"%i %i %i %i\n",tdim, xslices,yslices,zslices);
++ for(n=0;n<tdim;n++){
++ for(i=0;i<xslices;i++){
++ for (j=0;j<yslices;j++){
++ for (k=0;k<zslices;k++){
++ fprintf(fldH,"%i %i %i %f\n",i,j,k,Densmap[n][i][j][k]);
++ totdens+=(Densmap[n][i][j][k]);
++ }
++ }
++ }
++ }
++totdens/=(xslices*yslices*zslices*tdim);
++fprintf(stderr,"Total density [kg/m^3] %8f",totdens);
++ffclose(fldH);
++}
+
-
- fprintf(stderr,"Total density [kg/m^3] %8f",totdens);
++
++static void outputfield(char *fldfn, real ****Densmap, int xslices, int yslices, int zslices, int tdim)
++{
++/*Debug-filename and filehandle*/
++ FILE *fldH;
++ int n,i,j,k;
++ int dim[4]={tdim,xslices,yslices,zslices};
++ real totdens=0;
++ fldH=ffopen(fldfn,"w");
++ fwrite(dim,sizeof(int),4,fldH);
++ for(n=0;n<tdim;n++){
++ for(i=0;i<xslices;i++){
++ for (j=0;j<yslices;j++){
++ for (k=0;k<zslices;k++){
++ fwrite(&(Densmap[n][i][j][k]),sizeof(real),1,fldH);
++ totdens+=(Densmap[n][i][j][k]);
+ }
+ }
- int tblocks, real binwidth,int method, real dens1, real dens2, t_interf ****intf1, t_interf ****intf2,output_env_t oenv)
++ }
++ }
++totdens/=(xslices*yslices*zslices*tdim);
++fprintf(stderr,"Total density [kg/m^3] %8f",totdens);
++ffclose(fldH);
++}
+
++
++static void filterdensmap(real ****Densmap, int xslices, int yslices, int zslices,int tblocks,int ftsize)
++{
++ real *kernel;
++ real *output;
++ real std,var;
++ int i,j,k,n, order;
++ order=ftsize/2;
++ std=((real)order/2.0);
++ var=std*std;
++ snew(kernel,ftsize);
++ gausskernel(kernel,ftsize,var);
++ for(n=0;n<tblocks;n++){
++ for(i=0;i<xslices;i++){
++ for (j=0;j<yslices;j++){
++ snew(output,zslices);
++ convolve1D( Densmap[n][i][j],output,zslices,kernel,ftsize );
++ Densmap[n][i][j]=output;
++ }
++ }
++ }
+}
+
++
++
+
+static void interfaces_txy (real ****Densmap, int xslices, int yslices, int zslices,
- /*Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)*/
++ int tblocks, real binwidth,int method, real dens1, real dens2, t_interf ****intf1, t_interf ****intf2,const output_env_t oenv)
+{
+
- real *zDensavg; /* zDensavg[z]*/
++//Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)
+ FILE *xvg;
- t_interf ***int1=NULL,***int2=NULL; /*Interface matrices [t][x,y] - last index in row-major order*/
++ real *zDensavg; // zDensavg[z]
+ int i,j,k,n;
+ int xysize;
+ int ndx1, ndx2, deltandx, *zperm;
+ real densmid, densl, densr, alpha, pos, spread;
+ real splitpoint, startpoint, endpoint;
+ real *sigma1, *sigma2;
+ const real onehalf= 1.00/2.00;
- rangeArray(zperm,zslices); /*reset permutation array to identity*/
++ t_interf ***int1=NULL,***int2=NULL; //Interface matrices [t][x,y] - last index in row-major order
+ /*Create int1(t,xy) and int2(t,xy) arrays with correct number of interf_t elements*/
+ xysize=xslices*yslices;
+ snew(int1,tblocks);
+ snew(int2,tblocks);
+ for (i=0;i<tblocks;i++){
+ snew(int1[i],xysize);
+ snew(int2[i],xysize);
+ for(j=0;j<xysize;j++){
+ snew(int1[i][j],1);
+ snew(int2[i][j],1);
+ init_interf(int1[i][j]);
+ init_interf(int2[i][j]);
+ }
+ }
+
+if(method==methBISECT){
+ densmid= onehalf*(dens1+dens2);
+ snew(zperm,zslices);
+ for(n=0;n<tblocks;n++){
+ for(i=0;i<xslices;i++){
+ for (j=0;j<yslices;j++){
- /* rho_11= Densmap[n][i][j][zperm[ndx1]]
- * rho_12 =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
- * deltandx=zperm[ndx1+1]-zperm[ndx+1]
- * * alpha=(densmid-rho_11)/(rho_12-rho_11) */
++ rangeArray(zperm,zslices); //reset permutation array to identity
+ /*Binsearch returns slice-nr where the order param is <= setpoint sgmid*/
+ ndx1=start_binsearch(Densmap[n][i][j],zperm,0,zslices/2-1,densmid,1);
+ ndx2=start_binsearch(Densmap[n][i][j],zperm,zslices/2,zslices-1,densmid,-1);
- printf("Alpha, Deltandx %f %i\n", alpha,deltandx);
++
++ /* Linear interpolation (for use later if time allows)
++ * rho_1s= Densmap[n][i][j][zperm[ndx1]]
++ * rho_1e =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
++ * rho_2s =Densmap[n][i][j][zperm[ndx2+1]]
++ * rho_2e =Densmap[n][i][j][zperm[ndx2]]
++ * For 1st interface we have:
+ densl= Densmap[n][i][j][zperm[ndx1]];
+ densr= Densmap[n][i][j][zperm[ndx1+1]];
+ alpha=(densmid-densl)/(densr-densl);
+ deltandx=zperm[ndx1+1]-zperm[ndx1];
- int1[n][j+(i*yslices)]->Z=(pos+onehalf)*binwidth;
- int1[n][j+(i*yslices)]->t=spread;
- /*Leave the 2nd interface as test at the moment -needs to be made similar to 1st later.
- *We can use the same formulation, since alpha should become negative ie:
- alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/
- (Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
- deltandx=zperm[ndx2+1]-zperm[ndx2];
- pos=zperm[ndx2]+alpha*deltandx; */
++
++ if(debug){
++ printf("Alpha, Deltandx %f %i\n", alpha,deltandx);
++ }
+ if(abs(alpha)>1.0 || abs(deltandx)>3){
+ pos=zperm[ndx1];
+ spread=-1;
+ }
+ else {
+ pos=zperm[ndx1]+alpha*deltandx;
+ spread=binwidth*deltandx;
+ }
- /*Assume a box divided in 2 along midpoint of z for starters*/
++ * For the 2nd interface can use the same formulation, since alpha should become negative ie:
++ * alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/(Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
++ * deltandx=zperm[ndx2+1]-zperm[ndx2];
++ * pos=zperm[ndx2]+alpha*deltandx; */
++
++ //After filtering we use the direct approach
++ int1[n][j+(i*yslices)]->Z=(zperm[ndx1]+onehalf)*binwidth;
++ int1[n][j+(i*yslices)]->t=binwidth;
+ int2[n][j+(i*yslices)]->Z=(zperm[ndx2]+onehalf)*binwidth;
+ int2[n][j+(i*yslices)]->t=binwidth;
+ }
+ }
+ }
+}
+
+if(method==methFUNCFIT){
- /*Initial fit proposals*/
++ //Assume a box divided in 2 along midpoint of z for starters
+ startpoint=0.0;
+ endpoint = binwidth*zslices;
+ splitpoint = (startpoint+endpoint)/2.0;
- /*Calculate average density along z - avoid smoothing by using coarse-grained-mesh*/
++ //Initial fit proposals
+ real beginfit1[4] = {dens1, dens2, (splitpoint/2), 0.5};
+ real beginfit2[4] = {dens2, dens1, (3*splitpoint/2), 0.5};
+ real *fit1=NULL,*fit2=NULL;
+
+ snew(zDensavg,zslices);
+ snew(sigma1,zslices);
+ snew(sigma2,zslices);
+
+ for(k=0;k<zslices;k++) sigma1[k]=sigma2[k]=1;
- xvg=xvgropen("DensprofileonZ.xvg", "Averaged Densityprofile on Z","z[nm]","Density[kg/m^3]",oenv);
++ //Calculate average density along z - avoid smoothing by using coarse-grained-mesh
+ for(k=0;k<zslices;k++){
+ for(n=0;n<tblocks;n++){
+ for (i=0;i<xslices;i++){
+ for(j=0;j<yslices;j++){
+ zDensavg[k]+=(Densmap[n][i][j][k]/(xslices*yslices*tblocks));
+ }
+ }
+ }
+ }
+
+ if(debug){
- /*Fit 1st half of box*/
++ xvg=xvgropen("DensprofileonZ.xvg", "Averaged Densityprofile on Z","z[nm]","Density[kg/m^3]",oenv);
+ for(k=0;k<zslices;k++) fprintf(xvg, "%4f.3 %8f.4\n", k*binwidth,zDensavg[k]);
+ fclose(xvg);
+ }
+
+ /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
+
- /*Fit 2nd half of box*/
++ //Fit 1st half of box
+ do_lmfit(zslices, zDensavg,sigma1,binwidth,NULL,startpoint,splitpoint,oenv,FALSE,effnERF,beginfit1,3);
- /*Initialise the const arrays for storing the average fit parameters*/
++ //Fit 2nd half of box
+ do_lmfit(zslices ,zDensavg,sigma2,binwidth,NULL,splitpoint,endpoint,oenv,FALSE,effnERF,beginfit2,3);
+
- /*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*/
++ //Initialise the const arrays for storing the average fit parameters
+ const real * const avgfit1=beginfit1;
+ const real * const avgfit2=beginfit2;
+
+
+
- /*Reinitialise fit for each mesh-point*/
++ //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
+ for(n=0;n<tblocks;n++){
+ for(i=0;i<xslices;i++){
+ for (j=0;j<yslices;j++){
- /*Now fit and store in structures in row-major order int[n][i][j]*/
- do_lmfit(zslices,Densmap[n][i][j],sigma1,binwidth,NULL,startpoint,splitpoint,
- oenv,FALSE,effnERF,fit1,1);
++ //Reinitialise fit for each mesh-point
+ srenew(fit1,4);
+ srenew(fit2,4);
+ for (k=0;k<4;k++){
+ fit1[k]=avgfit1[k];
+ fit2[k]=avgfit2[k];
+ }
- do_lmfit(zslices,Densmap[n][i][j],sigma2,binwidth,NULL,splitpoint,endpoint,
- oenv,FALSE,effnERF,fit2,2);
++ //Now fit and store in structures in row-major order int[n][i][j]
++ do_lmfit(zslices,Densmap[n][i][j],sigma1,binwidth,NULL,startpoint,splitpoint,oenv,FALSE,effnERF,fit1,1);
+ int1[n][j+(yslices*i)]->Z=fit1[2];
+ int1[n][j+(yslices*i)]->t=fit1[3];
- void writesurftoxpms(t_interf ***surf1,t_interf ***surf2, int tblocks,int xbins, int ybins, int zbins,
- real bw,char **outfiles,int maplevels )
++ do_lmfit(zslices,Densmap[n][i][j],sigma2,binwidth,NULL,splitpoint,endpoint,oenv,FALSE,effnERF,fit2,2);
+ int2[n][j+(yslices*i)]->Z=fit2[2];
+ int2[n][j+(yslices*i)]->t=fit2[3];
+ }
+ }
+ }
+}
+
+
+*intf1=int1;
+*intf2=int2;
+
+}
+
- char numbuf[16];
++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 )
+{
- /*Allocate memory to tick's and matrices*/
++ char numbuf[12];
+ int n, i, j;
+ real **profile1, **profile2;
+ real max1, max2, min1, min2, *xticks, *yticks;
+ t_rgb lo={0,0,0};
+ t_rgb hi={1,1,1};
+ FILE *xpmfile1, *xpmfile2;
+
+/*Prepare xpm structures for output*/
+
- min1=min2=zbins*bw;
++//Allocate memory to tick's and matrices
+snew (xticks,xbins+1);
+snew (yticks,ybins+1);
+
+profile1=mk_matrix(xbins,ybins,FALSE);
+profile2=mk_matrix(xbins,ybins,FALSE);
+
+for (i=0;i<xbins+1;i++) xticks[i]+=bw;
+for (j=0;j<ybins+1;j++) yticks[j]+=bw;
+
+xpmfile1 = ffopen(outfiles[0],"w");
+xpmfile2 = ffopen(outfiles[1],"w");
+
+max1=max2=0.0;
- /*Filling matrices for inclusion in xpm-files*/
++min1=min2=zbins*bwz;
+
+for(n=0;n<tblocks;n++){
+sprintf(numbuf,"tblock: %4i",n);
- /*Finding max and min values*/
++//Filling matrices for inclusion in xpm-files
+ for(i=0;i<xbins;i++){
+ for(j=0;j<ybins;j++){
+ profile1[i][j]=(surf1[n][j+ybins*i])->Z;
+ profile2[i][j]=(surf2[n][j+ybins*i])->Z;
- void writeraw(t_interf ***int1, t_interf ***int2, int tblocks,int xbins, int ybins,char **fnms){
++ //Finding max and min values
+ if(profile1[i][j]>max1) max1=profile1[i][j];
+ if(profile1[i][j]<min1) min1=profile1[i][j];
+ if(profile2[i][j]>max2) max2=profile2[i][j];
+ if(profile2[i][j]<min2) min2=profile2[i][j];
+ }
+ }
+
+ write_xpm(xpmfile1,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile1,min1,max1,lo,hi,&maplevels);
+ write_xpm(xpmfile2,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile2,min2,max2,lo,hi,&maplevels);
+}
+
+ffclose(xpmfile1);
+ffclose(xpmfile2);
+
+
+sfree(profile1);
+sfree(profile2);
+sfree(xticks);
+sfree(yticks);
+}
+
- fprintf(raw1,"#Legend\n#TBlock\n#Xbin Ybin Z t\n");
- fprintf(raw2,"#Legend\n#TBlock\n#Xbin Ybin Z t\n");
++static void writeraw(t_interf ***int1, t_interf ***int2, int tblocks,int xbins, int ybins,char **fnms){
+ FILE *raw1, *raw2;
+ int i,j,n;
+
+ raw1=ffopen(fnms[0],"w");
+ raw2=ffopen(fnms[1],"w");
+ fprintf(raw1,"#Produced by: %s #\n", command_line());
+ fprintf(raw2,"#Produced by: %s #\n", command_line());
- fprintf(raw1,"%5d\n",n);
- fprintf(raw2,"%5d\n",n);
++ fprintf(raw1,"#Legend: nt nx ny\n#Xbin Ybin Z t\n");
++ fprintf(raw2,"#Legend: nt nx ny\n#Xbin Ybin Z t\n");
++ fprintf(raw1,"%i %i %i\n", tblocks,xbins,ybins);
++ fprintf(raw2,"%i %i %i\n", tblocks,xbins,ybins);
+ for (n=0;n<tblocks;n++){
- "A small program to reduce a two-phase density distribution",
+ for(i=0;i<xbins;i++){
+ for(j=0;j<ybins;j++){
+ fprintf(raw1,"%i %i %8.5f %6.4f\n",i,j,(int1[n][j+ybins*i])->Z,(int1[n][j+ybins*i])->t);
+ fprintf(raw2,"%i %i %8.5f %6.4f\n",i,j,(int2[n][j+ybins*i])->Z,(int2[n][j+ybins*i])->t);
+ }
+ }
+ }
+
+ ffclose(raw1);
+ ffclose(raw2);
+}
+
+
+
+int gmx_densorder(int argc,char *argv[])
+{
+ static const char *desc[] = {
-
++ "A small program to reduce a two-phase densitydistribution",
+ "along an axis, computed over a MD trajectory",
+ "to 2D surfaces fluctuating in time, by a fit to",
+ "a functional profile for interfacial densities",
+ "A time-averaged spatial representation of the" ,
+ "interfaces can be output with the option -tavg"
+ };
+
+ /* Extra arguments - but note how you always get the begin/end
+ * options when running the program, without mentioning them here!
+ */
+
- atom_id **index; /* Index list for single group*/
++ output_env_t oenv;
+ t_topology *top;
- static real binwidth=0.2;
+ char title[STRLEN],**grpname;
+ int ePBC, *ngx;
- static bool bCenter = FALSE;
- static bool bFourier=FALSE;
- static bool bAvgt = FALSE;
- static bool bRawOut=FALSE;
++ static real binw=0.2;
++ static real binwz=0.05;
+ static real dens1=0.00;
+ static real dens2=1000.00;
++ static int ftorder=1;
+ static int nsttblock=100;
+ static int axis= 2;
+ static char *axtitle="Z";
+ static int ngrps=1;
++ atom_id **index; // Index list for single group
+ int xslices, yslices, zslices, tblock;
- /*Densitymap - Densmap[t][x][y][z]*/
++ static gmx_bool bCenter = FALSE;
++ static gmx_bool bFourier=FALSE;
++ static gmx_bool bAvgt = FALSE;
++ static gmx_bool bRawOut=FALSE;
++ static gmx_bool bOut=FALSE;
+ static int nlevels=100;
- /* Surfaces surf[t][surf_x,surf_y]*/
++//Densitymap - Densmap[t][x][y][z]
+ real ****Densmap=NULL;
- static const char *meth[methNR+1]={NULL,"bisect","functional",NULL};
++// Surfaces surf[t][surf_x,surf_y]
+ t_interf ***surf1, ***surf2;
+
- { "-bw",FALSE,etREAL,{&binwidth},
- "Binwidth of density distribution along axis"},
++ static const char *meth[]={NULL,"bisect","functional",NULL};
+ int eMeth;
+
+ char **outfiles, **rawfiles, **spectra; /* Filenames for xpm-surface maps, rawdata and powerspectra */
+ int nfxpm,nfraw, nfspect; /* # files for interface maps and spectra = # interfaces */
+
+t_pargs pa[] = {
+ { "-tavg", FALSE, etBOOL, {&bAvgt},
+ "Plot time averaged interface profile"},
- { efXPM, "-o" ,"interface",ffWRMULT},/* This is for writing out the interface meshes - one xpm-file per tblock*/
- { efOUT, "-Spect","intfspect",ffOPTWRMULT}, /* This is for the trajectory averaged Fourier-spectra*/
- { efPDB, "-field","densfield",ffOPTWR}
++ { "-bw",FALSE,etREAL,{&binw},
++ "Binwidth of density distribution tangential to interface"},
++ { "-bwn", FALSE, etREAL, {&binwz},
++ "Binwidth of density distribution normal to interface"},
++ { "-order", FALSE, etINT,{&ftorder},
++ "Order of Gaussian filter"},
+ {"-axis", FALSE, etSTR,{&axtitle},
+ "Axis Direction - X, Y or Z"},
+ {"-method",FALSE ,etENUM,{meth},
+ "Interface location method"},
+ {"-d1", FALSE, etREAL,{&dens1},
+ "Bulk density phase 1 (at small z)"},
+ {"-d2", FALSE, etREAL,{&dens2},
+ "Bulk density phase 2 (at large z)"},
+ { "-tblock",FALSE,etINT,{&nsttblock},
+ "Number of frames in one time-block average"},
+ { "-nlevel", FALSE,etINT, {&nlevels},
+ "Number of Height levels in 2D - XPixMaps"}
+ };
+
+
+ t_filenm fnm[] = {
+ { efTPX, "-s", NULL, ffREAD }, /* this is for the topology */
+ { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
+ { efNDX, "-n", NULL, ffREAD}, /* this is to select groups */
+ { efOUT, "-or", NULL,ffOPTWRMULT}, /* This is for writing out the entire information in the t_interf arrays */
- output_env_t oenv;
-
++ { efXPM, "-o" ,"interface",ffOPTWRMULT},/* This is for writing out the interface meshes - one xpm-file per tblock*/
++ { efOUT, "-Spect","intfspect",ffOPTWRMULT}, /* This is for the trajectory averaged Fourier-spectra*/
+ };
+
+#define NFILE asize(fnm)
-
++
+ CopyRight(stderr,argv[0]);
+
+ /* This is the routine responsible for adding default options,
+ * calling the X/motif interface, etc. */
+ parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,
+ NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
+
+
+eMeth=nenum(meth);
+bFourier=opt2bSet("-Spect",NFILE,fnm);
+bRawOut=opt2bSet("-or",NFILE,fnm);
- density_in_time(ftp2fn(efTRX,NFILE,fnm),ftp2fn_null(efPDB,NFILE,fnm),
- index,ngx[0],0,binwidth,nsttblock,&Densmap,
- &xslices,&yslices,&zslices,&tblock, top,ePBC,axis,bCenter,oenv);
++bOut=opt2bSet("-o",NFILE,fnm);
+top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
+snew(grpname,ngrps);
+snew(index,ngrps);
+snew(ngx,ngrps);
+
+/* Calculate axis */
+ axis = toupper(axtitle[0]) - 'X';
+
+get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),ngrps,ngx,index,grpname);
+
- interfaces_txy(Densmap,xslices,yslices,zslices,tblock,binwidth,eMeth,dens1,dens2,&surf1,&surf2,oenv);
++density_in_time(ftp2fn(efTRX,NFILE,fnm),index,ngx,ngrps,binw,binwz,nsttblock,&Densmap,&xslices,&yslices,&zslices,&tblock,top,ePBC,axis,bCenter,oenv);
++if (debug){
++outputfield("Density_4D.dat",Densmap,xslices,yslices,zslices,tblock);
++}
+
- /*Output surface-xpms*/
- nfxpm=opt2fns(&outfiles,"-o",NFILE,fnm);
- if(nfxpm!=2){
- gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
++filterdensmap(Densmap,xslices,yslices,zslices,tblock,2*ftorder+1);
+
- writesurftoxpms(surf1, surf2, tblock,xslices,yslices,zslices,binwidth , outfiles,zslices);
++if (debug){
++outputfield("Density_4D_filtered.dat",Densmap,xslices,yslices,zslices,tblock);
+}
- if (debug){
- printdensavg(Densmap,binwidth,xslices,yslices,zslices,tblock);
+
++interfaces_txy(Densmap,xslices,yslices,zslices,tblock,binwz,eMeth,dens1,dens2,&surf1,&surf2,oenv);
+
- /*Output raw-data*/
++if(bOut){
++
++//Output surface-xpms
++ nfxpm=opt2fns(&outfiles,"-o",NFILE,fnm);
++ if(nfxpm!=2){
++ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
++ }
++writesurftoxpms(surf1, surf2, tblock,xslices,yslices,zslices,binw,binwz,outfiles,zslices);
+}
+
++
+
+
+
-
++//Output raw-data
+if (bRawOut){
+ nfraw=opt2fns(&rawfiles,"-or",NFILE,fnm);
+ if(nfraw!=2){
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
+ }
+ writeraw(surf1,surf2,tblock,xslices,yslices,rawfiles);
+}
+
+
+
+if(bFourier){
+ nfspect=opt2fns(&spectra,"-Spect",NFILE,fnm);
+ if(nfspect!=2){
+ gmx_fatal(FARGS,"No or not correct number (2) of output-file-series: %d"
+ ,nfspect);
+ }
+ powerspectavg_intf(surf1, surf2, tblock,xslices,yslices,spectra);
+}
+
+sfree(Densmap);
++if(bOut || bFourier || bRawOut){
+sfree(surf1);
+sfree(surf2);
++}
+
+ thanx(stderr);
+ return 0;
+}
--- /dev/null
- set_pbc(&pbc,ePBC,box);
- rm_pbc(&(top.idef),ePBC,natoms,box,x,x);
+/*
+ * $Id: gmx_order.c,v 1.11 2009/03/11 12:13:15 spoel Exp $
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Green Red Orange Magenta Azure Cyan Skyblue
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <ctype.h>
+
+#include "sysstuff.h"
+#include "string.h"
+#include "typedefs.h"
++#include "statutil.h"
+#include "smalloc.h"
+#include "macros.h"
+#include "gstat.h"
+#include "vec.h"
+#include "xvgr.h"
+#include "pbc.h"
+#include "copyrite.h"
+#include "futil.h"
+#include "statutil.h"
+#include "index.h"
+#include "tpxio.h"
+#include "matio.h"
+#include "binsearch.h"
+#include "powerspect.h"
+// #include "cmat.h"
+/****************************************************************************/
+/* This program calculates the order parameter per atom for an interface or */
+/* bilayer, averaged over time. */
+/* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
+/* It is assumed that the order parameter with respect to a box-axis */
+/* is calculated. In that case i = j = axis, and delta(ij) = 1. */
+/* */
+/* Peter Tieleman, April 1995 */
+/* P.J. van Maaren, November 2005 Added tetrahedral stuff */
+/****************************************************************************/
+
+
+
+
+
+
+/* Print name of first atom in all groups in index file */
+static void print_types(atom_id index[], atom_id a[], int ngrps,
+ char *groups[], t_topology *top)
+{
+ int i;
+
+ fprintf(stderr,"Using following groups: \n");
+ for(i = 0; i < ngrps; i++)
+ fprintf(stderr,"Groupname: %s First atomname: %s First atomnr %u\n",
+ groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
+ fprintf(stderr,"\n");
+}
+
+static void check_length(real length, int a, int b)
+{
+ if (length > 0.3)
+ fprintf(stderr,"WARNING: distance between atoms %d and "
+ "%d > 0.3 nm (%f). Index file might be corrupt.\n",
+ a, b, length);
+}
+
+static void find_tetra_order_grid(t_topology top, int ePBC,
+ int natoms, matrix box,
+ rvec x[],int maxidx,atom_id index[],
+ real time, real *sgmean, real *skmean,
+ int nslicex, int nslicey, int nslicez,
+ real ***sggrid,real ***skgrid)
+{
+ int ix,jx,i,j,k,l,n,*nn[4];
+ rvec dx,rj,rk,urk,urj;
+ real cost,cost2,*sgmol,*skmol,rmean,rmean2,r2,box2,*r_nn[4];
+ t_pbc pbc;
+ //t_mat *dmat;
+ //t_dist *d;
+ int slindex_x,slindex_y,slindex_z;
+ int ***sl_count;
+ real onethird=1.0/3.0;
++ gmx_rmpbc_t gpbc;
+
+ /* dmat = init_mat(maxidx, FALSE); */
+
+ box2 = box[XX][XX] * box[XX][XX];
+
+ // Initialize expanded sl_count array
+ snew(sl_count,nslicex);
+ for(i=0;i<nslicex;i++){
+ snew(sl_count[i],nslicey);
+ for(j=0;j<nslicey;j++){
+ snew(sl_count[i][j],nslicez);
+ }
+ }
+
+
+ for (i=0; (i<4); i++) {
+ snew(r_nn[i],natoms);
+ snew(nn[i],natoms);
+
+ for (j=0; (j<natoms); j++) {
+ r_nn[i][j] = box2;
+ }
+ }
+
+ snew(sgmol,maxidx);
+ snew(skmol,maxidx);
+
+ /* Must init pbc every step because of pressure coupling */
- static void calc_tetra_order_interface(const char *fnNDX,const char *fnTPS,const char *fnTRX,
- real binw, int tblock,
- int *nframes, int *nslicex, int *nslicey,
++ gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
++ gmx_rmpbc(gpbc,natoms,box,x);
+
+ *sgmean = 0.0;
+ *skmean = 0.0;
+ l=0;
+ for (i=0; (i<maxidx); i++) { /* loop over index file */
+ ix = index[i];
+ for (j=0; (j<maxidx); j++) {
+ if (i == j) continue;
+
+ jx = index[j];
+
+ pbc_dx(&pbc,x[ix],x[jx],dx);
+ r2=iprod(dx,dx);
+
+ /* set_mat_entry(dmat,i,j,r2); */
+
+ /* determine the nearest neighbours */
+ if (r2 < r_nn[0][i]) {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
+ r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
+ r_nn[0][i] = r2; nn[0][i] = j;
+ } else if (r2 < r_nn[1][i]) {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
+ r_nn[1][i] = r2; nn[1][i] = j;
+ } else if (r2 < r_nn[2][i]) {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r2; nn[2][i] = j;
+ } else if (r2 < r_nn[3][i]) {
+ r_nn[3][i] = r2; nn[3][i] = j;
+ }
+ }
+
+
+ /* calculate mean distance between nearest neighbours */
+ rmean = 0;
+ for (j=0; (j<4); j++) {
+ r_nn[j][i] = sqrt(r_nn[j][i]);
+ rmean += r_nn[j][i];
+ }
+ rmean /= 4;
+
+ n = 0;
+ sgmol[i] = 0.0;
+ skmol[i] = 0.0;
+
+ /* Chau1998a eqn 3 */
+ /* angular part tetrahedrality order parameter per atom */
+ for (j=0; (j<3); j++) {
+ for (k=j+1; (k<4); k++) {
+ pbc_dx(&pbc,x[ix],x[index[nn[k][i]]],rk);
+ pbc_dx(&pbc,x[ix],x[index[nn[j][i]]],rj);
+
+ unitv(rk,urk);
+ unitv(rj,urj);
+
+ cost = iprod(urk,urj) + onethird;
+ cost2 = cost * cost;
+
+ sgmol[i] += cost2;
+ l++;
+ n++;
+ }
+ }
+ /* normalize sgmol between 0.0 and 1.0 */
+ sgmol[i] = 3*sgmol[i]/32;
+ *sgmean += sgmol[i];
+
+ /* distance part tetrahedrality order parameter per atom */
+ rmean2 = 4 * 3 * rmean * rmean;
+ for (j=0; (j<4); j++) {
+ skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
+ /* printf("%d %f (%f %f %f %f) \n",
+ i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
+ */
+ }
+
+ *skmean += skmol[i];
+
+ /* Compute sliced stuff in x y z*/
+ slindex_x = gmx_nint((1+x[i][XX]/box[XX][XX])*nslicex) % nslicex;
+ slindex_y = gmx_nint((1+x[i][YY]/box[YY][YY])*nslicey) % nslicey;
+ slindex_z = gmx_nint((1+x[i][ZZ]/box[ZZ][ZZ])*nslicez) % nslicez;
+ sggrid[slindex_x][slindex_y][slindex_z] += sgmol[i];
+ skgrid[slindex_x][slindex_y][slindex_z] += skmol[i];
+ (sl_count[slindex_x][slindex_y][slindex_z])++;
+ } /* loop over entries in index file */
+
+ *sgmean /= maxidx;
+ *skmean /= maxidx;
+
+ for(i=0; (i<nslicex); i++) {
+ for(j=0; j<nslicey; j++){
+ for(k=0;k<nslicez;k++){
+ if (sl_count[i][j][k] > 0) {
+ sggrid[i][j][k] /= sl_count[i][j][k];
+ skgrid[i][j][k] /= sl_count[i][j][k];
+ }
+ }
+ }
+ }
+
+ sfree(sl_count);
+ sfree(sgmol);
+ sfree(skmol);
+ for (i=0; (i<4); i++) {
+ sfree(r_nn[i]);
+ sfree(nn[i]);
+ }
+}
+
+/*Determines interface from tetrahedral order parameter in box with specified binwidth. */
+/*Outputs interface positions(bins), the number of timeframes, and the number of surface-mesh points in xy*/
+
- int status;
++static void calc_tetra_order_interface(const char *fnNDX,const char *fnTPS,const char *fnTRX, real binw, int tblock,
++ int *nframes, int *nslicex, int *nslicey,
+ real sgang1,real sgang2,real ****intfpos,
+ output_env_t oenv)
+{
+ FILE *fpsg=NULL,*fpsk=NULL;
+ char *sgslfn="sg_ang_mesh"; //Hardcoded filenames for debugging
+ char *skslfn="sk_dist_mesh";
+ t_topology top;
+ int ePBC;
+ char title[STRLEN],subtitle[STRLEN];
- void writesurftoxpms(real ***surf, int tblocks,int xbins, int ybins, real bw,char **outfiles,int maplevels )
++ t_trxstatus *status;
+ int natoms;
+ real t;
+ rvec *xtop,*x;
+ matrix box;
+ real sg,sk, sgintf, pos;
+ atom_id **index;
+ char **grpname;
+ int i,j,k,n,*isize,ng, nslicez, framenr;
+ real ***sg_grid,***sk_grid, ***sg_fravg, ***sk_fravg, ****sk_4d, ****sg_4d;
+ int *perm;
+ int ndx1, ndx2;
+ int bins;
+ const real onehalf=1.0/2.0;
+ /* real ***intfpos[2]; pointers to arrays of two interface positions zcoord(framenr,xbin,ybin): intfpos[interface_index][t][nslicey*x+y]
+ * i.e 1D Row-major order in (t,x,y) */
+
+
+ read_tps_conf(fnTPS,title,&top,&ePBC,&xtop,NULL,box,FALSE);
+
+ *nslicex= (int)(box[XX][XX]/binw + onehalf); //Calculate slicenr from binwidth
+ *nslicey= (int)(box[YY][YY]/binw + onehalf);
+ nslicez= (int)(box[ZZ][ZZ]/binw + onehalf);
+
+
+
+ ng = 1;
+ /* get index groups */
+ printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
+ snew(grpname,ng);
+ snew(index,ng);
+ snew(isize,ng);
+ get_index(&top.atoms,fnNDX,ng,isize,index,grpname);
+
+ /* Analyze trajectory */
+ natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box);
+ if ( natoms > top.atoms.nr )
+ gmx_fatal(FARGS,"Topology (%d atoms) does not match trajectory (%d atoms)",
+ top.atoms.nr,natoms);
+ check_index(NULL,ng,index[0],NULL,natoms);
+
+
+ /*Prepare structures for temporary storage of frame info*/
+ snew(sg_grid,*nslicex);
+ snew(sk_grid,*nslicex);
+ for(i=0;i<*nslicex;i++){
+ snew(sg_grid[i],*nslicey);
+ snew(sk_grid[i],*nslicey);
+ for(j=0;j<*nslicey;j++){
+ snew(sg_grid[i][j],nslicez);
+ snew(sk_grid[i][j],nslicez);
+ }
+ }
+
+ sg_4d=NULL;
+ sk_4d=NULL;
+ *nframes = 0;
+ framenr=0;
+
+/* Loop over frames*/
+ do {
+
+//Initialize box meshes (temporary storage for each tblock frame -reinitialise every tblock steps
+ if(framenr%tblock ==0){
+ srenew(sk_4d,*nframes+1);
+ srenew(sg_4d,*nframes+1);
+ snew(sg_fravg,*nslicex);
+ snew(sk_fravg,*nslicex);
+ for(i=0;i<*nslicex;i++){
+ snew(sg_fravg[i],*nslicey);
+ snew(sk_fravg[i],*nslicey);
+ for(j=0;j<*nslicey;j++){
+ snew(sg_fravg[i][j],nslicez);
+ snew(sk_fravg[i][j],nslicez);
+ }
+ }
+ }
+
+ find_tetra_order_grid(top,ePBC,natoms,box,x,isize[0],index[0],t,
+ &sg,&sk,*nslicex,*nslicey,nslicez,sg_grid,sk_grid);
+ for(i=0;i<*nslicex;i++){
+ for(j=0;j<*nslicey;j++){
+ for(k=0;k<nslicez;k++){
+ sk_fravg[i][j][k]+=sk_grid[i][j][k]/tblock;
+ sg_fravg[i][j][k]+=sg_grid[i][j][k]/tblock;
+ }
+ }
+ }
+
+ framenr++;
+
+ if(framenr%tblock== 0){
+ sk_4d[*nframes] = sk_fravg;
+ sg_4d[*nframes] = sg_fravg;
+ (*nframes)++;
+ }
+
+ } while (read_next_x(oenv,status,&t,natoms,x,box));
+close_trj(status);
+
+ sfree(grpname);
+ sfree(index);
+ sfree(isize);
+
+//Debugging for printing out the entire order parameter meshes.
+ if(debug){
+ fpsg = xvgropen(sgslfn,"S\\sg\\N Angle Order Parameter / Meshpoint","(nm)","S\\sg\\N",oenv);
+ fpsk = xvgropen(skslfn,"S\\sk\\N Distance Order Parameter / Meshpoint","(nm)","S\\sk\\N",oenv);
+ for(n=0;n<(*nframes);n++){
+ fprintf(fpsg,"%i\n",n);
+ fprintf(fpsk,"%i\n",n);
+ for(i=0; (i<*nslicex); i++) {
+ for(j=0;j<*nslicey;j++){
+ for(k=0;k<nslicez;k++){
+ 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]);
+ 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]);
+ }
+ }
+ }
+ }
+fclose(fpsg);
+fclose(fpsk);
+ }
+
+
+ /* Find positions of interface z by scanning orderparam for each frame and for each xy-mesh cylinder along z*/
+
+//Simple trial: assume interface is in the middle of -sgang1 and sgang2
+ sgintf=0.5*(sgang1+sgang2);
+
+
+//Allocate memory for interface arrays;
+snew((*intfpos),2);
+snew((*intfpos)[0],*nframes);
+snew((*intfpos)[1],*nframes);
+
+bins=(*nslicex)*(*nslicey);
+
+
+snew(perm,nslicez); //permutation array for sorting along normal coordinate
+
+
+ for (n=0;n<*nframes;n++){
+ snew((*intfpos)[0][n],bins);
+ snew((*intfpos)[1][n],bins);
+ for(i=0;i<*nslicex;i++){
+ for(j=0;j<*nslicey;j++){
+ rangeArray(perm,nslicez); //reset permutation array to identity
+ /*Binsearch returns 2 bin-numbers where the order param is <= setpoint sgintf*/
+ ndx1=start_binsearch(sg_4d[n][i][j],perm,0,nslicez/2-1,sgintf,1);
+ ndx2=start_binsearch(sg_4d[n][i][j],perm,nslicez/2,nslicez-1,sgintf,-1);
+ /*Use linear interpolation to smooth out the interface position*/
+
+ //left interface (0)
+ /*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){
+ pos=( (sgintf-sg_4d[n][i][j][perm[ndx1]])*perm[ndx1+1]+(sg_4d[n][i][j][perm[ndx1+1]]-sgintf)*perm[ndx1 ])*/
+ (*intfpos)[0][n][j+*nslicey*i]=(perm[ndx1]+onehalf)*binw;
+ //right interface (1)
+ //alpha=(sgintf-sg_4d[n][i][j][perm[ndx2]])/(sg_4d[n][i][j][perm[ndx2]+1]-sg_4d[n][i][j][perm[ndx2]]);
+ //(*intfpos)[1][n][j+*nslicey*i]=((1-alpha)*perm[ndx2]+alpha*(perm[ndx2]+1)+onehalf)*box[ZZ][ZZ]/nslicez;
+ (*intfpos)[1][n][j+*nslicey*i]=(perm[ndx2]+onehalf)*binw;
+ }
+ }
+ }
+
+
+ //sfree(perm);
+ sfree(sk_4d);
+ sfree(sg_4d);
+ //sfree(sg_grid);
+ //sfree(sk_grid);
+
+
+}
+
- void writeraw(real ***surf, int tblocks,int xbins, int ybins,char **fnms){
++static void writesurftoxpms(real ***surf, int tblocks,int xbins, int ybins, real bw,char **outfiles,int maplevels )
+{
+
+ char numbuf[8];
+ int n, i, j;
+ real **profile1, **profile2;
+ real max1, max2, min1, min2, *xticks, *yticks;
+ t_rgb lo={1,1,1};
+ t_rgb hi={0,0,0};
+ FILE *xpmfile1, *xpmfile2;
+
+/*Prepare xpm structures for output*/
+
+//Allocate memory to tick's and matrices
+snew (xticks,xbins+1);
+snew (yticks,ybins+1);
+
+profile1=mk_matrix(xbins,ybins,FALSE);
+profile2=mk_matrix(xbins,ybins,FALSE);
+
+for (i=0;i<xbins+1;i++) xticks[i]+=bw;
+for (j=0;j<ybins+1;j++) yticks[j]+=bw;
+
+xpmfile1 = ffopen(outfiles[0],"w");
+xpmfile2 = ffopen(outfiles[1],"w");
+
+max1=max2=0.0;
+min1=min2=1000.00;
+
+for(n=0;n<tblocks;n++){
+sprintf(numbuf,"%5d",n);
+//Filling matrices for inclusion in xpm-files
+ for(i=0;i<xbins;i++){
+ for(j=0;j<ybins;j++){
+ profile1[i][j]=(surf[0][n][j+ybins*i]);
+ profile2[i][j]=(surf[1][n][j+ybins*i]);
+ //Finding max and min values
+ if(profile1[i][j]>max1) max1=profile1[i][j];
+ if(profile1[i][j]<min1) min1=profile1[i][j];
+ if(profile2[i][j]>max2) max2=profile2[i][j];
+ if(profile2[i][j]<min2) min2=profile2[i][j];
+ }
+ }
+
+ write_xpm(xpmfile1,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile1,min1,max1,lo,hi,&maplevels);
+ write_xpm(xpmfile2,3,numbuf,"Height","x[nm]","y[nm]",xbins,ybins,xticks,yticks,profile2,min2,max2,lo,hi,&maplevels);
+}
+
+ffclose(xpmfile1);
+ffclose(xpmfile2);
+
+
+
+sfree(profile1);
+sfree(profile2);
+sfree(xticks);
+sfree(yticks);
+}
+
+
- static bool bFourier = FALSE;
- static bool bRawOut = FALSE;
++static void writeraw(real ***surf, int tblocks,int xbins, int ybins,char **fnms){
+ FILE *raw1, *raw2;
+ int i,j,n;
+
+ raw1=ffopen(fnms[0],"w");
+ raw2=ffopen(fnms[1],"w");
+ fprintf(raw1,"#Legend\n#TBlock\n#Xbin Ybin Z t\n");
+ fprintf(raw2,"#Legend\n#TBlock\n#Xbin Ybin Z t\n");
+ for (n=0;n<tblocks;n++){
+ fprintf(raw1,"%5d\n",n);
+ fprintf(raw2,"%5d\n",n);
+ for(i=0;i<xbins;i++){
+ for(j=0;j<ybins;j++){
+ fprintf(raw1,"%i %i %8.5f\n",i,j,(surf[0][n][j+ybins*i]));
+ fprintf(raw2,"%i %i %8.5f\n",i,j,(surf[1][n][j+ybins*i]));
+ }
+ }
+ }
+
+ ffclose(raw1);
+ ffclose(raw2);
+}
+
+
+
+int gmx_hydorder(int argc,char *argv[])
+{
+ static const char *desc[] = {
+ "The tetrahedrality order parameters can be determined",
+ "around an atom. Both angle an distance order parameters are calculated. See",
+ "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
+ "for more details.[BR]"
+ "This application calculates the orderparameter in a 3d-mesh in the box, and",
+ "with 2 phases in the box gives the user the option to define a 2D interface in time",
+ "separating the faces by specifying parameters -sgang1 and -sgang2 (It is important",
+ "to select these judiciously)"};
+
+ int axis = 0;
+ static int nsttblock=1;
+ static int nlevels=100;
+ static real binwidth = 1.0; /* binwidth in mesh */
+ static real sg1 = 1;
+ static real sg2 = 1; /* order parameters for bulk phases */
- int nfspect,nfxpm, nfraw;
- output_env_t oenv;
++ static gmx_bool bFourier = FALSE;
++ static gmx_bool bRawOut = FALSE;
+ int frames,xslices,yslices; /* Dimensions of interface arrays*/
+ real ***intfpos; /* Interface arrays (intfnr,t,xy) -potentially large */
+ static char *normal_axis[] = { NULL, "z", "x", "y", NULL };
+
+ t_pargs pa[] = {
+ { "-d", FALSE, etENUM, {normal_axis},
+ "Direction of the normal on the membrane" },
+ { "-bw", FALSE, etREAL, {&binwidth},
+ "Binwidth of box mesh" },
+ { "-sgang1",FALSE,etREAL, {&sg1},
+ "tetrahedral angle parameter in Phase 1 (bulk)" },
+ { "-sgang2",FALSE,etREAL, {&sg2},
+ "tetrahedral angle parameter in Phase 2 (bulk)" },
+ { "-tblock",FALSE,etINT,{&nsttblock},
+ "Number of frames in one time-block average"},
+ { "-nlevel", FALSE,etINT, {&nlevels},
+ "Number of Height levels in 2D - XPixMaps"}
+ };
+
+ t_filenm fnm[] = { /* files for g_order */
+ { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
+ { efNDX, "-n", NULL, ffREAD }, /* index file */
+ { efTPX, "-s", NULL, ffREAD }, /* topology file */
+ { efXPM, "-o", "intf", ffWRMULT}, /* XPM- surface maps */
+ { efOUT,"-or","raw", ffOPTWRMULT }, /* xvgr output file */
+ { efOUT,"-Spect","intfspect",ffOPTWRMULT}, /* Fourier spectrum interfaces */
+ };
+#define NFILE asize(fnm)
+
+ /*Filenames*/
+ const char *ndxfnm,*tpsfnm,*trxfnm;
+ char **spectra,**intfn, **raw;
- NFILE,fnm,asize(pa),pa,asize(desc),desc,0, NULL,&oenv);
++ int nfspect,nfxpm, nfraw;
++ output_env_t oenv;
+
+ CopyRight(stderr,argv[0]);
+
+ parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
- calc_tetra_order_interface(ndxfnm,tpsfnm,trxfnm,binwidth,nsttblock,&frames,&xslices,&yslices,
- sg1,sg2,&intfpos,oenv);
++ NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
+ bFourier= opt2bSet("-Spect",NFILE,fnm);
+ bRawOut = opt2bSet("-or",NFILE,fnm);
+
+ if (binwidth < 0.0)
+ gmx_fatal(FARGS,"Can not have binwidth < 0");
+
+ ndxfnm = ftp2fn(efNDX,NFILE,fnm);
+ tpsfnm = ftp2fn(efTPX,NFILE,fnm);
+ trxfnm = ftp2fn(efTRX,NFILE,fnm);
+
+ /* Calculate axis */
+ if (strcmp(normal_axis[0],"x") == 0) axis = XX;
+ else if (strcmp(normal_axis[0],"y") == 0) axis = YY;
+ else if (strcmp(normal_axis[0],"z") == 0) axis = ZZ;
+ else gmx_fatal(FARGS,"Invalid axis, use x, y or z");
+
+ switch (axis) {
+ case 0:
+ fprintf(stderr,"Taking x axis as normal to the membrane\n");
+ break;
+ case 1:
+ fprintf(stderr,"Taking y axis as normal to the membrane\n");
+ break;
+ case 2:
+ fprintf(stderr,"Taking z axis as normal to the membrane\n");
+ break;
+ }
+
+ /* tetraheder order parameter */
+ /* If either of the options is set we compute both */
+ nfxpm=opt2fns(&intfn,"-o",NFILE,fnm);
+ if(nfxpm!=2){
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfxpm);
+ }
++ calc_tetra_order_interface(ndxfnm,tpsfnm,trxfnm,binwidth,nsttblock,&frames,&xslices,&yslices,sg1,sg2,&intfpos,oenv);
+ writesurftoxpms(intfpos,frames,xslices,yslices,binwidth,intfn,nlevels);
+
+ if(bFourier){
+ nfspect=opt2fns(&spectra,"-Spect",NFILE,fnm);
+ if(nfspect!=2){
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfspect);
+ }
+ powerspectavg(intfpos, frames,xslices,yslices,spectra);
+ }
+
+ if (bRawOut){
+ nfraw=opt2fns(&raw,"-or",NFILE,fnm);
+ if(nfraw!=2){
+ gmx_fatal(FARGS,"No or not correct number (2) of output-files: %d",nfraw);
+ }
+ writeraw(intfpos,frames,xslices,yslices,raw);
+ }
+
+
+
+ thanx(stderr);
+
+ return 0;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff