Add replacements for pbc enumerations
[gromacs.git] / src / gromacs / gmxana / gmx_densorder.cpp
blob61803865554e8f803b38c24eb3178ca7ce1eb1c0
1 /*
2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,2018,2019, 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.
35 #include "gmxpre.h"
37 #include <cctype>
38 #include <cmath>
39 #include <cstring>
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/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/pbcutil/rmpbc.h"
57 #include "gromacs/topology/index.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/utility/arraysize.h"
60 #include "gromacs/utility/binaryinformation.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 enum {
69 methSEL, methBISECT, methFUNCFIT, methNR
72 static void center_coords(const t_atoms *atoms, matrix box, rvec x0[], int axis)
74 int i, m;
75 real tmass, mm;
76 rvec com, shift, box_center;
78 tmass = 0;
79 clear_rvec(com);
80 for (i = 0; (i < atoms->nr); i++)
82 mm = atoms->atom[i].m;
83 tmass += mm;
84 for (m = 0; (m < DIM); m++)
86 com[m] += mm*x0[i][m];
89 for (m = 0; (m < DIM); m++)
91 com[m] /= tmass;
93 calc_box_center(ecenterDEF, box, box_center);
94 rvec_sub(box_center, com, shift);
95 shift[axis] -= box_center[axis];
97 for (i = 0; (i < atoms->nr); i++)
99 rvec_dec(x0[i], shift);
104 static void density_in_time (const char *fn, int **index, const int gnx[], real bw, real bwz, int nsttblock, real *****Densdevel, int *xslices, int *yslices, int *zslices, int *tblock, const t_topology *top, int ePBC, int axis, gmx_bool bCenter, gmx_bool bps1d, const gmx_output_env_t *oenv)
108 * *****Densdevel pointer to array of density values in slices and frame-blocks Densdevel[*nsttblock][*xslices][*yslices][*zslices]
109 * Densslice[x][y][z]
110 * nsttblock - nr of frames in each time-block
111 * bw widths of normal slices
113 * axis - axis direction (normal to slices)
114 * nndx - number ot atoms in **index
115 * grpn - group number in index
117 t_trxstatus *status;
118 gmx_rmpbc_t gpbc = nullptr;
119 matrix box; /* Box - 3x3 -each step*/
120 rvec *x0; /* List of Coord without PBC*/
121 int i, j, /* loop indices, checks etc*/
122 ax1 = 0, ax2 = 0, /* tangent directions */
123 framenr = 0, /* frame number in trajectory*/
124 slicex, slicey, slicez; /*slice # of x y z position */
125 real ***Densslice = nullptr; /* Density-slice in one frame*/
126 real dscale; /*physical scaling factor*/
127 real t, x, y, z; /* time and coordinates*/
128 rvec bbww;
130 *tblock = 0; /* blocknr in block average - initialise to 0*/
131 /* Axis: X=0, Y=1,Z=2 */
132 switch (axis)
134 case 0:
135 ax1 = YY; ax2 = ZZ; /*Surface: YZ*/
136 break;
137 case 1:
138 ax1 = ZZ; ax2 = XX; /* Surface : XZ*/
139 break;
140 case 2:
141 ax1 = XX; ax2 = YY; /* Surface XY*/
142 break;
143 default:
144 gmx_fatal(FARGS, "Invalid axes. Terminating\n");
147 if (read_first_x(oenv, &status, fn, &t, &x0, box) == 0)
149 gmx_fatal(FARGS, "Could not read coordinates from file"); /* Open trajectory for read*/
153 *zslices = 1+static_cast<int>(std::floor(box[axis][axis]/bwz));
154 *yslices = 1+static_cast<int>(std::floor(box[ax2][ax2]/bw));
155 *xslices = 1+static_cast<int>(std::floor(box[ax1][ax1]/bw));
156 if (bps1d)
158 if (*xslices < *yslices)
160 *xslices = 1;
162 else
164 *yslices = 1;
167 fprintf(stderr,
168 "\nDividing the box in %5d x %5d x %5d slices with binw %f along axis %d\n", *xslices, *yslices, *zslices, bw, axis );
171 /****Start trajectory processing***/
173 /*Initialize Densdevel and PBC-remove*/
174 gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr);
176 *Densdevel = nullptr;
180 bbww[XX] = box[ax1][ax1]/ *xslices;
181 bbww[YY] = box[ax2][ax2]/ *yslices;
182 bbww[ZZ] = box[axis][axis]/ *zslices;
183 gmx_rmpbc(gpbc, top->atoms.nr, box, x0);
184 /*Reset Densslice every nsttblock steps*/
185 /* The first conditional is for clang to understand that this branch is
186 * always taken the first time. */
187 if (Densslice == nullptr || framenr % nsttblock == 0)
189 snew(Densslice, *xslices);
190 for (i = 0; i < *xslices; i++)
192 snew(Densslice[i], *yslices);
193 for (j = 0; j < *yslices; j++)
195 snew(Densslice[i][j], *zslices);
199 /* Allocate Memory to extra frame in Densdevel - rather stupid approach:
200 * A single frame each time, although only every nsttblock steps.
202 srenew(*Densdevel, *tblock+1);
203 (*Densdevel)[*tblock] = Densslice;
206 dscale = (*xslices)*(*yslices)*(*zslices)*AMU/ (box[ax1][ax1]*box[ax2][ax2]*box[axis][axis]*nsttblock*(NANO*NANO*NANO));
208 if (bCenter)
210 center_coords(&top->atoms, box, x0, axis);
214 for (j = 0; j < gnx[0]; j++)
215 { /*Loop over all atoms in selected index*/
216 x = x0[index[0][j]][ax1];
217 y = x0[index[0][j]][ax2];
218 z = x0[index[0][j]][axis];
219 while (x < 0)
221 x += box[ax1][ax1];
223 while (x > box[ax1][ax1])
225 x -= box[ax1][ax1];
228 while (y < 0)
230 y += box[ax2][ax2];
232 while (y > box[ax2][ax2])
234 y -= box[ax2][ax2];
237 while (z < 0)
239 z += box[axis][axis];
241 while (z > box[axis][axis])
243 z -= box[axis][axis];
246 slicex = static_cast<int>(x/bbww[XX]) % *xslices;
247 slicey = static_cast<int>(y/bbww[YY]) % *yslices;
248 slicez = static_cast<int>(z/bbww[ZZ]) % *zslices;
249 Densslice[slicex][slicey][slicez] += (top->atoms.atom[index[0][j]].m*dscale);
252 framenr++;
254 if (framenr % nsttblock == 0)
256 /*Implicit incrementation of Densdevel via renewal of Densslice*/
257 /*only every nsttblock steps*/
258 (*tblock)++;
262 while (read_next_x(oenv, status, &t, x0, box));
265 /*Free memory we no longer need and exit.*/
266 gmx_rmpbc_done(gpbc);
267 close_trx(status);
269 if (/* DISABLES CODE */ (false))
271 FILE *fp;
272 fp = fopen("koko.xvg", "w");
273 for (j = 0; (j < *zslices); j++)
275 fprintf(fp, "%5d", j);
276 for (i = 0; (i < *tblock); i++)
278 fprintf(fp, " %10g", (*Densdevel)[i][9][1][j]);
280 fprintf(fp, "\n");
282 fclose(fp);
287 static void outputfield(const char *fldfn, real ****Densmap,
288 int xslices, int yslices, int zslices, int tdim)
290 /*Debug-filename and filehandle*/
291 FILE *fldH;
292 int n, i, j, k;
293 int dim[4];
294 real totdens = 0;
296 dim[0] = tdim;
297 dim[1] = xslices;
298 dim[2] = yslices;
299 dim[3] = zslices;
301 fldH = gmx_ffopen(fldfn, "w");
302 fwrite(dim, sizeof(int), 4, fldH);
303 for (n = 0; n < tdim; n++)
305 for (i = 0; i < xslices; i++)
307 for (j = 0; j < yslices; j++)
309 for (k = 0; k < zslices; k++)
311 fwrite(&(Densmap[n][i][j][k]), sizeof(real), 1, fldH);
312 totdens += (Densmap[n][i][j][k]);
317 totdens /= (xslices*yslices*zslices*tdim);
318 fprintf(stderr, "Total density [kg/m^3] %8f", totdens);
319 gmx_ffclose(fldH);
322 static void filterdensmap(real ****Densmap, int xslices, int yslices, int zslices, int tblocks, int ftsize)
324 real *kernel;
325 real std, var;
326 int i, j, n, order;
327 order = ftsize/2;
328 std = order/2.0;
329 var = std*std;
330 snew(kernel, ftsize);
331 gausskernel(kernel, ftsize, var);
332 for (n = 0; n < tblocks; n++)
334 for (i = 0; i < xslices; i++)
336 for (j = 0; j < yslices; j++)
338 periodic_convolution(zslices, Densmap[n][i][j], ftsize, kernel);
347 static void interfaces_txy (real ****Densmap, int xslices, int yslices, int zslices,
348 int tblocks, real binwidth, int method,
349 real dens1, real dens2, t_interf ****intf1,
350 t_interf ****intf2, const gmx_output_env_t *oenv)
352 /*Returns two pointers to 3D arrays of t_interf structs containing (position,thickness) of the interface(s)*/
353 FILE *xvg;
354 real *zDensavg; /* zDensavg[z]*/
355 int i, j, k, n;
356 int xysize;
357 int ndx1, ndx2, *zperm;
358 real densmid;
359 real splitpoint, startpoint, endpoint;
360 real *sigma1, *sigma2;
361 double beginfit1[4];
362 double beginfit2[4];
363 double *fit1 = nullptr, *fit2 = nullptr;
364 const double *avgfit1;
365 const double *avgfit2;
366 const real onehalf = 1.00/2.00;
367 t_interf ***int1 = nullptr, ***int2 = nullptr; /*Interface matrices [t][x,y] - last index in row-major order*/
368 /*Create int1(t,xy) and int2(t,xy) arrays with correct number of interf_t elements*/
369 xysize = xslices*yslices;
370 snew(int1, tblocks);
371 snew(int2, tblocks);
372 for (i = 0; i < tblocks; i++)
374 snew(int1[i], xysize);
375 snew(int2[i], xysize);
376 for (j = 0; j < xysize; j++)
378 snew(int1[i][j], 1);
379 snew(int2[i][j], 1);
380 init_interf(int1[i][j]);
381 init_interf(int2[i][j]);
385 if (method == methBISECT)
387 densmid = onehalf*(dens1+dens2);
388 snew(zperm, zslices);
389 for (n = 0; n < tblocks; n++)
391 for (i = 0; i < xslices; i++)
393 for (j = 0; j < yslices; j++)
395 rangeArray(zperm, zslices); /*reset permutation array to identity*/
396 /*Binsearch returns slice-nr where the order param is <= setpoint sgmid*/
397 ndx1 = start_binsearch(Densmap[n][i][j], zperm, 0, zslices/2-1, densmid, 1);
398 ndx2 = start_binsearch(Densmap[n][i][j], zperm, zslices/2, zslices-1, densmid, -1);
400 /* Linear interpolation (for use later if time allows)
401 * rho_1s= Densmap[n][i][j][zperm[ndx1]]
402 * rho_1e =Densmap[n][i][j][zperm[ndx1+1]] - in worst case might be far off
403 * rho_2s =Densmap[n][i][j][zperm[ndx2+1]]
404 * rho_2e =Densmap[n][i][j][zperm[ndx2]]
405 * For 1st interface we have:
406 densl= Densmap[n][i][j][zperm[ndx1]];
407 densr= Densmap[n][i][j][zperm[ndx1+1]];
408 alpha=(densmid-densl)/(densr-densl);
409 deltandx=zperm[ndx1+1]-zperm[ndx1];
411 if(debug){
412 printf("Alpha, Deltandx %f %i\n", alpha,deltandx);
414 if(abs(alpha)>1.0 || abs(deltandx)>3){
415 pos=zperm[ndx1];
416 spread=-1;
418 else {
419 pos=zperm[ndx1]+alpha*deltandx;
420 spread=binwidth*deltandx;
422 * For the 2nd interface can use the same formulation, since alpha should become negative ie:
423 * alpha=(densmid-Densmap[n][i][j][zperm[ndx2]])/(Densmap[n][i][j][zperm[nxd2+1]]-Densmap[n][i][j][zperm[ndx2]]);
424 * deltandx=zperm[ndx2+1]-zperm[ndx2];
425 * pos=zperm[ndx2]+alpha*deltandx; */
427 /*After filtering we use the direct approach */
428 int1[n][j+(i*yslices)]->Z = (zperm[ndx1]+onehalf)*binwidth;
429 int1[n][j+(i*yslices)]->t = binwidth;
430 int2[n][j+(i*yslices)]->Z = (zperm[ndx2]+onehalf)*binwidth;
431 int2[n][j+(i*yslices)]->t = binwidth;
437 if (method == methFUNCFIT)
439 /*Assume a box divided in 2 along midpoint of z for starters*/
440 startpoint = 0.0;
441 endpoint = binwidth*zslices;
442 splitpoint = (startpoint+endpoint)/2.0;
443 /*Initial fit proposals*/
444 beginfit1[0] = dens1;
445 beginfit1[1] = dens2;
446 beginfit1[2] = (splitpoint/2);
447 beginfit1[3] = 0.5;
449 beginfit2[0] = dens2;
450 beginfit2[1] = dens1;
451 beginfit2[2] = (3*splitpoint/2);
452 beginfit2[3] = 0.5;
454 snew(zDensavg, zslices);
455 snew(sigma1, zslices);
456 snew(sigma2, zslices);
458 for (k = 0; k < zslices; k++)
460 sigma1[k] = sigma2[k] = 1;
462 /*Calculate average density along z - avoid smoothing by using coarse-grained-mesh*/
463 for (k = 0; k < zslices; k++)
465 for (n = 0; n < tblocks; n++)
467 for (i = 0; i < xslices; i++)
469 for (j = 0; j < yslices; j++)
471 zDensavg[k] += (Densmap[n][i][j][k]/(xslices*yslices*tblocks));
477 if (debug)
479 xvg = xvgropen("DensprofileonZ.xvg", "Averaged Densityprofile on Z", "z[nm]", "Density[kg/m^3]", oenv);
480 for (k = 0; k < zslices; k++)
482 fprintf(xvg, "%4f.3 %8f.4\n", k*binwidth, zDensavg[k]);
484 xvgrclose(xvg);
487 /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/
489 /*Fit 1st half of box*/
490 do_lmfit(zslices, zDensavg, sigma1, binwidth, nullptr, startpoint, splitpoint, oenv, FALSE, effnERF, beginfit1, 8, nullptr);
491 /*Fit 2nd half of box*/
492 do_lmfit(zslices, zDensavg, sigma2, binwidth, nullptr, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 8, nullptr);
494 /*Initialise the const arrays for storing the average fit parameters*/
495 avgfit1 = beginfit1;
496 avgfit2 = beginfit2;
500 /*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*/
501 for (n = 0; n < tblocks; n++)
503 for (i = 0; i < xslices; i++)
505 for (j = 0; j < yslices; j++)
507 /*Reinitialise fit for each mesh-point*/
508 srenew(fit1, 4);
509 srenew(fit2, 4);
510 for (k = 0; k < 4; k++)
512 fit1[k] = avgfit1[k];
513 fit2[k] = avgfit2[k];
515 /*Now fit and store in structures in row-major order int[n][i][j]*/
516 do_lmfit(zslices, Densmap[n][i][j], sigma1, binwidth, nullptr, startpoint, splitpoint, oenv, FALSE, effnERF, fit1, 0, nullptr);
517 int1[n][j+(yslices*i)]->Z = fit1[2];
518 int1[n][j+(yslices*i)]->t = fit1[3];
519 do_lmfit(zslices, Densmap[n][i][j], sigma2, binwidth, nullptr, splitpoint, endpoint, oenv, FALSE, effnERF, fit2, 0, nullptr);
520 int2[n][j+(yslices*i)]->Z = fit2[2];
521 int2[n][j+(yslices*i)]->t = fit2[3];
528 *intf1 = int1;
529 *intf2 = int2;
533 static void writesurftoxpms(t_interf ***surf1, t_interf ***surf2, int tblocks, int xbins, int ybins, int zbins, real bw, real bwz, gmx::ArrayRef<const std::string> outfiles, int maplevels )
535 char numbuf[STRLEN];
536 int n, i, j;
537 real **profile1, **profile2;
538 real max1, max2, min1, min2, *xticks, *yticks;
539 t_rgb lo = {0, 0, 0};
540 t_rgb hi = {1, 1, 1};
541 FILE *xpmfile1, *xpmfile2;
543 /*Prepare xpm structures for output*/
545 /*Allocate memory to tick's and matrices*/
546 snew (xticks, xbins+1);
547 snew (yticks, ybins+1);
549 profile1 = mk_matrix(xbins, ybins, FALSE);
550 profile2 = mk_matrix(xbins, ybins, FALSE);
552 for (i = 0; i < xbins+1; i++)
554 xticks[i] += bw;
556 for (j = 0; j < ybins+1; j++)
558 yticks[j] += bw;
561 xpmfile1 = gmx_ffopen(outfiles[0], "w");
562 xpmfile2 = gmx_ffopen(outfiles[1], "w");
564 max1 = max2 = 0.0;
565 min1 = min2 = zbins*bwz;
567 for (n = 0; n < tblocks; n++)
569 sprintf(numbuf, "tblock: %4i", n);
570 /*Filling matrices for inclusion in xpm-files*/
571 for (i = 0; i < xbins; i++)
573 for (j = 0; j < ybins; j++)
575 profile1[i][j] = (surf1[n][j+ybins*i])->Z;
576 profile2[i][j] = (surf2[n][j+ybins*i])->Z;
577 /*Finding max and min values*/
578 if (profile1[i][j] > max1)
580 max1 = profile1[i][j];
582 if (profile1[i][j] < min1)
584 min1 = profile1[i][j];
586 if (profile2[i][j] > max2)
588 max2 = profile2[i][j];
590 if (profile2[i][j] < min2)
592 min2 = profile2[i][j];
597 write_xpm(xpmfile1, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile1, min1, max1, lo, hi, &maplevels);
598 write_xpm(xpmfile2, 3, numbuf, "Height", "x[nm]", "y[nm]", xbins, ybins, xticks, yticks, profile2, min2, max2, lo, hi, &maplevels);
601 gmx_ffclose(xpmfile1);
602 gmx_ffclose(xpmfile2);
605 sfree(profile1);
606 sfree(profile2);
607 sfree(xticks);
608 sfree(yticks);
611 static void writeraw(t_interf ***int1, t_interf ***int2, int tblocks,
612 int xbins, int ybins, gmx::ArrayRef<const std::string> fnms,
613 const gmx_output_env_t *oenv)
615 FILE *raw1, *raw2;
616 int i, j, n;
618 raw1 = gmx_ffopen(fnms[0], "w");
619 raw2 = gmx_ffopen(fnms[1], "w");
622 gmx::BinaryInformationSettings settings;
623 settings.generatedByHeader(true);
624 settings.linePrefix("# ");
625 gmx::printBinaryInformation(raw1, output_env_get_program_context(oenv),
626 settings);
627 gmx::printBinaryInformation(raw2, output_env_get_program_context(oenv),
628 settings);
630 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
631 fprintf(raw1, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
632 fprintf(raw2, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
633 fprintf(raw1, "%i %i %i\n", tblocks, xbins, ybins);
634 fprintf(raw2, "%i %i %i\n", tblocks, xbins, ybins);
635 for (n = 0; n < tblocks; n++)
637 for (i = 0; i < xbins; i++)
639 for (j = 0; j < ybins; j++)
641 fprintf(raw1, "%i %i %8.5f %6.4f\n", i, j, (int1[n][j+ybins*i])->Z, (int1[n][j+ybins*i])->t);
642 fprintf(raw2, "%i %i %8.5f %6.4f\n", i, j, (int2[n][j+ybins*i])->Z, (int2[n][j+ybins*i])->t);
647 gmx_ffclose(raw1);
648 gmx_ffclose(raw2);
653 int gmx_densorder(int argc, char *argv[])
655 static const char *desc[] = {
656 "[THISMODULE] reduces a two-phase density distribution",
657 "along an axis, computed over a MD trajectory,",
658 "to 2D surfaces fluctuating in time, by a fit to",
659 "a functional profile for interfacial densities.",
660 "A time-averaged spatial representation of the",
661 "interfaces can be output with the option [TT]-tavg[tt]."
664 /* Extra arguments - but note how you always get the begin/end
665 * options when running the program, without mentioning them here!
668 gmx_output_env_t *oenv;
669 t_topology *top;
670 char **grpname;
671 int ePBC, *ngx;
672 static real binw = 0.2;
673 static real binwz = 0.05;
674 static real dens1 = 0.00;
675 static real dens2 = 1000.00;
676 static int ftorder = 0;
677 static int nsttblock = 100;
678 static int axis = 2;
679 static const char *axtitle = "Z";
680 int **index; /* Index list for single group*/
681 int xslices, yslices, zslices, tblock;
682 static gmx_bool bGraph = FALSE;
683 static gmx_bool bCenter = FALSE;
684 static gmx_bool bFourier = FALSE;
685 static gmx_bool bRawOut = FALSE;
686 static gmx_bool bOut = FALSE;
687 static gmx_bool b1d = FALSE;
688 static int nlevels = 100;
689 /*Densitymap - Densmap[t][x][y][z]*/
690 real ****Densmap = nullptr;
691 /* Surfaces surf[t][surf_x,surf_y]*/
692 t_interf ***surf1, ***surf2;
694 static const char *meth[] = {nullptr, "bisect", "functional", nullptr};
695 int eMeth;
697 t_pargs pa[] = {
698 { "-1d", FALSE, etBOOL, {&b1d},
699 "Pseudo-1d interface geometry"},
700 { "-bw", FALSE, etREAL, {&binw},
701 "Binwidth of density distribution tangential to interface"},
702 { "-bwn", FALSE, etREAL, {&binwz},
703 "Binwidth of density distribution normal to interface"},
704 { "-order", FALSE, etINT, {&ftorder},
705 "Order of Gaussian filter, order 0 equates to NO filtering"},
706 {"-axis", FALSE, etSTR, {&axtitle},
707 "Axis Direction - X, Y or Z"},
708 {"-method", FALSE, etENUM, {meth},
709 "Interface location method"},
710 {"-d1", FALSE, etREAL, {&dens1},
711 "Bulk density phase 1 (at small z)"},
712 {"-d2", FALSE, etREAL, {&dens2},
713 "Bulk density phase 2 (at large z)"},
714 { "-tblock", FALSE, etINT, {&nsttblock},
715 "Number of frames in one time-block average"},
716 { "-nlevel", FALSE, etINT, {&nlevels},
717 "Number of Height levels in 2D - XPixMaps"}
721 t_filenm fnm[] = {
722 { efTPR, "-s", nullptr, ffREAD }, /* this is for the topology */
723 { efTRX, "-f", nullptr, ffREAD }, /* and this for the trajectory */
724 { efNDX, "-n", nullptr, ffREAD}, /* this is to select groups */
725 { efDAT, "-o", "Density4D", ffOPTWR}, /* This is for outputting the entire 4D densityfield in binary format */
726 { efOUT, "-or", nullptr, ffOPTWRMULT}, /* This is for writing out the entire information in the t_interf arrays */
727 { efXPM, "-og", "interface", ffOPTWRMULT}, /* This is for writing out the interface meshes - one xpm-file per tblock*/
728 { efOUT, "-Spect", "intfspect", ffOPTWRMULT}, /* This is for the trajectory averaged Fourier-spectra*/
731 #define NFILE asize(fnm)
733 /* This is the routine responsible for adding default options,
734 * calling the X/motif interface, etc. */
735 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
736 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
738 return 0;
742 eMeth = nenum(meth);
743 bFourier = opt2bSet("-Spect", NFILE, fnm);
744 bRawOut = opt2bSet("-or", NFILE, fnm);
745 bGraph = opt2bSet("-og", NFILE, fnm);
746 bOut = opt2bSet("-o", NFILE, fnm);
747 top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
748 snew(grpname, 1);
749 snew(index, 1);
750 snew(ngx, 1);
752 /* Calculate axis */
753 axis = toupper(axtitle[0]) - 'X';
755 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, ngx, index, grpname);
757 density_in_time(ftp2fn(efTRX, NFILE, fnm), index, ngx, binw, binwz, nsttblock, &Densmap, &xslices, &yslices, &zslices, &tblock, top, ePBC, axis, bCenter, b1d, oenv);
759 if (ftorder > 0)
761 filterdensmap(Densmap, xslices, yslices, zslices, tblock, 2*ftorder+1);
764 if (bOut)
766 outputfield(opt2fn("-o", NFILE, fnm), Densmap, xslices, yslices, zslices, tblock);
769 interfaces_txy(Densmap, xslices, yslices, zslices, tblock, binwz, eMeth, dens1, dens2, &surf1, &surf2, oenv);
771 if (bGraph)
774 /*Output surface-xpms*/
775 gmx::ArrayRef<const std::string> graphFiles = opt2fns("-og", NFILE, fnm);
776 if (graphFiles.size() != 2)
778 gmx_fatal(FARGS, "No or not correct number (2) of output-files: %td", graphFiles.ssize());
780 writesurftoxpms(surf1, surf2, tblock, xslices, yslices, zslices, binw, binwz, graphFiles, zslices);
787 /*Output raw-data*/
788 if (bRawOut)
790 gmx::ArrayRef<const std::string> rawFiles = opt2fns("-or", NFILE, fnm);
791 if (rawFiles.size() != 2)
793 gmx_fatal(FARGS, "No or not correct number (2) of output-files: %td", rawFiles.ssize());
795 writeraw(surf1, surf2, tblock, xslices, yslices, rawFiles, oenv);
800 if (bFourier)
802 gmx::ArrayRef<const std::string> spectra = opt2fns("-Spect", NFILE, fnm);
803 if (spectra.size() != 2)
805 gmx_fatal(FARGS, "No or not correct number (2) of output-file-series: %td",
806 spectra.ssize());
808 powerspectavg_intf(surf1, surf2, tblock, xslices, yslices, spectra);
811 sfree(Densmap);
812 if (bGraph || bFourier || bRawOut)
814 sfree(surf1);
815 sfree(surf2);
818 return 0;