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