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