1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 * $Id: densorder.c,v 0.9
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-2001
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
34 * Gyas ROwers Mature At Cryogenic Speed
51 #include "gromacs/fileio/futil.h"
52 #include "gromacs/commandline/pargs.h"
54 #include "gromacs/fileio/tpxio.h"
55 #include "gromacs/fileio/trxio.h"
57 #include "gromacs/fileio/matio.h"
58 #include "dens_filter.h"
59 #include "binsearch.h"
60 #include "powerspect.h"
64 #include "gromacs/utility/exceptions.h"
65 #include "gromacs/utility/programinfo.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 = 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 real *fit1 = NULL, *fit2 = NULL;
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, 3);
496 /*Fit 2nd half of box*/
497 do_lmfit(zslices, zDensavg, sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 3);
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, 1);
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, 2);
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 = ffopen(outfiles[0], "w");
567 xpmfile2 = 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);
616 static void writeraw(t_interf ***int1, t_interf ***int2, int tblocks, int xbins, int ybins, char **fnms)
621 raw1 = ffopen(fnms[0], "w");
622 raw2 = ffopen(fnms[1], "w");
625 gmx::BinaryInformationSettings settings;
626 settings.generatedByHeader(true);
627 settings.linePrefix("# ");
628 gmx::printBinaryInformation(raw1, gmx::ProgramInfo::getInstance(),
630 gmx::printBinaryInformation(raw2, gmx::ProgramInfo::getInstance(),
633 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
634 fprintf(raw1, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
635 fprintf(raw2, "# Legend: nt nx ny\n# Xbin Ybin Z t\n");
636 fprintf(raw1, "%i %i %i\n", tblocks, xbins, ybins);
637 fprintf(raw2, "%i %i %i\n", tblocks, xbins, ybins);
638 for (n = 0; n < tblocks; n++)
640 for (i = 0; i < xbins; i++)
642 for (j = 0; j < ybins; j++)
644 fprintf(raw1, "%i %i %8.5f %6.4f\n", i, j, (int1[n][j+ybins*i])->Z, (int1[n][j+ybins*i])->t);
645 fprintf(raw2, "%i %i %8.5f %6.4f\n", i, j, (int2[n][j+ybins*i])->Z, (int2[n][j+ybins*i])->t);
656 int gmx_densorder(int argc, char *argv[])
658 static const char *desc[] = {
659 "[THISMODULE] reduces a two-phase density distribution",
660 "along an axis, computed over a MD trajectory,",
661 "to 2D surfaces fluctuating in time, by a fit to",
662 "a functional profile for interfacial densities.",
663 "A time-averaged spatial representation of the",
664 "interfaces can be output with the option [TT]-tavg[tt]."
667 /* Extra arguments - but note how you always get the begin/end
668 * options when running the program, without mentioning them here!
675 static real binw = 0.2;
676 static real binwz = 0.05;
677 static real dens1 = 0.00;
678 static real dens2 = 1000.00;
679 static int ftorder = 0;
680 static int nsttblock = 100;
682 static const char *axtitle = "Z";
683 atom_id **index; /* Index list for single group*/
684 int xslices, yslices, zslices, tblock;
685 static gmx_bool bGraph = FALSE;
686 static gmx_bool bCenter = FALSE;
687 static gmx_bool bFourier = FALSE;
688 static gmx_bool bRawOut = FALSE;
689 static gmx_bool bOut = FALSE;
690 static gmx_bool b1d = FALSE;
691 static int nlevels = 100;
692 /*Densitymap - Densmap[t][x][y][z]*/
693 real ****Densmap = NULL;
694 /* Surfaces surf[t][surf_x,surf_y]*/
695 t_interf ***surf1, ***surf2;
697 static const char *meth[] = {NULL, "bisect", "functional", NULL};
700 char **graphfiles, **rawfiles, **spectra; /* Filenames for xpm-surface maps, rawdata and powerspectra */
701 int nfxpm = -1, nfraw, nfspect; /* # files for interface maps and spectra = # interfaces */
704 { "-1d", FALSE, etBOOL, {&b1d},
705 "Pseudo-1d interface geometry"},
706 { "-bw", FALSE, etREAL, {&binw},
707 "Binwidth of density distribution tangential to interface"},
708 { "-bwn", FALSE, etREAL, {&binwz},
709 "Binwidth of density distribution normal to interface"},
710 { "-order", FALSE, etINT, {&ftorder},
711 "Order of Gaussian filter, order 0 equates to NO filtering"},
712 {"-axis", FALSE, etSTR, {&axtitle},
713 "Axis Direction - X, Y or Z"},
714 {"-method", FALSE, etENUM, {meth},
715 "Interface location method"},
716 {"-d1", FALSE, etREAL, {&dens1},
717 "Bulk density phase 1 (at small z)"},
718 {"-d2", FALSE, etREAL, {&dens2},
719 "Bulk density phase 2 (at large z)"},
720 { "-tblock", FALSE, etINT, {&nsttblock},
721 "Number of frames in one time-block average"},
722 { "-nlevel", FALSE, etINT, {&nlevels},
723 "Number of Height levels in 2D - XPixMaps"}
728 { efTPX, "-s", NULL, ffREAD }, /* this is for the topology */
729 { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */
730 { efNDX, "-n", NULL, ffREAD}, /* this is to select groups */
731 { efDAT, "-o", "Density4D", ffOPTWR}, /* This is for outputting the entire 4D densityfield in binary format */
732 { efOUT, "-or", NULL, ffOPTWRMULT}, /* This is for writing out the entire information in the t_interf arrays */
733 { efXPM, "-og", "interface", ffOPTWRMULT}, /* This is for writing out the interface meshes - one xpm-file per tblock*/
734 { efOUT, "-Spect", "intfspect", ffOPTWRMULT}, /* This is for the trajectory averaged Fourier-spectra*/
737 #define NFILE asize(fnm)
739 /* This is the routine responsible for adding default options,
740 * calling the X/motif interface, etc. */
741 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
742 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
749 bFourier = opt2bSet("-Spect", NFILE, fnm);
750 bRawOut = opt2bSet("-or", NFILE, fnm);
751 bGraph = opt2bSet("-og", NFILE, fnm);
752 bOut = opt2bSet("-o", NFILE, fnm);
753 top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
759 axis = toupper(axtitle[0]) - 'X';
761 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, ngx, index, grpname);
763 density_in_time(ftp2fn(efTRX, NFILE, fnm), index, ngx, binw, binwz, nsttblock, &Densmap, &xslices, &yslices, &zslices, &tblock, top, ePBC, axis, bCenter, b1d, oenv);
767 filterdensmap(Densmap, xslices, yslices, zslices, tblock, 2*ftorder+1);
772 outputfield(opt2fn("-o", NFILE, fnm), Densmap, xslices, yslices, zslices, tblock);
775 interfaces_txy(Densmap, xslices, yslices, zslices, tblock, binwz, eMeth, dens1, dens2, &surf1, &surf2, oenv);
780 /*Output surface-xpms*/
781 nfxpm = opt2fns(&graphfiles, "-og", NFILE, fnm);
784 gmx_fatal(FARGS, "No or not correct number (2) of output-files: %d", nfxpm);
786 writesurftoxpms(surf1, surf2, tblock, xslices, yslices, zslices, binw, binwz, graphfiles, zslices);
796 nfraw = opt2fns(&rawfiles, "-or", NFILE, fnm);
799 gmx_fatal(FARGS, "No or not correct number (2) of output-files: %d", nfxpm);
801 writeraw(surf1, surf2, tblock, xslices, yslices, rawfiles);
808 nfspect = opt2fns(&spectra, "-Spect", NFILE, fnm);
811 gmx_fatal(FARGS, "No or not correct number (2) of output-file-series: %d",
814 powerspectavg_intf(surf1, surf2, tblock, xslices, yslices, spectra);
818 if (bGraph || bFourier || bRawOut)