3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
57 int gmx_filter(int argc, char *argv[])
59 const char *desc[] = {
60 "[TT]g_filter[tt] performs frequency filtering on a trajectory.",
61 "The filter shape is cos([GRK]pi[grk] t/A) + 1 from -A to +A, where A is given",
62 "by the option [TT]-nf[tt] times the time step in the input trajectory.",
63 "This filter reduces fluctuations with period A by 85%, with period",
64 "2*A by 50% and with period 3*A by 17% for low-pass filtering.",
65 "Both a low-pass and high-pass filtered trajectory can be written.[PAR]",
67 "Option [TT]-ol[tt] writes a low-pass filtered trajectory.",
68 "A frame is written every [TT]-nf[tt] input frames.",
69 "This ratio of filter length and output interval ensures a good",
70 "suppression of aliasing of high-frequency motion, which is useful for",
71 "making smooth movies. Also averages of properties which are linear",
72 "in the coordinates are preserved, since all input frames are weighted",
73 "equally in the output.",
74 "When all frames are needed, use the [TT]-all[tt] option.[PAR]",
76 "Option [TT]-oh[tt] writes a high-pass filtered trajectory.",
77 "The high-pass filtered coordinates are added to the coordinates",
78 "from the structure file. When using high-pass filtering use [TT]-fit[tt]",
79 "or make sure you use a trajectory that has been fitted on",
80 "the coordinates in the structure file."
84 static gmx_bool bNoJump = TRUE, bFit = FALSE, bLowAll = FALSE;
86 { "-nf", FALSE, etINT, {&nf},
87 "Sets the filter length as well as the output interval for low-pass filtering" },
88 { "-all", FALSE, etBOOL, {&bLowAll},
89 "Write all low-pass filtered frames" },
90 { "-nojump", FALSE, etBOOL, {&bNoJump},
91 "Remove jumps of atoms across the box" },
92 { "-fit", FALSE, etBOOL, {&bFit},
93 "Fit all frames to a reference structure" }
95 const char *topfile, *lowfile, *highfile;
96 gmx_bool bTop = FALSE;
100 matrix topbox, *box, boxf;
101 char title[256], *grpname;
106 t_trxstatus *outl, *outh;
107 int nffr, i, fr, nat, j, d, m;
109 real flen, *filt, sum, *t;
110 rvec xcmtop, xcm, **x, *ptr, *xf, *xn, *xp, hbox;
112 gmx_rmpbc_t gpbc = NULL;
114 #define NLEG asize(leg)
116 { efTRX, "-f", NULL, ffREAD },
117 { efTPS, NULL, NULL, ffOPTRD },
118 { efNDX, NULL, NULL, ffOPTRD },
119 { efTRO, "-ol", "lowpass", ffOPTWR },
120 { efTRO, "-oh", "highpass", ffOPTWR }
122 #define NFILE asize(fnm)
124 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
125 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
127 highfile = opt2fn_null("-oh", NFILE, fnm);
130 topfile = ftp2fn(efTPS, NFILE, fnm);
131 lowfile = opt2fn_null("-ol", NFILE, fnm);
135 topfile = ftp2fn_null(efTPS, NFILE, fnm);
136 lowfile = opt2fn("-ol", NFILE, fnm);
140 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC,
141 &xtop, NULL, topbox, TRUE);
144 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr, topbox);
145 gmx_rmpbc(gpbc, top.atoms.nr, topbox, xtop);
152 fprintf(stderr, "Select group for least squares fit\n");
153 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
155 snew(w_rls, top.atoms.nr);
156 for (i = 0; i < isize; i++)
158 w_rls[index[i]] = top.atoms.atom[index[i]].m;
160 calc_xcm(xtop, isize, index, top.atoms.atom, xcmtop, FALSE);
161 for (j = 0; j < top.atoms.nr; j++)
163 rvec_dec(xtop[j], xcmtop);
167 /* The actual filter length flen can actually be any real number */
169 /* nffr is the number of frames that we filter over */
173 for (i = 0; i < nffr; i++)
175 filt[i] = cos(2*M_PI*(i - nf + 1)/(real)flen) + 1;
178 fprintf(stdout, "filter weights:");
179 for (i = 0; i < nffr; i++)
182 fprintf(stdout, " %5.3f", filt[i]);
184 fprintf(stdout, "\n");
190 nat = read_first_x(oenv, &in, opt2fn("-f", NFILE, fnm),
191 &(t[nffr - 1]), &(x[nffr - 1]), box[nffr - 1]);
193 for (i = 0; i < nat; i++)
197 /* x[nffr - 1] was already allocated by read_first_x */
198 for (i = 0; i < nffr-1; i++)
205 outl = open_trx(lowfile, "w");
213 outh = open_trx(highfile, "w");
224 if (bNoJump && fr > 0)
227 for (j = 0; j < nat; j++)
229 for (d = 0; d < DIM; d++)
231 hbox[d] = 0.5*box[nffr - 1][d][d];
234 for (i = 0; i < nat; i++)
236 for (m = DIM-1; m >= 0; m--)
240 while (xn[i][m] - xp[i][m] <= -hbox[m])
242 for (d = 0; d <= m; d++)
244 xn[i][d] += box[nffr - 1][m][d];
247 while (xn[i][m] - xp[i][m] > hbox[m])
249 for (d = 0; d <= m; d++)
251 xn[i][d] -= box[nffr - 1][m][d];
260 gmx_rmpbc(gpbc, nat, box[nffr - 1], xn);
264 calc_xcm(xn, isize, index, top.atoms.atom, xcm, FALSE);
265 for (j = 0; j < nat; j++)
267 rvec_dec(xn[j], xcm);
269 do_fit(nat, w_rls, xtop, xn);
270 for (j = 0; j < nat; j++)
272 rvec_inc(xn[j], xcmtop);
275 if (fr >= nffr && (outh || bLowAll || fr % nf == nf - 1))
277 /* Lowpass filtering */
278 for (j = 0; j < nat; j++)
283 for (i = 0; i < nffr; i++)
285 for (j = 0; j < nat; j++)
287 for (d = 0; d < DIM; d++)
289 xf[j][d] += filt[i]*x[i][j][d];
292 for (j = 0; j < DIM; j++)
294 for (d = 0; d < DIM; d++)
296 boxf[j][d] += filt[i]*box[i][j][d];
300 if (outl && (bLowAll || fr % nf == nf - 1))
302 write_trx(outl, nat, ind, topfile ? &(top.atoms) : NULL,
303 0, t[nf - 1], bFit ? topbox : boxf, xf, NULL, NULL);
307 /* Highpass filtering */
308 for (j = 0; j < nat; j++)
310 for (d = 0; d < DIM; d++)
312 xf[j][d] = xtop[j][d] + x[nf - 1][j][d] - xf[j][d];
317 for (j = 0; j < nat; j++)
319 rvec_inc(xf[j], xcmtop);
322 for (j = 0; j < DIM; j++)
324 for (d = 0; d < DIM; d++)
326 boxf[j][d] = topbox[j][d] + box[nf - 1][j][d] - boxf[j][d];
329 write_trx(outh, nat, ind, topfile ? &(top.atoms) : NULL,
330 0, t[nf - 1], bFit ? topbox : boxf, xf, NULL, NULL);
333 /* Cycle all the pointer and the box by one */
335 for (i = 0; i < nffr-1; i++)
339 copy_mat(box[i+1], box[i]);
344 while (read_next_x(oenv, in, &(t[nffr - 1]), nat, x[nffr - 1], box[nffr - 1]));
348 gmx_rmpbc_done(gpbc);