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