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