2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
41 #include "gromacs/commandline/pargs.h"
42 #include "gromacs/correlationfunctions/autocorr.h"
43 #include "gromacs/correlationfunctions/expfit.h"
44 #include "gromacs/fileio/matio.h"
45 #include "gromacs/fileio/tpxio.h"
46 #include "gromacs/fileio/trxio.h"
47 #include "gromacs/fileio/xvgr.h"
48 #include "gromacs/gmxana/binsearch.h"
49 #include "gromacs/gmxana/dens_filter.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/gmxana/gstat.h"
52 #include "gromacs/gmxana/powerspect.h"
53 #include "gromacs/legacyheaders/copyrite.h"
54 #include "gromacs/legacyheaders/macros.h"
55 #include "gromacs/legacyheaders/typedefs.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/pbcutil/pbc.h"
59 #include "gromacs/pbcutil/rmpbc.h"
60 #include "gromacs/topology/index.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/exceptions.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/smalloc.h"
68 #define FLOOR(x) ((int) floor(x))
70 #define FLOOR(x) ((int) floorf(x))
74 methSEL
, methBISECT
, methFUNCFIT
, methNR
77 static void center_coords(t_atoms
*atoms
, matrix box
, rvec x0
[], int axis
)
81 rvec com
, shift
, box_center
;
85 for (i
= 0; (i
< atoms
->nr
); i
++)
87 mm
= atoms
->atom
[i
].m
;
89 for (m
= 0; (m
< DIM
); m
++)
91 com
[m
] += mm
*x0
[i
][m
];
94 for (m
= 0; (m
< DIM
); m
++)
98 calc_box_center(ecenterDEF
, box
, box_center
);
99 rvec_sub(box_center
, com
, shift
);
100 shift
[axis
] -= box_center
[axis
];
102 for (i
= 0; (i
< atoms
->nr
); i
++)
104 rvec_dec(x0
[i
], shift
);
109 static void density_in_time (const char *fn
, atom_id
**index
, int gnx
[], 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
)
113 * *****Densdevel pointer to array of density values in slices and frame-blocks Densdevel[*nsttblock][*xslices][*yslices][*zslices]
115 * nsttblock - nr of frames in each time-block
116 * bw widths of normal slices
118 * axis - axis direction (normal to slices)
119 * nndx - number ot atoms in **index
120 * grpn - group number in index
123 gmx_rmpbc_t gpbc
= NULL
;
124 matrix box
; /* Box - 3x3 -each step*/
125 rvec
*x0
; /* List of Coord without PBC*/
126 int i
, j
, /* loop indices, checks etc*/
127 ax1
= 0, ax2
= 0, /* tangent directions */
128 framenr
= 0, /* frame number in trajectory*/
129 slicex
, slicey
, slicez
; /*slice # of x y z position */
130 real
***Densslice
= NULL
; /* Density-slice in one frame*/
131 real dscale
; /*physical scaling factor*/
132 real t
, x
, y
, z
; /* time and coordinates*/
135 *tblock
= 0; /* blocknr in block average - initialise to 0*/
136 /* Axis: X=0, Y=1,Z=2 */
140 ax1
= YY
; ax2
= ZZ
; /*Surface: YZ*/
143 ax1
= ZZ
; ax2
= XX
; /* Surface : XZ*/
146 ax1
= XX
; ax2
= YY
; /* Surface XY*/
149 gmx_fatal(FARGS
, "Invalid axes. Terminating\n");
152 if (read_first_x(oenv
, &status
, fn
, &t
, &x0
, box
) == 0)
154 gmx_fatal(FARGS
, "Could not read coordinates from file"); /* Open trajectory for read*/
158 *zslices
= 1+FLOOR(box
[axis
][axis
]/bwz
);
159 *yslices
= 1+FLOOR(box
[ax2
][ax2
]/bw
);
160 *xslices
= 1+FLOOR(box
[ax1
][ax1
]/bw
);
163 if (*xslices
< *yslices
)
173 "\nDividing the box in %5d x %5d x %5d slices with binw %f along axis %d\n", *xslices
, *yslices
, *zslices
, bw
, axis
);
176 /****Start trajectory processing***/
178 /*Initialize Densdevel and PBC-remove*/
179 gpbc
= gmx_rmpbc_init(&top
->idef
, ePBC
, top
->atoms
.nr
);
185 bbww
[XX
] = box
[ax1
][ax1
]/ *xslices
;
186 bbww
[YY
] = box
[ax2
][ax2
]/ *yslices
;
187 bbww
[ZZ
] = box
[axis
][axis
]/ *zslices
;
188 gmx_rmpbc(gpbc
, top
->atoms
.nr
, box
, x0
);
189 /*Reset Densslice every nsttblock steps*/
190 /* The first conditional is for clang to understand that this branch is
191 * always taken the first time. */
192 if (Densslice
== NULL
|| framenr
% nsttblock
== 0)
194 snew(Densslice
, *xslices
);
195 for (i
= 0; i
< *xslices
; i
++)
197 snew(Densslice
[i
], *yslices
);
198 for (j
= 0; j
< *yslices
; j
++)
200 snew(Densslice
[i
][j
], *zslices
);
204 /* Allocate Memory to extra frame in Densdevel - rather stupid approach:
205 * A single frame each time, although only every nsttblock steps.
207 srenew(*Densdevel
, *tblock
+1);
208 (*Densdevel
)[*tblock
] = Densslice
;
211 dscale
= (*xslices
)*(*yslices
)*(*zslices
)*AMU
/ (box
[ax1
][ax1
]*box
[ax2
][ax2
]*box
[axis
][axis
]*nsttblock
*(NANO
*NANO
*NANO
));
215 center_coords(&top
->atoms
, box
, x0
, axis
);
219 for (j
= 0; j
< gnx
[0]; j
++)
220 { /*Loop over all atoms in selected index*/
221 x
= x0
[index
[0][j
]][ax1
];
222 y
= x0
[index
[0][j
]][ax2
];
223 z
= x0
[index
[0][j
]][axis
];
228 while (x
> box
[ax1
][ax1
])
237 while (y
> box
[ax2
][ax2
])
244 z
+= box
[axis
][axis
];
246 while (z
> box
[axis
][axis
])
248 z
-= box
[axis
][axis
];
251 slicex
= ((int) (x
/bbww
[XX
])) % *xslices
;
252 slicey
= ((int) (y
/bbww
[YY
])) % *yslices
;
253 slicez
= ((int) (z
/bbww
[ZZ
])) % *zslices
;
254 Densslice
[slicex
][slicey
][slicez
] += (top
->atoms
.atom
[index
[0][j
]].m
*dscale
);
259 if (framenr
% nsttblock
== 0)
261 /*Implicit incrementation of Densdevel via renewal of Densslice*/
262 /*only every nsttblock steps*/
267 while (read_next_x(oenv
, status
, &t
, x0
, box
));
270 /*Free memory we no longer need and exit.*/
271 gmx_rmpbc_done(gpbc
);
277 fp
= fopen("koko.xvg", "w");
278 for (j
= 0; (j
< *zslices
); j
++)
280 fprintf(fp
, "%5d", j
);
281 for (i
= 0; (i
< *tblock
); i
++)
283 fprintf(fp
, " %10g", (*Densdevel
)[i
][9][1][j
]);
292 static void outputfield(const char *fldfn
, real
****Densmap
,
293 int xslices
, int yslices
, int zslices
, int tdim
)
295 /*Debug-filename and filehandle*/
306 fldH
= gmx_ffopen(fldfn
, "w");
307 fwrite(dim
, sizeof(int), 4, fldH
);
308 for (n
= 0; n
< tdim
; n
++)
310 for (i
= 0; i
< xslices
; i
++)
312 for (j
= 0; j
< yslices
; j
++)
314 for (k
= 0; k
< zslices
; k
++)
316 fwrite(&(Densmap
[n
][i
][j
][k
]), sizeof(real
), 1, fldH
);
317 totdens
+= (Densmap
[n
][i
][j
][k
]);
322 totdens
/= (xslices
*yslices
*zslices
*tdim
);
323 fprintf(stderr
, "Total density [kg/m^3] %8f", totdens
);
327 static void filterdensmap(real
****Densmap
, int xslices
, int yslices
, int zslices
, int tblocks
, int ftsize
)
333 std
= ((real
)order
/2.0);
335 snew(kernel
, ftsize
);
336 gausskernel(kernel
, ftsize
, var
);
337 for (n
= 0; n
< tblocks
; n
++)
339 for (i
= 0; i
< xslices
; i
++)
341 for (j
= 0; j
< yslices
; j
++)
343 periodic_convolution(zslices
, Densmap
[n
][i
][j
], ftsize
, kernel
);
352 static void interfaces_txy (real
****Densmap
, int xslices
, int yslices
, int zslices
,
353 int tblocks
, real binwidth
, int method
,
354 real dens1
, real dens2
, t_interf
****intf1
,
355 t_interf
****intf2
, const output_env_t oenv
)
357 /*Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)*/
359 real
*zDensavg
; /* zDensavg[z]*/
362 int ndx1
, ndx2
, *zperm
;
364 real splitpoint
, startpoint
, endpoint
;
365 real
*sigma1
, *sigma2
;
368 double *fit1
= NULL
, *fit2
= NULL
;
369 const double *avgfit1
;
370 const double *avgfit2
;
371 const real onehalf
= 1.00/2.00;
372 t_interf
***int1
= NULL
, ***int2
= NULL
; /*Interface matrices [t][x,y] - last index in row-major order*/
373 /*Create int1(t,xy) and int2(t,xy) arrays with correct number of interf_t elements*/
374 xysize
= xslices
*yslices
;
377 for (i
= 0; i
< tblocks
; i
++)
379 snew(int1
[i
], xysize
);
380 snew(int2
[i
], xysize
);
381 for (j
= 0; j
< xysize
; j
++)
385 init_interf(int1
[i
][j
]);
386 init_interf(int2
[i
][j
]);
390 if (method
== methBISECT
)
392 densmid
= onehalf
*(dens1
+dens2
);
393 snew(zperm
, zslices
);
394 for (n
= 0; n
< tblocks
; n
++)
396 for (i
= 0; i
< xslices
; i
++)
398 for (j
= 0; j
< yslices
; j
++)
400 rangeArray(zperm
, zslices
); /*reset permutation array to identity*/
401 /*Binsearch returns slice-nr where the order param is <= setpoint sgmid*/
402 ndx1
= start_binsearch(Densmap
[n
][i
][j
], zperm
, 0, zslices
/2-1, densmid
, 1);
403 ndx2
= start_binsearch(Densmap
[n
][i
][j
], zperm
, zslices
/2, zslices
-1, densmid
, -1);
405 /* Linear interpolation (for use later if time allows)
406 * rho_1s= Densmap[n][i][j][zperm[ndx1]]
407 * rho_1e =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
408 * rho_2s =Densmap[n][i][j][zperm[ndx2+1]]
409 * rho_2e =Densmap[n][i][j][zperm[ndx2]]
410 * For 1st interface we have:
411 densl= Densmap[n][i][j][zperm[ndx1]];
412 densr= Densmap[n][i][j][zperm[ndx1+1]];
413 alpha=(densmid-densl)/(densr-densl);
414 deltandx=zperm[ndx1+1]-zperm[ndx1];
417 printf("Alpha, Deltandx %f %i\n", alpha,deltandx);
419 if(abs(alpha)>1.0 || abs(deltandx)>3){
424 pos=zperm[ndx1]+alpha*deltandx;
425 spread=binwidth*deltandx;
427 * For the 2nd interface can use the same formulation, since alpha should become negative ie:
428 * alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/(Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
429 * deltandx=zperm[ndx2+1]-zperm[ndx2];
430 * pos=zperm[ndx2]+alpha*deltandx; */
432 /*After filtering we use the direct approach */
433 int1
[n
][j
+(i
*yslices
)]->Z
= (zperm
[ndx1
]+onehalf
)*binwidth
;
434 int1
[n
][j
+(i
*yslices
)]->t
= binwidth
;
435 int2
[n
][j
+(i
*yslices
)]->Z
= (zperm
[ndx2
]+onehalf
)*binwidth
;
436 int2
[n
][j
+(i
*yslices
)]->t
= binwidth
;
442 if (method
== methFUNCFIT
)
444 /*Assume a box divided in 2 along midpoint of z for starters*/
446 endpoint
= binwidth
*zslices
;
447 splitpoint
= (startpoint
+endpoint
)/2.0;
448 /*Initial fit proposals*/
449 beginfit1
[0] = dens1
;
450 beginfit1
[1] = dens2
;
451 beginfit1
[2] = (splitpoint
/2);
454 beginfit2
[0] = dens2
;
455 beginfit2
[1] = dens1
;
456 beginfit2
[2] = (3*splitpoint
/2);
459 snew(zDensavg
, zslices
);
460 snew(sigma1
, zslices
);
461 snew(sigma2
, zslices
);
463 for (k
= 0; k
< zslices
; k
++)
465 sigma1
[k
] = sigma2
[k
] = 1;
467 /*Calculate average density along z - avoid smoothing by using coarse-grained-mesh*/
468 for (k
= 0; k
< zslices
; k
++)
470 for (n
= 0; n
< tblocks
; n
++)
472 for (i
= 0; i
< xslices
; i
++)
474 for (j
= 0; j
< yslices
; j
++)
476 zDensavg
[k
] += (Densmap
[n
][i
][j
][k
]/(xslices
*yslices
*tblocks
));
484 xvg
= xvgropen("DensprofileonZ.xvg", "Averaged Densityprofile on Z", "z[nm]", "Density[kg/m^3]", oenv
);
485 for (k
= 0; k
< zslices
; k
++)
487 fprintf(xvg
, "%4f.3 %8f.4\n", k
*binwidth
, zDensavg
[k
]);
492 /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
494 /*Fit 1st half of box*/
495 do_lmfit(zslices
, zDensavg
, sigma1
, binwidth
, NULL
, startpoint
, splitpoint
, oenv
, FALSE
, effnERF
, beginfit1
, 8, NULL
);
496 /*Fit 2nd half of box*/
497 do_lmfit(zslices
, zDensavg
, sigma2
, binwidth
, NULL
, splitpoint
, endpoint
, oenv
, FALSE
, effnERF
, beginfit2
, 8, NULL
);
499 /*Initialise the const arrays for storing the average fit parameters*/
505 /*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*/
506 for (n
= 0; n
< tblocks
; n
++)
508 for (i
= 0; i
< xslices
; i
++)
510 for (j
= 0; j
< yslices
; j
++)
512 /*Reinitialise fit for each mesh-point*/
515 for (k
= 0; k
< 4; k
++)
517 fit1
[k
] = avgfit1
[k
];
518 fit2
[k
] = avgfit2
[k
];
520 /*Now fit and store in structures in row-major order int[n][i][j]*/
521 do_lmfit(zslices
, Densmap
[n
][i
][j
], sigma1
, binwidth
, NULL
, startpoint
, splitpoint
, oenv
, FALSE
, effnERF
, fit1
, 0, NULL
);
522 int1
[n
][j
+(yslices
*i
)]->Z
= fit1
[2];
523 int1
[n
][j
+(yslices
*i
)]->t
= fit1
[3];
524 do_lmfit(zslices
, Densmap
[n
][i
][j
], sigma2
, binwidth
, NULL
, splitpoint
, endpoint
, oenv
, FALSE
, effnERF
, fit2
, 0, NULL
);
525 int2
[n
][j
+(yslices
*i
)]->Z
= fit2
[2];
526 int2
[n
][j
+(yslices
*i
)]->t
= fit2
[3];
538 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
)
542 real
**profile1
, **profile2
;
543 real max1
, max2
, min1
, min2
, *xticks
, *yticks
;
544 t_rgb lo
= {0, 0, 0};
545 t_rgb hi
= {1, 1, 1};
546 FILE *xpmfile1
, *xpmfile2
;
548 /*Prepare xpm structures for output*/
550 /*Allocate memory to tick's and matrices*/
551 snew (xticks
, xbins
+1);
552 snew (yticks
, ybins
+1);
554 profile1
= mk_matrix(xbins
, ybins
, FALSE
);
555 profile2
= mk_matrix(xbins
, ybins
, FALSE
);
557 for (i
= 0; i
< xbins
+1; i
++)
561 for (j
= 0; j
< ybins
+1; j
++)
566 xpmfile1
= gmx_ffopen(outfiles
[0], "w");
567 xpmfile2
= gmx_ffopen(outfiles
[1], "w");
570 min1
= min2
= zbins
*bwz
;
572 for (n
= 0; n
< tblocks
; n
++)
574 sprintf(numbuf
, "tblock: %4i", n
);
575 /*Filling matrices for inclusion in xpm-files*/
576 for (i
= 0; i
< xbins
; i
++)
578 for (j
= 0; j
< ybins
; j
++)
580 profile1
[i
][j
] = (surf1
[n
][j
+ybins
*i
])->Z
;
581 profile2
[i
][j
] = (surf2
[n
][j
+ybins
*i
])->Z
;
582 /*Finding max and min values*/
583 if (profile1
[i
][j
] > max1
)
585 max1
= profile1
[i
][j
];
587 if (profile1
[i
][j
] < min1
)
589 min1
= profile1
[i
][j
];
591 if (profile2
[i
][j
] > max2
)
593 max2
= profile2
[i
][j
];
595 if (profile2
[i
][j
] < min2
)
597 min2
= profile2
[i
][j
];
602 write_xpm(xpmfile1
, 3, numbuf
, "Height", "x[nm]", "y[nm]", xbins
, ybins
, xticks
, yticks
, profile1
, min1
, max1
, lo
, hi
, &maplevels
);
603 write_xpm(xpmfile2
, 3, numbuf
, "Height", "x[nm]", "y[nm]", xbins
, ybins
, xticks
, yticks
, profile2
, min2
, max2
, lo
, hi
, &maplevels
);
606 gmx_ffclose(xpmfile1
);
607 gmx_ffclose(xpmfile2
);
616 static void writeraw(t_interf
***int1
, t_interf
***int2
, int tblocks
,
617 int xbins
, int ybins
, char **fnms
,
618 const output_env_t oenv
)
623 raw1
= gmx_ffopen(fnms
[0], "w");
624 raw2
= gmx_ffopen(fnms
[1], "w");
627 gmx::BinaryInformationSettings settings
;
628 settings
.generatedByHeader(true);
629 settings
.linePrefix("# ");
630 gmx::printBinaryInformation(raw1
, output_env_get_program_context(oenv
),
632 gmx::printBinaryInformation(raw2
, output_env_get_program_context(oenv
),
635 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
636 fprintf(raw1
, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
637 fprintf(raw2
, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
638 fprintf(raw1
, "%i %i %i\n", tblocks
, xbins
, ybins
);
639 fprintf(raw2
, "%i %i %i\n", tblocks
, xbins
, ybins
);
640 for (n
= 0; n
< tblocks
; n
++)
642 for (i
= 0; i
< xbins
; i
++)
644 for (j
= 0; j
< ybins
; j
++)
646 fprintf(raw1
, "%i %i %8.5f %6.4f\n", i
, j
, (int1
[n
][j
+ybins
*i
])->Z
, (int1
[n
][j
+ybins
*i
])->t
);
647 fprintf(raw2
, "%i %i %8.5f %6.4f\n", i
, j
, (int2
[n
][j
+ybins
*i
])->Z
, (int2
[n
][j
+ybins
*i
])->t
);
658 int gmx_densorder(int argc
, char *argv
[])
660 static const char *desc
[] = {
661 "[THISMODULE] reduces a two-phase density distribution",
662 "along an axis, computed over a MD trajectory,",
663 "to 2D surfaces fluctuating in time, by a fit to",
664 "a functional profile for interfacial densities.",
665 "A time-averaged spatial representation of the",
666 "interfaces can be output with the option [TT]-tavg[tt]."
669 /* Extra arguments - but note how you always get the begin/end
670 * options when running the program, without mentioning them here!
677 static real binw
= 0.2;
678 static real binwz
= 0.05;
679 static real dens1
= 0.00;
680 static real dens2
= 1000.00;
681 static int ftorder
= 0;
682 static int nsttblock
= 100;
684 static const char *axtitle
= "Z";
685 atom_id
**index
; /* Index list for single group*/
686 int xslices
, yslices
, zslices
, tblock
;
687 static gmx_bool bGraph
= FALSE
;
688 static gmx_bool bCenter
= FALSE
;
689 static gmx_bool bFourier
= FALSE
;
690 static gmx_bool bRawOut
= FALSE
;
691 static gmx_bool bOut
= FALSE
;
692 static gmx_bool b1d
= FALSE
;
693 static int nlevels
= 100;
694 /*Densitymap - Densmap[t][x][y][z]*/
695 real
****Densmap
= NULL
;
696 /* Surfaces surf[t][surf_x,surf_y]*/
697 t_interf
***surf1
, ***surf2
;
699 static const char *meth
[] = {NULL
, "bisect", "functional", NULL
};
702 char **graphfiles
, **rawfiles
, **spectra
; /* Filenames for xpm-surface maps, rawdata and powerspectra */
703 int nfxpm
= -1, nfraw
, nfspect
; /* # files for interface maps and spectra = # interfaces */
706 { "-1d", FALSE
, etBOOL
, {&b1d
},
707 "Pseudo-1d interface geometry"},
708 { "-bw", FALSE
, etREAL
, {&binw
},
709 "Binwidth of density distribution tangential to interface"},
710 { "-bwn", FALSE
, etREAL
, {&binwz
},
711 "Binwidth of density distribution normal to interface"},
712 { "-order", FALSE
, etINT
, {&ftorder
},
713 "Order of Gaussian filter, order 0 equates to NO filtering"},
714 {"-axis", FALSE
, etSTR
, {&axtitle
},
715 "Axis Direction - X, Y or Z"},
716 {"-method", FALSE
, etENUM
, {meth
},
717 "Interface location method"},
718 {"-d1", FALSE
, etREAL
, {&dens1
},
719 "Bulk density phase 1 (at small z)"},
720 {"-d2", FALSE
, etREAL
, {&dens2
},
721 "Bulk density phase 2 (at large z)"},
722 { "-tblock", FALSE
, etINT
, {&nsttblock
},
723 "Number of frames in one time-block average"},
724 { "-nlevel", FALSE
, etINT
, {&nlevels
},
725 "Number of Height levels in 2D - XPixMaps"}
730 { efTPR
, "-s", NULL
, ffREAD
}, /* this is for the topology */
731 { efTRX
, "-f", NULL
, ffREAD
}, /* and this for the trajectory */
732 { efNDX
, "-n", NULL
, ffREAD
}, /* this is to select groups */
733 { efDAT
, "-o", "Density4D", ffOPTWR
}, /* This is for outputting the entire 4D densityfield in binary format */
734 { efOUT
, "-or", NULL
, ffOPTWRMULT
}, /* This is for writing out the entire information in the t_interf arrays */
735 { efXPM
, "-og", "interface", ffOPTWRMULT
}, /* This is for writing out the interface meshes - one xpm-file per tblock*/
736 { efOUT
, "-Spect", "intfspect", ffOPTWRMULT
}, /* This is for the trajectory averaged Fourier-spectra*/
739 #define NFILE asize(fnm)
741 /* This is the routine responsible for adding default options,
742 * calling the X/motif interface, etc. */
743 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_CAN_VIEW
,
744 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, 0, NULL
, &oenv
))
751 bFourier
= opt2bSet("-Spect", NFILE
, fnm
);
752 bRawOut
= opt2bSet("-or", NFILE
, fnm
);
753 bGraph
= opt2bSet("-og", NFILE
, fnm
);
754 bOut
= opt2bSet("-o", NFILE
, fnm
);
755 top
= read_top(ftp2fn(efTPR
, NFILE
, fnm
), &ePBC
);
761 axis
= toupper(axtitle
[0]) - 'X';
763 get_index(&top
->atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1, ngx
, index
, grpname
);
765 density_in_time(ftp2fn(efTRX
, NFILE
, fnm
), index
, ngx
, binw
, binwz
, nsttblock
, &Densmap
, &xslices
, &yslices
, &zslices
, &tblock
, top
, ePBC
, axis
, bCenter
, b1d
, oenv
);
769 filterdensmap(Densmap
, xslices
, yslices
, zslices
, tblock
, 2*ftorder
+1);
774 outputfield(opt2fn("-o", NFILE
, fnm
), Densmap
, xslices
, yslices
, zslices
, tblock
);
777 interfaces_txy(Densmap
, xslices
, yslices
, zslices
, tblock
, binwz
, eMeth
, dens1
, dens2
, &surf1
, &surf2
, oenv
);
782 /*Output surface-xpms*/
783 nfxpm
= opt2fns(&graphfiles
, "-og", NFILE
, fnm
);
786 gmx_fatal(FARGS
, "No or not correct number (2) of output-files: %d", nfxpm
);
788 writesurftoxpms(surf1
, surf2
, tblock
, xslices
, yslices
, zslices
, binw
, binwz
, graphfiles
, zslices
);
798 nfraw
= opt2fns(&rawfiles
, "-or", NFILE
, fnm
);
801 gmx_fatal(FARGS
, "No or not correct number (2) of output-files: %d", nfxpm
);
803 writeraw(surf1
, surf2
, tblock
, xslices
, yslices
, rawfiles
, oenv
);
810 nfspect
= opt2fns(&spectra
, "-Spect", NFILE
, fnm
);
813 gmx_fatal(FARGS
, "No or not correct number (2) of output-file-series: %d",
816 powerspectavg_intf(surf1
, surf2
, tblock
, xslices
, yslices
, spectra
);
820 if (bGraph
|| bFourier
|| bRawOut
)