src/tools/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     "g_filter performs frequency filtering on a trajectory.",
  61     "The filter shape is cos(pi 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 which has been fitted on",
  80     "the coordinates in the structure file."
  81   };
  82
  83   static int nf=10;
  84   static 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   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   int        in,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
 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   CopyRight(stderr,argv[0]);
 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     topfile = ftp2fn(efTPS,NFILE,fnm);
 129     lowfile = opt2fn_null("-ol",NFILE,fnm);
 130   } else {
 131     topfile = ftp2fn_null(efTPS,NFILE,fnm);
 132     lowfile = opt2fn("-ol",NFILE,fnm);
 133   }
 134
 135   if (topfile) {
 136     bTop = read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,
 137                          &xtop,NULL,topbox,TRUE);
 138     if (bTop) {
 139       rm_pbc(&(top.idef),ePBC,top.atoms.nr,topbox,xtop,xtop);
 140     }
 141   }
 142
 143   clear_rvec(xcmtop);
 144   if (bFit) {
 145     fprintf(stderr,"Select group for least squares fit\n");
 146     get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 147     /* Set the weight */
 148     snew(w_rls,top.atoms.nr);
 149     for(i=0; i<isize; i++)
 150       w_rls[index[i]] = top.atoms.atom[index[i]].m;
 151     calc_xcm(xtop,isize,index,top.atoms.atom,xcmtop,FALSE);
 152     for(j=0; j<top.atoms.nr; j++)
 153       rvec_dec(xtop[j],xcmtop);
 154   }
 155
 156   /* The actual filter length flen can actually be any real number */
 157   flen = 2*nf;
 158   /* nffr is the number of frames that we filter over */
 159   nffr = 2*nf - 1;
 160   snew(filt,nffr);
 161   sum = 0;
 162   for(i=0; i<nffr; i++) {
 163     filt[i] = cos(2*M_PI*(i - nf + 1)/(real)flen) + 1;
 164     sum += filt[i];
 165   }
 166   fprintf(stdout,"filter weights:");
 167   for(i=0; i<nffr; i++) {
 168     filt[i] /= sum;
 169     fprintf(stdout," %5.3f",filt[i]);
 170   }
 171   fprintf(stdout,"\n");
 172
 173   snew(t,nffr);
 174   snew(x,nffr);
 175   snew(box,nffr);
 176
 177   nat = read_first_x(oenv,&in,opt2fn("-f",NFILE,fnm),
 178                      &(t[nffr - 1]),&(x[nffr - 1]),box[nffr - 1]);
 179   snew(ind,nat);
 180   for(i=0; i<nat; i++)
 181     ind[i] = i;
 182   /* x[nffr - 1] was already allocated by read_first_x */
 183   for(i=0; i<nffr-1; i++)
 184     snew(x[i],nat);
 185   snew(xf,nat);
 186   if (lowfile)
 187     outl = open_trx(lowfile,"w");
 188   else
 189     outl = 0;
 190   if (highfile)
 191     outh = open_trx(highfile,"w");
 192   else
 193     outh = 0;
 194
 195   fr = 0;
 196   do {
 197     xn = x[nffr - 1];
 198     if (bNoJump && fr > 0) {
 199       xp = x[nffr - 2];
 200       for(j=0; j<nat; j++)
 201         for(d=0; d<DIM; d++)
 202           hbox[d] = 0.5*box[nffr - 1][d][d];
 203       for(i=0; i<nat; i++)
 204         for(m=DIM-1; m>=0; m--)
 205           if (hbox[m] > 0) {
 206             while (xn[i][m] - xp[i][m] <= -hbox[m])
 207               for(d=0; d<=m; d++)
 208                 xn[i][d] += box[nffr - 1][m][d];
 209             while (xn[i][m] - xp[i][m] > hbox[m])
 210               for(d=0; d<=m; d++)
 211                 xn[i][d] -= box[nffr - 1][m][d];
 212           }
 213     }
 214     if (bTop) {
 215       rm_pbc(&(top.idef),ePBC,nat,box[nffr - 1],xn,xn);
 216     }
 217     if (bFit) {
 218       calc_xcm(xn,isize,index,top.atoms.atom,xcm,FALSE);
 219       for(j=0; j<nat; j++)
 220         rvec_dec(xn[j],xcm);
 221       do_fit(nat,w_rls,xtop,xn);
 222       for(j=0; j<nat; j++)
 223         rvec_inc(xn[j],xcmtop);
 224     }
 225     if (fr >= nffr && (outh || bLowAll || fr % nf == nf - 1)) {
 226       /* Lowpass filtering */
 227       for(j=0; j<nat; j++)
 228         clear_rvec(xf[j]);
 229       clear_mat(boxf);
 230       for(i=0; i<nffr; i++) {
 231         for(j=0; j<nat; j++)
 232           for(d=0; d<DIM; d++)
 233             xf[j][d] += filt[i]*x[i][j][d];
 234         for(j=0; j<DIM; j++)
 235           for(d=0; d<DIM; d++)
 236             boxf[j][d] += filt[i]*box[i][j][d];
 237       }
 238       if (outl && (bLowAll || fr % nf == nf - 1))
 239         write_trx(outl,nat,ind,topfile ? &(top.atoms) : NULL,
 240                   0,t[nf - 1],bFit ? topbox : boxf,xf,NULL,NULL);
 241       if (outh) {
 242         /* Highpass filtering */
 243         for(j=0; j<nat; j++)
 244           for(d=0; d<DIM; d++)
 245             xf[j][d] = xtop[j][d] + x[nf - 1][j][d] - xf[j][d];
 246         if (bFit)
 247           for(j=0; j<nat; j++)
 248             rvec_inc(xf[j],xcmtop);
 249         for(j=0; j<DIM; j++)
 250           for(d=0; d<DIM; d++)
 251             boxf[j][d] = topbox[j][d] + box[nf - 1][j][d] - boxf[j][d];
 252         write_trx(outh,nat,ind,topfile ? &(top.atoms) : NULL,
 253                   0,t[nf - 1],bFit ? topbox : boxf,xf,NULL,NULL);
 254       }
 255     }
 256     /* Cycle all the pointer and the box by one */
 257     ptr = x[0];
 258     for(i=0; i<nffr-1; i++) {
 259       t[i] = t[i+1];
 260       x[i] = x[i+1];
 261       copy_mat(box[i+1],box[i]);
 262     }
 263     x[nffr - 1] = ptr;
 264     fr++;
 265   } while (read_next_x(oenv,in,&(t[nffr - 1]),nat,x[nffr - 1],box[nffr - 1]));
 266
 267   if (outh)
 268     close_trx(outh);
 269   if (outl)
 270     close_trx(outl);
 271   close_trx(in);
 272
 273   return 0;
 274 }