src/gromacs/gmxana/gmx_filter.c

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