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