src/gromacs/gmxana/gmx_filter.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <math.h>
  39 #include <string.h>
  40
  41 #include "statutil.h"
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "vec.h"
  47 #include "statutil.h"
  48 #include "index.h"
  49 #include "tpxio.h"
  50 #include "princ.h"
  51 #include "do_fit.h"
  52 #include "rmpbc.h"
  53 #include "gmx_ana.h"
  54
  55
  56 int gmx_filter(int argc, char *argv[])
  57 {
  58     const char     *desc[] = {
  59         "[TT]g_filter[tt] 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 | PCA_BE_NICE,
 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 }