src/tools/proptim.c

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