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