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