src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "statutil.h"
  45 #include "sysstuff.h"
  46 #include "typedefs.h"
  47 #include "smalloc.h"
  48 #include "macros.h"
  49 #include "vec.h"
  50 #include "statutil.h"
  51 #include "index.h"
  52 #include "tpxio.h"
  53 #include "princ.h"
  54 #include "do_fit.h"
  55 #include "copyrite.h"
  56 #include "rmpbc.h"
  57 #include "gmx_ana.h"
  58
  59
  60 int gmx_filter(int argc,char *argv[])
  61 {
  62   const char *desc[] = {
  63     "[TT]g_filter[tt] performs frequency filtering on a trajectory.",
  64     "The filter shape is cos([GRK]pi[grk] t/A) + 1 from -A to +A, where A is given",
  65     "by the option [TT]-nf[tt] times the time step in the input trajectory.",
  66     "This filter reduces fluctuations with period A by 85%, with period",
  67     "2*A by 50% and with period 3*A by 17% for low-pass filtering.",
  68     "Both a low-pass and high-pass filtered trajectory can be written.[PAR]",
  69
  70     "Option [TT]-ol[tt] writes a low-pass filtered trajectory.",
  71     "A frame is written every [TT]-nf[tt] input frames.",
  72     "This ratio of filter length and output interval ensures a good",
  73     "suppression of aliasing of high-frequency motion, which is useful for",
  74     "making smooth movies. Also averages of properties which are linear",
  75     "in the coordinates are preserved, since all input frames are weighted",
  76     "equally in the output.",
  77     "When all frames are needed, use the [TT]-all[tt] option.[PAR]",
  78
  79     "Option [TT]-oh[tt] writes a high-pass filtered trajectory.",
  80     "The high-pass filtered coordinates are added to the coordinates",
  81     "from the structure file. When using high-pass filtering use [TT]-fit[tt]",
  82     "or make sure you use a trajectory that has been fitted on",
  83     "the coordinates in the structure file."
  84   };
  85
  86   static int nf=10;
  87   static gmx_bool bNoJump = TRUE,bFit = FALSE,bLowAll = FALSE;
  88   t_pargs pa[] = {
  89     { "-nf", FALSE, etINT, {&nf},
  90       "Sets the filter length as well as the output interval for low-pass filtering" },
  91     { "-all", FALSE, etBOOL, {&bLowAll},
  92       "Write all low-pass filtered frames" },
  93     { "-nojump", FALSE, etBOOL, {&bNoJump},
  94       "Remove jumps of atoms across the box" },
  95     { "-fit", FALSE, etBOOL, {&bFit},
  96       "Fit all frames to a reference structure" }
  97   };
  98   const char *topfile,*lowfile,*highfile;
  99   gmx_bool       bTop=FALSE;
 100   t_topology top;
 101   int        ePBC=-1;
 102   rvec       *xtop;
 103   matrix     topbox,*box,boxf;
 104   char       title[256],*grpname;
 105   int        isize;
 106   atom_id    *index;
 107   real       *w_rls=NULL;
 108   t_trxstatus *in;
 109   t_trxstatus *outl,*outh;
 110   int        nffr,i,fr,nat,j,d,m;
 111   atom_id    *ind;
 112   real       flen,*filt,sum,*t;
 113   rvec       xcmtop,xcm,**x,*ptr,*xf,*xn,*xp,hbox;
 114   output_env_t oenv;
 115   gmx_rmpbc_t  gpbc=NULL;
 116
 117 #define NLEG asize(leg)
 118   t_filenm fnm[] = {
 119     { efTRX, "-f", NULL, ffREAD  },
 120     { efTPS, NULL, NULL, ffOPTRD },
 121     { efNDX, NULL, NULL, ffOPTRD },
 122     { efTRO, "-ol", "lowpass",  ffOPTWR },
 123     { efTRO, "-oh", "highpass", ffOPTWR }
 124   };
 125 #define NFILE asize(fnm)
 126
 127   CopyRight(stderr,argv[0]);
 128   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 129                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 130
 131   highfile = opt2fn_null("-oh",NFILE,fnm);
 132   if (highfile) {
 133     topfile = ftp2fn(efTPS,NFILE,fnm);
 134     lowfile = opt2fn_null("-ol",NFILE,fnm);
 135   } else {
 136     topfile = ftp2fn_null(efTPS,NFILE,fnm);
 137     lowfile = opt2fn("-ol",NFILE,fnm);
 138   }
 139   if (topfile) {
 140     bTop = read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,
 141                          &xtop,NULL,topbox,TRUE);
 142     if (bTop) {
 143       gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,topbox);
 144       gmx_rmpbc(gpbc,top.atoms.nr,topbox,xtop);
 145     }
 146   }
 147
 148   clear_rvec(xcmtop);
 149   if (bFit) {
 150     fprintf(stderr,"Select group for least squares fit\n");
 151     get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 152     /* Set the weight */
 153     snew(w_rls,top.atoms.nr);
 154     for(i=0; i<isize; i++)
 155       w_rls[index[i]] = top.atoms.atom[index[i]].m;
 156     calc_xcm(xtop,isize,index,top.atoms.atom,xcmtop,FALSE);
 157     for(j=0; j<top.atoms.nr; j++)
 158       rvec_dec(xtop[j],xcmtop);
 159   }
 160
 161   /* The actual filter length flen can actually be any real number */
 162   flen = 2*nf;
 163   /* nffr is the number of frames that we filter over */
 164   nffr = 2*nf - 1;
 165   snew(filt,nffr);
 166   sum = 0;
 167   for(i=0; i<nffr; i++) {
 168     filt[i] = cos(2*M_PI*(i - nf + 1)/(real)flen) + 1;
 169     sum += filt[i];
 170   }
 171   fprintf(stdout,"filter weights:");
 172   for(i=0; i<nffr; i++) {
 173     filt[i] /= sum;
 174     fprintf(stdout," %5.3f",filt[i]);
 175   }
 176   fprintf(stdout,"\n");
 177
 178   snew(t,nffr);
 179   snew(x,nffr);
 180   snew(box,nffr);
 181
 182   nat = read_first_x(oenv,&in,opt2fn("-f",NFILE,fnm),
 183                      &(t[nffr - 1]),&(x[nffr - 1]),box[nffr - 1]);
 184   snew(ind,nat);
 185   for(i=0; i<nat; i++)
 186     ind[i] = i;
 187   /* x[nffr - 1] was already allocated by read_first_x */
 188   for(i=0; i<nffr-1; i++)
 189     snew(x[i],nat);
 190   snew(xf,nat);
 191   if (lowfile)
 192     outl = open_trx(lowfile,"w");
 193   else
 194     outl = 0;
 195   if (highfile)
 196     outh = open_trx(highfile,"w");
 197   else
 198     outh = 0;
 199
 200   fr = 0;
 201   do {
 202     xn = x[nffr - 1];
 203     if (bNoJump && fr > 0) {
 204       xp = x[nffr - 2];
 205       for(j=0; j<nat; j++)
 206         for(d=0; d<DIM; d++)
 207           hbox[d] = 0.5*box[nffr - 1][d][d];
 208       for(i=0; i<nat; i++)
 209         for(m=DIM-1; m>=0; m--)
 210           if (hbox[m] > 0) {
 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             while (xn[i][m] - xp[i][m] > hbox[m])
 215               for(d=0; d<=m; d++)
 216                 xn[i][d] -= box[nffr - 1][m][d];
 217           }
 218     }
 219     if (bTop) {
 220       gmx_rmpbc(gpbc,nat,box[nffr - 1],xn);
 221     }
 222     if (bFit) {
 223       calc_xcm(xn,isize,index,top.atoms.atom,xcm,FALSE);
 224       for(j=0; j<nat; j++)
 225         rvec_dec(xn[j],xcm);
 226       do_fit(nat,w_rls,xtop,xn);
 227       for(j=0; j<nat; j++)
 228         rvec_inc(xn[j],xcmtop);
 229     }
 230     if (fr >= nffr && (outh || bLowAll || fr % nf == nf - 1)) {
 231       /* Lowpass filtering */
 232       for(j=0; j<nat; j++)
 233         clear_rvec(xf[j]);
 234       clear_mat(boxf);
 235       for(i=0; i<nffr; i++) {
 236         for(j=0; j<nat; j++)
 237           for(d=0; d<DIM; d++)
 238             xf[j][d] += filt[i]*x[i][j][d];
 239         for(j=0; j<DIM; j++)
 240           for(d=0; d<DIM; d++)
 241             boxf[j][d] += filt[i]*box[i][j][d];
 242       }
 243       if (outl && (bLowAll || fr % nf == nf - 1))
 244         write_trx(outl,nat,ind,topfile ? &(top.atoms) : NULL,
 245                   0,t[nf - 1],bFit ? topbox : boxf,xf,NULL,NULL);
 246       if (outh) {
 247         /* Highpass filtering */
 248         for(j=0; j<nat; j++)
 249           for(d=0; d<DIM; d++)
 250             xf[j][d] = xtop[j][d] + x[nf - 1][j][d] - xf[j][d];
 251         if (bFit)
 252           for(j=0; j<nat; j++)
 253             rvec_inc(xf[j],xcmtop);
 254         for(j=0; j<DIM; j++)
 255           for(d=0; d<DIM; d++)
 256             boxf[j][d] = topbox[j][d] + box[nf - 1][j][d] - boxf[j][d];
 257         write_trx(outh,nat,ind,topfile ? &(top.atoms) : NULL,
 258                   0,t[nf - 1],bFit ? topbox : boxf,xf,NULL,NULL);
 259       }
 260     }
 261     /* Cycle all the pointer and the box by one */
 262     ptr = x[0];
 263     for(i=0; i<nffr-1; i++) {
 264       t[i] = t[i+1];
 265       x[i] = x[i+1];
 266       copy_mat(box[i+1],box[i]);
 267     }
 268     x[nffr - 1] = ptr;
 269     fr++;
 270   } while (read_next_x(oenv,in,&(t[nffr - 1]),nat,x[nffr - 1],box[nffr - 1]));
 271
 272   if (bTop)
 273     gmx_rmpbc_done(gpbc);
 274
 275   if (outh)
 276     close_trx(outh);
 277   if (outl)
 278     close_trx(outl);
 279   close_trx(in);
 280
 281   return 0;
 282 }