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