src/gromacs/gmxana/gmx_filter.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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, 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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstring>
  41
  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/utilities.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/rmpbc.h"
  51 #include "gromacs/topology/index.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/arraysize.h"
  54 #include "gromacs/utility/smalloc.h"
  55
  56 int gmx_filter(int argc, char* argv[])
  57 {
  58     const char* desc[] = {
  59         "[THISMODULE] performs frequency filtering on a trajectory.",
  60         "The filter shape is cos([GRK]pi[grk] t/A) + 1 from -A to +A, where A is given",
  61         "by the option [TT]-nf[tt] times the time step in the input trajectory.",
  62         "This filter reduces fluctuations with period A by 85%, with period",
  63         "2*A by 50% and with period 3*A by 17% for low-pass filtering.",
  64         "Both a low-pass and high-pass filtered trajectory can be written.[PAR]",
  65
  66         "Option [TT]-ol[tt] writes a low-pass filtered trajectory.",
  67         "A frame is written every [TT]-nf[tt] input frames.",
  68         "This ratio of filter length and output interval ensures a good",
  69         "suppression of aliasing of high-frequency motion, which is useful for",
  70         "making smooth movies. Also averages of properties which are linear",
  71         "in the coordinates are preserved, since all input frames are weighted",
  72         "equally in the output.",
  73         "When all frames are needed, use the [TT]-all[tt] option.[PAR]",
  74
  75         "Option [TT]-oh[tt] writes a high-pass filtered trajectory.",
  76         "The high-pass filtered coordinates are added to the coordinates",
  77         "from the structure file. When using high-pass filtering use [TT]-fit[tt]",
  78         "or make sure you use a trajectory that has been fitted on",
  79         "the coordinates in the structure file."
  80     };
  81
  82     static int      nf      = 10;
  83     static gmx_bool bNoJump = TRUE, bFit = FALSE, bLowAll = FALSE;
  84     t_pargs         pa[] = {
  85         { "-nf",
  86           FALSE,
  87           etINT,
  88           { &nf },
  89           "Sets the filter length as well as the output interval for low-pass filtering" },
  90         { "-all", FALSE, etBOOL, { &bLowAll }, "Write all low-pass filtered frames" },
  91         { "-nojump", FALSE, etBOOL, { &bNoJump }, "Remove jumps of atoms across the box" },
  92         { "-fit", FALSE, etBOOL, { &bFit }, "Fit all frames to a reference structure" }
  93     };
  94     const char *      topfile, *lowfile, *highfile;
  95     gmx_bool          bTop = FALSE;
  96     t_topology        top;
  97     int               ePBC = -1;
  98     rvec*             xtop;
  99     matrix            topbox, *box, boxf;
 100     char*             grpname;
 101     int               isize;
 102     int*              index;
 103     real*             w_rls = nullptr;
 104     t_trxstatus*      in;
 105     t_trxstatus *     outl, *outh;
 106     int               nffr, i, fr, nat, j, d, m;
 107     int*              ind;
 108     real              flen, *filt, sum, *t;
 109     rvec              xcmtop, xcm, **x, *ptr, *xf, *xn, *xp, hbox;
 110     gmx_output_env_t* oenv;
 111     gmx_rmpbc_t       gpbc = nullptr;
 112
 113 #define NLEG asize(leg)
 114     t_filenm fnm[] = { { efTRX, "-f", nullptr, ffREAD },
 115                        { efTPS, nullptr, nullptr, ffOPTRD },
 116                        { efNDX, nullptr, nullptr, ffOPTRD },
 117                        { efTRO, "-ol", "lowpass", ffOPTWR },
 118                        { efTRO, "-oh", "highpass", ffOPTWR } };
 119 #define NFILE asize(fnm)
 120
 121     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
 122                            asize(desc), desc, 0, nullptr, &oenv))
 123     {
 124         return 0;
 125     }
 126
 127     highfile = opt2fn_null("-oh", NFILE, fnm);
 128     if (highfile)
 129     {
 130         topfile = ftp2fn(efTPS, NFILE, fnm);
 131         lowfile = opt2fn_null("-ol", NFILE, fnm);
 132     }
 133     else
 134     {
 135         topfile = ftp2fn_null(efTPS, NFILE, fnm);
 136         lowfile = opt2fn("-ol", NFILE, fnm);
 137     }
 138     if (topfile)
 139     {
 140         bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, nullptr, topbox, TRUE);
 141         if (bTop)
 142         {
 143             gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
 144             gmx_rmpbc(gpbc, top.atoms.nr, topbox, xtop);
 145         }
 146     }
 147
 148     clear_rvec(xcmtop);
 149     if (bFit)
 150     {
 151         fprintf(stderr, "Select group for least squares fit\n");
 152         get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 153         /* Set the weight */
 154         snew(w_rls, top.atoms.nr);
 155         for (i = 0; i < isize; i++)
 156         {
 157             w_rls[index[i]] = top.atoms.atom[index[i]].m;
 158         }
 159         calc_xcm(xtop, isize, index, top.atoms.atom, xcmtop, FALSE);
 160         for (j = 0; j < top.atoms.nr; j++)
 161         {
 162             rvec_dec(xtop[j], xcmtop);
 163         }
 164     }
 165
 166     /* The actual filter length flen can actually be any real number */
 167     flen = 2 * nf;
 168     /* nffr is the number of frames that we filter over */
 169     nffr = 2 * nf - 1;
 170     snew(filt, nffr);
 171     sum = 0;
 172     for (i = 0; i < nffr; i++)
 173     {
 174         filt[i] = std::cos(2 * M_PI * (i - nf + 1) / static_cast<real>(flen)) + 1;
 175         sum += filt[i];
 176     }
 177     fprintf(stdout, "filter weights:");
 178     for (i = 0; i < nffr; i++)
 179     {
 180         filt[i] /= sum;
 181         fprintf(stdout, " %5.3f", filt[i]);
 182     }
 183     fprintf(stdout, "\n");
 184
 185     snew(t, nffr);
 186     snew(x, nffr);
 187     snew(box, nffr);
 188
 189     nat = read_first_x(oenv, &in, opt2fn("-f", NFILE, fnm), &(t[nffr - 1]), &(x[nffr - 1]), box[nffr - 1]);
 190     snew(ind, nat);
 191     for (i = 0; i < nat; i++)
 192     {
 193         ind[i] = i;
 194     }
 195     /* x[nffr - 1] was already allocated by read_first_x */
 196     for (i = 0; i < nffr - 1; i++)
 197     {
 198         snew(x[i], nat);
 199     }
 200     snew(xf, nat);
 201     if (lowfile)
 202     {
 203         outl = open_trx(lowfile, "w");
 204     }
 205     else
 206     {
 207         outl = nullptr;
 208     }
 209     if (highfile)
 210     {
 211         outh = open_trx(highfile, "w");
 212     }
 213     else
 214     {
 215         outh = nullptr;
 216     }
 217
 218     fr = 0;
 219     do
 220     {
 221         xn = x[nffr - 1];
 222         if (bNoJump && fr > 0)
 223         {
 224             xp = x[nffr - 2];
 225             for (j = 0; j < nat; j++)
 226             {
 227                 for (d = 0; d < DIM; d++)
 228                 {
 229                     hbox[d] = 0.5 * box[nffr - 1][d][d];
 230                 }
 231             }
 232             for (i = 0; i < nat; i++)
 233             {
 234                 for (m = DIM - 1; m >= 0; m--)
 235                 {
 236                     if (hbox[m] > 0)
 237                     {
 238                         while (xn[i][m] - xp[i][m] <= -hbox[m])
 239                         {
 240                             for (d = 0; d <= m; d++)
 241                             {
 242                                 xn[i][d] += box[nffr - 1][m][d];
 243                             }
 244                         }
 245                         while (xn[i][m] - xp[i][m] > hbox[m])
 246                         {
 247                             for (d = 0; d <= m; d++)
 248                             {
 249                                 xn[i][d] -= box[nffr - 1][m][d];
 250                             }
 251                         }
 252                     }
 253                 }
 254             }
 255         }
 256         if (bTop)
 257         {
 258             gmx_rmpbc(gpbc, nat, box[nffr - 1], xn);
 259         }
 260         if (bFit)
 261         {
 262             calc_xcm(xn, isize, index, top.atoms.atom, xcm, FALSE);
 263             for (j = 0; j < nat; j++)
 264             {
 265                 rvec_dec(xn[j], xcm);
 266             }
 267             do_fit(nat, w_rls, xtop, xn);
 268             for (j = 0; j < nat; j++)
 269             {
 270                 rvec_inc(xn[j], xcmtop);
 271             }
 272         }
 273         if (fr >= nffr && (outh || bLowAll || fr % nf == nf - 1))
 274         {
 275             /* Lowpass filtering */
 276             for (j = 0; j < nat; j++)
 277             {
 278                 clear_rvec(xf[j]);
 279             }
 280             clear_mat(boxf);
 281             for (i = 0; i < nffr; i++)
 282             {
 283                 for (j = 0; j < nat; j++)
 284                 {
 285                     for (d = 0; d < DIM; d++)
 286                     {
 287                         xf[j][d] += filt[i] * x[i][j][d];
 288                     }
 289                 }
 290                 for (j = 0; j < DIM; j++)
 291                 {
 292                     for (d = 0; d < DIM; d++)
 293                     {
 294                         boxf[j][d] += filt[i] * box[i][j][d];
 295                     }
 296                 }
 297             }
 298             if (outl && (bLowAll || fr % nf == nf - 1))
 299             {
 300                 write_trx(outl, nat, ind, topfile ? &(top.atoms) : nullptr, 0, t[nf - 1],
 301                           bFit ? topbox : boxf, xf, nullptr, nullptr);
 302             }
 303             if (outh)
 304             {
 305                 /* Highpass filtering */
 306                 for (j = 0; j < nat; j++)
 307                 {
 308                     for (d = 0; d < DIM; d++)
 309                     {
 310                         xf[j][d] = xtop[j][d] + x[nf - 1][j][d] - xf[j][d];
 311                     }
 312                 }
 313                 if (bFit)
 314                 {
 315                     for (j = 0; j < nat; j++)
 316                     {
 317                         rvec_inc(xf[j], xcmtop);
 318                     }
 319                 }
 320                 for (j = 0; j < DIM; j++)
 321                 {
 322                     for (d = 0; d < DIM; d++)
 323                     {
 324                         boxf[j][d] = topbox[j][d] + box[nf - 1][j][d] - boxf[j][d];
 325                     }
 326                 }
 327                 write_trx(outh, nat, ind, topfile ? &(top.atoms) : nullptr, 0, t[nf - 1],
 328                           bFit ? topbox : boxf, xf, nullptr, nullptr);
 329             }
 330         }
 331         /* Cycle all the pointer and the box by one */
 332         ptr = x[0];
 333         for (i = 0; i < nffr - 1; i++)
 334         {
 335             t[i] = t[i + 1];
 336             x[i] = x[i + 1];
 337             copy_mat(box[i + 1], box[i]);
 338         }
 339         x[nffr - 1] = ptr;
 340         fr++;
 341     } while (read_next_x(oenv, in, &(t[nffr - 1]), x[nffr - 1], box[nffr - 1]));
 342
 343     if (bTop)
 344     {
 345         gmx_rmpbc_done(gpbc);
 346     }
 347
 348     if (outh)
 349     {
 350         close_trx(outh);
 351     }
 352     if (outl)
 353     {
 354         close_trx(outl);
 355     }
 356     close_trx(in);
 357
 358     return 0;
 359 }