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/fileio/tpxio.h"
  44 #include "gromacs/fileio/trxio.h"
  45 #include "gromacs/gmxana/gmx_ana.h"
  46 #include "gromacs/gmxana/princ.h"
  47 #include "gromacs/legacyheaders/macros.h"
  48 #include "gromacs/legacyheaders/typedefs.h"
  49 #include "gromacs/math/do_fit.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/pbcutil/rmpbc.h"
  52 #include "gromacs/topology/index.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 int gmx_filter(int argc, char *argv[])
  56 {
  57     const char     *desc[] = {
  58         "[THISMODULE] performs frequency filtering on a trajectory.",
  59         "The filter shape is cos([GRK]pi[grk] t/A) + 1 from -A to +A, where A is given",
  60         "by the option [TT]-nf[tt] times the time step in the input trajectory.",
  61         "This filter reduces fluctuations with period A by 85%, with period",
  62         "2*A by 50% and with period 3*A by 17% for low-pass filtering.",
  63         "Both a low-pass and high-pass filtered trajectory can be written.[PAR]",
  64
  65         "Option [TT]-ol[tt] writes a low-pass filtered trajectory.",
  66         "A frame is written every [TT]-nf[tt] input frames.",
  67         "This ratio of filter length and output interval ensures a good",
  68         "suppression of aliasing of high-frequency motion, which is useful for",
  69         "making smooth movies. Also averages of properties which are linear",
  70         "in the coordinates are preserved, since all input frames are weighted",
  71         "equally in the output.",
  72         "When all frames are needed, use the [TT]-all[tt] option.[PAR]",
  73
  74         "Option [TT]-oh[tt] writes a high-pass filtered trajectory.",
  75         "The high-pass filtered coordinates are added to the coordinates",
  76         "from the structure file. When using high-pass filtering use [TT]-fit[tt]",
  77         "or make sure you use a trajectory that has been fitted on",
  78         "the coordinates in the structure file."
  79     };
  80
  81     static int      nf      = 10;
  82     static gmx_bool bNoJump = TRUE, bFit = FALSE, bLowAll = FALSE;
  83     t_pargs         pa[]    = {
  84         { "-nf", FALSE, etINT, {&nf},
  85           "Sets the filter length as well as the output interval for low-pass filtering" },
  86         { "-all", FALSE, etBOOL, {&bLowAll},
  87           "Write all low-pass filtered frames" },
  88         { "-nojump", FALSE, etBOOL, {&bNoJump},
  89           "Remove jumps of atoms across the box" },
  90         { "-fit", FALSE, etBOOL, {&bFit},
  91           "Fit all frames to a reference structure" }
  92     };
  93     const char     *topfile, *lowfile, *highfile;
  94     gmx_bool        bTop = FALSE;
  95     t_topology      top;
  96     int             ePBC = -1;
  97     rvec           *xtop;
  98     matrix          topbox, *box, boxf;
  99     char            title[256], *grpname;
 100     int             isize;
 101     atom_id        *index;
 102     real           *w_rls = NULL;
 103     t_trxstatus    *in;
 104     t_trxstatus    *outl, *outh;
 105     int             nffr, i, fr, nat, j, d, m;
 106     atom_id        *ind;
 107     real            flen, *filt, sum, *t;
 108     rvec            xcmtop, xcm, **x, *ptr, *xf, *xn, *xp, hbox;
 109     output_env_t    oenv;
 110     gmx_rmpbc_t     gpbc = NULL;
 111
 112 #define NLEG asize(leg)
 113     t_filenm fnm[] = {
 114         { efTRX, "-f", NULL, ffREAD  },
 115         { efTPS, NULL, NULL, ffOPTRD },
 116         { efNDX, NULL, NULL, ffOPTRD },
 117         { efTRO, "-ol", "lowpass",  ffOPTWR },
 118         { efTRO, "-oh", "highpass", ffOPTWR }
 119     };
 120 #define NFILE asize(fnm)
 121
 122     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 123                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &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), title, &top, &ePBC,
 142                              &xtop, NULL, topbox, TRUE);
 143         if (bTop)
 144         {
 145             gpbc = gmx_rmpbc_init(&top.idef, ePBC, 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] = cos(2*M_PI*(i - nf + 1)/(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),
 192                        &(t[nffr - 1]), &(x[nffr - 1]), box[nffr - 1]);
 193     snew(ind, nat);
 194     for (i = 0; i < nat; i++)
 195     {
 196         ind[i] = i;
 197     }
 198     /* x[nffr - 1] was already allocated by read_first_x */
 199     for (i = 0; i < nffr-1; i++)
 200     {
 201         snew(x[i], nat);
 202     }
 203     snew(xf, nat);
 204     if (lowfile)
 205     {
 206         outl = open_trx(lowfile, "w");
 207     }
 208     else
 209     {
 210         outl = 0;
 211     }
 212     if (highfile)
 213     {
 214         outh = open_trx(highfile, "w");
 215     }
 216     else
 217     {
 218         outh = 0;
 219     }
 220
 221     fr = 0;
 222     do
 223     {
 224         xn = x[nffr - 1];
 225         if (bNoJump && fr > 0)
 226         {
 227             xp = x[nffr - 2];
 228             for (j = 0; j < nat; j++)
 229             {
 230                 for (d = 0; d < DIM; d++)
 231                 {
 232                     hbox[d] = 0.5*box[nffr - 1][d][d];
 233                 }
 234             }
 235             for (i = 0; i < nat; i++)
 236             {
 237                 for (m = DIM-1; m >= 0; m--)
 238                 {
 239                     if (hbox[m] > 0)
 240                     {
 241                         while (xn[i][m] - xp[i][m] <= -hbox[m])
 242                         {
 243                             for (d = 0; d <= m; d++)
 244                             {
 245                                 xn[i][d] += box[nffr - 1][m][d];
 246                             }
 247                         }
 248                         while (xn[i][m] - xp[i][m] > hbox[m])
 249                         {
 250                             for (d = 0; d <= m; d++)
 251                             {
 252                                 xn[i][d] -= box[nffr - 1][m][d];
 253                             }
 254                         }
 255                     }
 256                 }
 257             }
 258         }
 259         if (bTop)
 260         {
 261             gmx_rmpbc(gpbc, nat, box[nffr - 1], xn);
 262         }
 263         if (bFit)
 264         {
 265             calc_xcm(xn, isize, index, top.atoms.atom, xcm, FALSE);
 266             for (j = 0; j < nat; j++)
 267             {
 268                 rvec_dec(xn[j], xcm);
 269             }
 270             do_fit(nat, w_rls, xtop, xn);
 271             for (j = 0; j < nat; j++)
 272             {
 273                 rvec_inc(xn[j], xcmtop);
 274             }
 275         }
 276         if (fr >= nffr && (outh || bLowAll || fr % nf == nf - 1))
 277         {
 278             /* Lowpass filtering */
 279             for (j = 0; j < nat; j++)
 280             {
 281                 clear_rvec(xf[j]);
 282             }
 283             clear_mat(boxf);
 284             for (i = 0; i < nffr; i++)
 285             {
 286                 for (j = 0; j < nat; j++)
 287                 {
 288                     for (d = 0; d < DIM; d++)
 289                     {
 290                         xf[j][d] += filt[i]*x[i][j][d];
 291                     }
 292                 }
 293                 for (j = 0; j < DIM; j++)
 294                 {
 295                     for (d = 0; d < DIM; d++)
 296                     {
 297                         boxf[j][d] += filt[i]*box[i][j][d];
 298                     }
 299                 }
 300             }
 301             if (outl && (bLowAll || fr % nf == nf - 1))
 302             {
 303                 write_trx(outl, nat, ind, topfile ? &(top.atoms) : NULL,
 304                           0, t[nf - 1], bFit ? topbox : boxf, xf, NULL, NULL);
 305             }
 306             if (outh)
 307             {
 308                 /* Highpass filtering */
 309                 for (j = 0; j < nat; j++)
 310                 {
 311                     for (d = 0; d < DIM; d++)
 312                     {
 313                         xf[j][d] = xtop[j][d] + x[nf - 1][j][d] - xf[j][d];
 314                     }
 315                 }
 316                 if (bFit)
 317                 {
 318                     for (j = 0; j < nat; j++)
 319                     {
 320                         rvec_inc(xf[j], xcmtop);
 321                     }
 322                 }
 323                 for (j = 0; j < DIM; j++)
 324                 {
 325                     for (d = 0; d < DIM; d++)
 326                     {
 327                         boxf[j][d] = topbox[j][d] + box[nf - 1][j][d] - boxf[j][d];
 328                     }
 329                 }
 330                 write_trx(outh, nat, ind, topfile ? &(top.atoms) : NULL,
 331                           0, t[nf - 1], bFit ? topbox : boxf, xf, NULL, NULL);
 332             }
 333         }
 334         /* Cycle all the pointer and the box by one */
 335         ptr = x[0];
 336         for (i = 0; i < nffr-1; i++)
 337         {
 338             t[i] = t[i+1];
 339             x[i] = x[i+1];
 340             copy_mat(box[i+1], box[i]);
 341         }
 342         x[nffr - 1] = ptr;
 343         fr++;
 344     }
 345     while (read_next_x(oenv, in, &(t[nffr - 1]), x[nffr - 1], box[nffr - 1]));
 346
 347     if (bTop)
 348     {
 349         gmx_rmpbc_done(gpbc);
 350     }
 351
 352     if (outh)
 353     {
 354         close_trx(outh);
 355     }
 356     if (outl)
 357     {
 358         close_trx(outl);
 359     }
 360     close_trx(in);
 361
 362     return 0;
 363 }