src/tools/recomb.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_recomb_c = "$Id$";
  30
  31 #include "recomb.h"
  32 #include "futil.h"
  33 #include "wgms.h"
  34 #include "smalloc.h"
  35
  36 real *read_gammaf(char *fn,int nframes)
  37 {
  38   FILE   *in;
  39   real   *gf;
  40   double y;
  41   int    i;
  42
  43   snew(gf,nframes);
  44   in=ffopen(fn,"r");
  45   for(i=0; (i<nframes); i++) {
  46     fscanf(in,"%lf",&y);
  47     gf[i]=y;
  48   }
  49   fclose(in);
  50   fprintf(stderr,"Succesfully read gamma\n");
  51   return gf;
  52 }
  53
  54 void recombine(char *base,char *gammaf,int nskip,
  55                int nframes,int nev,int natoms,
  56                rvec *ev[],real *evprj[],
  57                rvec yav[],atom_id all_index[])
  58 {
  59   static char *format=
  60     "Recombined projection of Gamma trj (EV %d) in Cartesian Space\n";
  61   FILE *out;
  62   rvec *xxx,*evptr;
  63   real *gamma;
  64   real prj;
  65   char buf[256];
  66   int  i,j,n;
  67   real gt;
  68
  69   gamma=read_gammaf(gammaf,nframes);
  70   snew(xxx,natoms);
  71   for(n=0; (n<nev); n++) {
  72     sprintf(buf,"%s%d",base,n+1);
  73     out=ffopen(buf,"w");
  74     fprintf(out,format,n+1);
  75     fprintf(stderr,format,n+1);
  76     evptr=ev[n];
  77
  78     for(j=0; (j<nframes); j++) {
  79       if ((j % 50) == 0)
  80         fprintf(stderr,"\r frame %d",j);
  81       if ((nskip == 0) || ((j % nskip) == 0)) {
  82         gt=1.0/gamma[j];
  83         prj=evprj[n][j];
  84         for(i=0; (i<natoms); i++) {
  85           xxx[i][XX]=(yav[i][XX]+prj*evptr[i][XX])*gt;
  86           xxx[i][YY]=(yav[i][YY]+prj*evptr[i][YY])*gt;
  87           xxx[i][ZZ]=(yav[i][ZZ]+prj*evptr[i][ZZ])*gt;
  88         }
  89         write_gms_ndx(out,natoms,all_index,xxx,NULL);
  90       }
  91     }
  92     fclose(out);
  93     fprintf(stderr,"\r");
  94   }
  95   fprintf(stderr,"\n");
  96   sfree(xxx);
  97   sfree(gamma);
  98 }