2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017,2019,2020,2021, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/fileio/confio.h"
44 #include "gromacs/fileio/trxio.h"
45 #include "gromacs/gmxana/gmx_ana.h"
46 #include "gromacs/gmxana/princ.h"
47 #include "gromacs/math/do_fit.h"
48 #include "gromacs/math/units.h"
49 #include "gromacs/math/utilities.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/pbcutil/pbc.h"
52 #include "gromacs/pbcutil/rmpbc.h"
53 #include "gromacs/topology/index.h"
54 #include "gromacs/topology/topology.h"
55 #include "gromacs/utility/arraysize.h"
56 #include "gromacs/utility/smalloc.h"
58 int gmx_filter(int argc, char* argv[])
60 const char* desc[] = {
61 "[THISMODULE] performs frequency filtering on a trajectory.",
62 "The filter shape is cos([GRK]pi[grk] t/A) + 1 from -A to +A, where A is given",
63 "by the option [TT]-nf[tt] times the time step in the input trajectory.",
64 "This filter reduces fluctuations with period A by 85%, with period",
65 "2*A by 50% and with period 3*A by 17% for low-pass filtering.",
66 "Both a low-pass and high-pass filtered trajectory can be written.[PAR]",
68 "Option [TT]-ol[tt] writes a low-pass filtered trajectory.",
69 "A frame is written every [TT]-nf[tt] input frames.",
70 "This ratio of filter length and output interval ensures a good",
71 "suppression of aliasing of high-frequency motion, which is useful for",
72 "making smooth movies. Also averages of properties which are linear",
73 "in the coordinates are preserved, since all input frames are weighted",
74 "equally in the output.",
75 "When all frames are needed, use the [TT]-all[tt] option.[PAR]",
77 "Option [TT]-oh[tt] writes a high-pass filtered trajectory.",
78 "The high-pass filtered coordinates are added to the coordinates",
79 "from the structure file. When using high-pass filtering use [TT]-fit[tt]",
80 "or make sure you use a trajectory that has been fitted on",
81 "the coordinates in the structure file."
85 static gmx_bool bNoJump = TRUE, bFit = FALSE, bLowAll = FALSE;
91 "Sets the filter length as well as the output interval for low-pass filtering" },
92 { "-all", FALSE, etBOOL, { &bLowAll }, "Write all low-pass filtered frames" },
93 { "-nojump", FALSE, etBOOL, { &bNoJump }, "Remove jumps of atoms across the box" },
94 { "-fit", FALSE, etBOOL, { &bFit }, "Fit all frames to a reference structure" }
96 const char * topfile, *lowfile, *highfile;
97 gmx_bool bTop = FALSE;
99 PbcType pbcType = PbcType::Unset;
101 matrix topbox, *box, boxf;
105 real* w_rls = nullptr;
107 t_trxstatus * outl, *outh;
108 int nffr, i, fr, nat, j, d, m;
110 real flen, *filt, sum, *t;
111 rvec xcmtop, xcm, **x, *ptr, *xf, *xn, *xp, hbox;
112 gmx_output_env_t* oenv;
113 gmx_rmpbc_t gpbc = nullptr;
115 #define NLEG asize(leg)
116 t_filenm fnm[] = { { efTRX, "-f", nullptr, ffREAD },
117 { efTPS, nullptr, nullptr, ffOPTRD },
118 { efNDX, nullptr, nullptr, ffOPTRD },
119 { efTRO, "-ol", "lowpass", ffOPTWR },
120 { efTRO, "-oh", "highpass", ffOPTWR } };
121 #define NFILE asize(fnm)
123 if (!parse_common_args(
124 &argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
129 highfile = opt2fn_null("-oh", NFILE, fnm);
132 topfile = ftp2fn(efTPS, NFILE, fnm);
133 lowfile = opt2fn_null("-ol", NFILE, fnm);
137 topfile = ftp2fn_null(efTPS, NFILE, fnm);
138 lowfile = opt2fn("-ol", NFILE, fnm);
142 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xtop, nullptr, topbox, TRUE);
145 gpbc = gmx_rmpbc_init(&top.idef, pbcType, top.atoms.nr);
146 gmx_rmpbc(gpbc, top.atoms.nr, topbox, xtop);
153 fprintf(stderr, "Select group for least squares fit\n");
154 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
156 snew(w_rls, top.atoms.nr);
157 for (i = 0; i < isize; i++)
159 w_rls[index[i]] = top.atoms.atom[index[i]].m;
161 calc_xcm(xtop, isize, index, top.atoms.atom, xcmtop, FALSE);
162 for (j = 0; j < top.atoms.nr; j++)
164 rvec_dec(xtop[j], xcmtop);
168 /* The actual filter length flen can actually be any real number */
170 /* nffr is the number of frames that we filter over */
174 for (i = 0; i < nffr; i++)
176 filt[i] = std::cos(2 * M_PI * (i - nf + 1) / static_cast<real>(flen)) + 1;
179 fprintf(stdout, "filter weights:");
180 for (i = 0; i < nffr; i++)
183 fprintf(stdout, " %5.3f", filt[i]);
185 fprintf(stdout, "\n");
191 nat = read_first_x(oenv, &in, opt2fn("-f", NFILE, fnm), &(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))
305 topfile ? &(top.atoms) : nullptr,
308 bFit ? topbox : boxf,
315 /* Highpass filtering */
316 for (j = 0; j < nat; j++)
318 for (d = 0; d < DIM; d++)
320 xf[j][d] = xtop[j][d] + x[nf - 1][j][d] - xf[j][d];
325 for (j = 0; j < nat; j++)
327 rvec_inc(xf[j], xcmtop);
330 for (j = 0; j < DIM; j++)
332 for (d = 0; d < DIM; d++)
334 boxf[j][d] = topbox[j][d] + box[nf - 1][j][d] - boxf[j][d];
340 topfile ? &(top.atoms) : nullptr,
343 bFit ? topbox : boxf,
349 /* Cycle all the pointer and the box by one */
351 for (i = 0; i < nffr - 1; i++)
355 copy_mat(box[i + 1], box[i]);
359 } while (read_next_x(oenv, in, &(t[nffr - 1]), x[nffr - 1], box[nffr - 1]));
363 gmx_rmpbc_done(gpbc);