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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Green Red Orange Magenta Azure Cyan Skyblue
  31  */
  32
  33 #include "typedefs.h"
  34 #include "maths.h"
  35 #include "string2.h"
  36 #include "hxprops.h"
  37 #include "fatal.h"
  38 #include "futil.h"
  39 #include "smalloc.h"
  40 #include "readev.h"
  41 #include "macros.h"
  42 #include "confio.h"
  43 #include "copyrite.h"
  44 #include "statutil.h"
  45 #include "orise.h"
  46 #include "pinput.h"
  47 #include "recomb.h"
  48
  49 void calc_prj(int natoms,rvec xav[],rvec x[],rvec ev[],real eprj)
  50 {
  51   int  i,m;
  52
  53   for(i=0; (i<natoms); i++) {
  54     for(m=0; (m<DIM); m++)
  55       x[i][m]=xav[i][m]+ev[i][m]*eprj;
  56   }
  57 }
  58
  59 void mkptrj(char *prop,int nSel,
  60             int natoms,rvec xav[],int nframes,
  61             int nev,rvec **EV,real **evprj,
  62             int nca,atom_id ca_index[],atom_id bb_index[],
  63             t_atom atom[],matrix box)
  64 {
  65   FILE       *out;
  66   char       buf[256];
  67   double     propje,*pav,*pav2;
  68   int        i,j;
  69   rvec       *xxx;
  70
  71   snew(pav,nev);
  72   snew(pav2,nev);
  73   snew(xxx,natoms);
  74
  75   sprintf(buf,"%s.trj1",prop);
  76   out=ffopen(buf,"w");
  77   fprintf(out,"Projection of %s on EigenVectors\n",prop);
  78   for(j=0; (j<nframes); j++) {
  79     if ((j % 10) == 0)
  80       fprintf(stderr,"\rFrame %d",j);
  81     for(i=0; (i<nev); i++) {
  82       calc_prj(natoms,xav,xxx,EV[i],evprj[i][j]);
  83       switch (nSel) {
  84       case efhRAD:
  85         propje=radius(NULL,nca,ca_index,xxx);
  86         break;
  87       case efhLEN:
  88         propje=ahx_len(nca,ca_index,xxx,box);
  89         break;
  90       case efhTWIST:
  91         propje=twist(NULL,nca,ca_index,xxx);
  92         break;
  93       case efhCPHI:
  94         propje=ca_phi(nca,ca_index,xxx);
  95         break;
  96       case efhDIP:
  97         propje=dip(natoms,bb_index,xxx,atom);
  98         break;
  99       default:
 100         fatal_error(0,"Not implemented");
 101       }
 102       pav[i]+=propje;
 103       pav2[i]+=propje*propje;
 104       fprintf(out,"%8.3f",propje);
 105     }
 106     fprintf(out,"\n");
 107   }
 108   fclose(out);
 109   fprintf(stderr,"\n");
 110   for(i=0; (i<nev); i++) {
 111     printf("ev %2d, average: %8.3f  rms: %8.3f\n",
 112            i+1,pav[i]/nframes,sqrt(pav2[i]/nframes-sqr(pav[i]/nframes)));
 113   }
 114   sfree(pav);
 115   sfree(pav2);
 116   sfree(xxx);
 117 }
 118
 119 void proptrj(char *fngro,char *fndat,t_topology *top,t_pinp *p)
 120 {
 121   static char *ppp[efhNR] = {
 122     "RAD", "TWIST", "RISE", "LEN", "NHX", "DIP", "RMS", "CPHI",
 123     "RMSA", "PHI", "PSI", "HB3", "HB4", "HB5"
 124   };
 125   FILE       *out;
 126   rvec       **EV;
 127   real       **evprj;
 128   atom_id    *index;
 129   int        natoms,nca,nSel,nframes,nev;
 130   rvec       *xav,*vav;
 131   atom_id    *ca_index,*bb_index;
 132   matrix     box;
 133   char       buf[256],*prop;
 134   double     x;
 135   int        i,j,d;
 136
 137   nframes = p->nframes;
 138   nSel    = p->nSel;
 139   nev     = p->nev;
 140
 141   prop=ppp[nSel];
 142   evprj=read_proj(nev,nframes,p->base);
 143
 144   get_coordnum(fngro,&natoms);
 145   snew(xav,natoms);
 146   snew(vav,natoms);
 147   read_conf(fngro,buf,&natoms,xav,vav,box);
 148   fprintf(stderr,"Succesfully read average positions (%s)\n",buf);
 149
 150   EV=read_ev(fndat,natoms);
 151
 152   fprintf(stderr,"Succesfully read eigenvectors\n");
 153
 154   snew(index,nev);
 155   for(i=0; (i<nev); i++)
 156     index[i]=i;
 157   snew(bb_index,natoms);
 158   for(i=0; (i<natoms); i++)
 159     bb_index[i]=i;
 160   snew(ca_index,natoms);
 161   for(i=nca=0; (i<natoms); i++)
 162     if ((strcmp("CA",*(top->atoms.atomname[i])) == 0))
 163       ca_index[nca++]=i;
 164
 165   switch (p->funct) {
 166   case ptMC:
 167     switch(nSel) {
 168     case efhRAD:
 169       optim_radius(nev,xav,EV,evprj,natoms,nca,ca_index,p);
 170       break;
 171     case efhRISE:
 172       optim_rise(nev,xav,EV,evprj,natoms,nca,ca_index,p);
 173       break;
 174     default:
 175       break;
 176     }
 177     break;
 178   case ptREC:
 179     recombine(p->recomb,p->gamma,p->nskip,nframes,nev,natoms,
 180               EV,evprj,xav,bb_index);
 181     break;
 182   case ptPTRJ:
 183     mkptrj(prop,nSel,natoms,xav,nframes,
 184            nev,EV,evprj,nca,ca_index,bb_index,
 185            top->atoms.atom,box);
 186     break;
 187   default:
 188     fatal_error(0,"I Don't Know What to Do");
 189   }
 190 }
 191
 192 int main(int argc,char *argv[])
 193 {
 194   static char *desc[] = {
 195     "proptrj"
 196   };
 197   t_manual man = { asize(desc),desc,0,NULL,NULL,0,NULL};
 198
 199   t_filenm  fnm[] = {
 200     { efGRO, "-c", "aver",FALSE },
 201     { efDAT, "-d", "eigenvec", FALSE },
 202     { efTPX, NULL, NULL, FALSE },
 203     { efDAT, "-pi","pinp", FALSE },
 204     { efDAT, "-po","poutp", FALSE }
 205   };
 206 #define NFILE asize(fnm)
 207   t_topology *top;
 208   t_pinp     *p;
 209
 210   CopyRight(stderr,argv[0]);
 211   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,
 212                     NFILE,fnm,TRUE,&man);
 213
 214   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 215   init_debug("proptim.dbg",0);
 216   snew(p,1);
 217   read_inp(opt2fn("-pi",NFILE,fnm),opt2fn("-po",NFILE,fnm),p);
 218
 219   proptrj(ftp2fn(efGRO,NFILE,fnm),ftp2fn(efDAT,NFILE,fnm),top,p);
 220
 221   thanx(stderr);
 222
 223   return 0;
 224 }
 225