src/tools/proptim.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_proptim_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "maths.h"
  33 #include "string2.h"
  34 #include "hxprops.h"
  35 #include "fatal.h"
  36 #include "futil.h"
  37 #include "smalloc.h"
  38 #include "readev.h"
  39 #include "macros.h"
  40 #include "confio.h"
  41 #include "copyrite.h"
  42 #include "statutil.h"
  43 #include "orise.h"
  44 #include "pinput.h"
  45 #include "recomb.h"
  46
  47 void calc_prj(int natoms,rvec xav[],rvec x[],rvec ev[],real eprj)
  48 {
  49   int  i,m;
  50
  51   for(i=0; (i<natoms); i++) {
  52     for(m=0; (m<DIM); m++)
  53       x[i][m]=xav[i][m]+ev[i][m]*eprj;
  54   }
  55 }
  56
  57 void mkptrj(char *prop,int nSel,
  58             int natoms,rvec xav[],int nframes,
  59             int nev,rvec **EV,real **evprj,
  60             int nca,atom_id ca_index[],atom_id bb_index[],
  61             t_atom atom[],matrix box)
  62 {
  63   FILE       *out;
  64   char       buf[256];
  65   double     propje,*pav,*pav2;
  66   int        i,j;
  67   rvec       *xxx;
  68
  69   snew(pav,nev);
  70   snew(pav2,nev);
  71   snew(xxx,natoms);
  72
  73   sprintf(buf,"%s.trj1",prop);
  74   out=ffopen(buf,"w");
  75   fprintf(out,"Projection of %s on EigenVectors\n",prop);
  76   for(j=0; (j<nframes); j++) {
  77     if ((j % 10) == 0)
  78       fprintf(stderr,"\rFrame %d",j);
  79     for(i=0; (i<nev); i++) {
  80       calc_prj(natoms,xav,xxx,EV[i],evprj[i][j]);
  81       switch (nSel) {
  82       case efhRAD:
  83         propje=radius(NULL,nca,ca_index,xxx);
  84         break;
  85       case efhLEN:
  86         propje=ahx_len(nca,ca_index,xxx,box);
  87         break;
  88       case efhTWIST:
  89         propje=twist(NULL,nca,ca_index,xxx);
  90         break;
  91       case efhCPHI:
  92         propje=ca_phi(nca,ca_index,xxx);
  93         break;
  94       case efhDIP:
  95         propje=dip(natoms,bb_index,xxx,atom);
  96         break;
  97       default:
  98         fatal_error(0,"Not implemented");
  99       }
 100       pav[i]+=propje;
 101       pav2[i]+=propje*propje;
 102       fprintf(out,"%8.3f",propje);
 103     }
 104     fprintf(out,"\n");
 105   }
 106   fclose(out);
 107   fprintf(stderr,"\n");
 108   for(i=0; (i<nev); i++) {
 109     printf("ev %2d, average: %8.3f  rms: %8.3f\n",
 110            i+1,pav[i]/nframes,sqrt(pav2[i]/nframes-sqr(pav[i]/nframes)));
 111   }
 112   sfree(pav);
 113   sfree(pav2);
 114   sfree(xxx);
 115 }
 116
 117 void proptrj(char *fngro,char *fndat,t_topology *top,t_pinp *p)
 118 {
 119   static char *ppp[efhNR] = {
 120     "RAD", "TWIST", "RISE", "LEN", "NHX", "DIP", "RMS", "CPHI",
 121     "RMSA", "PHI", "PSI", "HB3", "HB4", "HB5"
 122   };
 123   FILE       *out;
 124   rvec       **EV;
 125   real       **evprj;
 126   atom_id    *index;
 127   int        natoms,nca,nSel,nframes,nev;
 128   rvec       *xav,*vav;
 129   atom_id    *ca_index,*bb_index;
 130   matrix     box;
 131   char       buf[256],*prop;
 132   double     x;
 133   int        i,j,d;
 134
 135   nframes = p->nframes;
 136   nSel    = p->nSel;
 137   nev     = p->nev;
 138
 139   prop=ppp[nSel];
 140   evprj=read_proj(nev,nframes,p->base);
 141
 142   get_coordnum(fngro,&natoms);
 143   snew(xav,natoms);
 144   snew(vav,natoms);
 145   read_conf(fngro,buf,&natoms,xav,vav,box);
 146   fprintf(stderr,"Succesfully read average positions (%s)\n",buf);
 147
 148   EV=read_ev(fndat,natoms);
 149
 150   fprintf(stderr,"Succesfully read eigenvectors\n");
 151
 152   snew(index,nev);
 153   for(i=0; (i<nev); i++)
 154     index[i]=i;
 155   snew(bb_index,natoms);
 156   for(i=0; (i<natoms); i++)
 157     bb_index[i]=i;
 158   snew(ca_index,natoms);
 159   for(i=nca=0; (i<natoms); i++)
 160     if ((strcmp("CA",*(top->atoms.atomname[i])) == 0))
 161       ca_index[nca++]=i;
 162
 163   switch (p->funct) {
 164   case ptMC:
 165     switch(nSel) {
 166     case efhRAD:
 167       optim_radius(nev,xav,EV,evprj,natoms,nca,ca_index,p);
 168       break;
 169     case efhRISE:
 170       optim_rise(nev,xav,EV,evprj,natoms,nca,ca_index,p);
 171       break;
 172     default:
 173       break;
 174     }
 175     break;
 176   case ptREC:
 177     recombine(p->recomb,p->gamma,p->nskip,nframes,nev,natoms,
 178               EV,evprj,xav,bb_index);
 179     break;
 180   case ptPTRJ:
 181     mkptrj(prop,nSel,natoms,xav,nframes,
 182            nev,EV,evprj,nca,ca_index,bb_index,
 183            top->atoms.atom,box);
 184     break;
 185   default:
 186     fatal_error(0,"I Don't Know What to Do");
 187   }
 188 }
 189
 190 int main(int argc,char *argv[])
 191 {
 192   static char *desc[] = {
 193     "proptrj"
 194   };
 195   t_manual man = { asize(desc),desc,0,NULL,NULL,0,NULL};
 196
 197   t_filenm  fnm[] = {
 198     { efGRO, "-c", "aver",FALSE },
 199     { efDAT, "-d", "eigenvec", FALSE },
 200     { efTPX, NULL, NULL, FALSE },
 201     { efDAT, "-pi","pinp", FALSE },
 202     { efDAT, "-po","poutp", FALSE }
 203   };
 204 #define NFILE asize(fnm)
 205   t_topology *top;
 206   t_pinp     *p;
 207
 208   CopyRight(stderr,argv[0]);
 209   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,
 210                     NFILE,fnm,TRUE,&man);
 211
 212   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 213   init_debug("proptim.dbg",0);
 214   snew(p,1);
 215   read_inp(opt2fn("-pi",NFILE,fnm),opt2fn("-po",NFILE,fnm),p);
 216
 217   proptrj(ftp2fn(efGRO,NFILE,fnm),ftp2fn(efDAT,NFILE,fnm),top,p);
 218
 219   thanx(stdout);
 220
 221   return 0;
 222 }
 223