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