src/tools/mk_angndx.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_mk_angndx_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "smalloc.h"
  33 #include "copyrite.h"
  34 #include "statutil.h"
  35 #include "macros.h"
  36 #include "string2.h"
  37 #include "futil.h"
  38 #include "assert.h"
  39 #include "fatal.h"
  40
  41 static int calc_nftype(int FTYPE,t_idef *idef)
  42 {
  43   int i,nf=0;
  44
  45   for(i=0; (i<idef->ntypes); i++)
  46     if (idef->functype[i] == FTYPE)
  47       nf++;
  48
  49   return nf;
  50 }
  51
  52 static void fill_ft_ind(int FTYPE,int ft_ind[],t_idef *idef,char *grpnames[])
  53 {
  54   char buf[125];
  55   int  i,ft,ind=0;
  56
  57   for(i=0; (i<idef->ntypes); i++) {
  58     ft=idef->functype[i];
  59     if (ft == FTYPE) {
  60       ft_ind[i]=ind;
  61       switch (FTYPE) {
  62       case F_G96ANGLES:
  63         sprintf(buf,"Theta=%.1f_%g",idef->iparams[i].harmonic.rA,
  64                 idef->iparams[i].harmonic.krA);
  65         break;
  66       case F_ANGLES:
  67         sprintf(buf,"Theta=%.1f_%g",idef->iparams[i].harmonic.rA,
  68                 idef->iparams[i].harmonic.krA);
  69         break;
  70       case F_PDIHS:
  71         sprintf(buf,"Phi=%.1f_%d_%g",idef->iparams[i].pdihs.phiA,
  72                 idef->iparams[i].pdihs.mult,idef->iparams[i].pdihs.cpA);
  73         break;
  74       case F_IDIHS:
  75         sprintf(buf,"Xi=%.1f_%g",idef->iparams[i].harmonic.rA,
  76                 idef->iparams[i].harmonic.krA);
  77         break;
  78       case F_RBDIHS:
  79         sprintf(buf,"RB-Dihs");
  80         break;
  81       default:
  82         fatal_error(0,"kjdshfgkajhgajytgajtgasuyf");
  83       }
  84       grpnames[ind]=strdup(buf);
  85       ind++;
  86     }
  87     else
  88       ft_ind[i]=-1;
  89   }
  90 }
  91
  92 static void fill_ang(int FTYPE,int fac,
  93                      int nr[],int *index[],int ft_ind[],t_idef *idef)
  94 {
  95   int     i,j,ft,fft,nr_fac;
  96   t_iatom *ia;
  97
  98   ia=idef->il[FTYPE].iatoms;
  99   for(i=0; (i<idef->il[FTYPE].nr); ) {
 100     ft  = idef->functype[ia[0]];
 101     fft = ft_ind[ia[0]];
 102     assert(fft != -1);
 103     nr_fac=fac*nr[fft];
 104     for(j=0; (j<fac); j++)
 105       index[fft][nr_fac+j]=ia[j+1];
 106     nr[fft]++;
 107     ia += interaction_function[ft].nratoms+1;
 108     i  += interaction_function[ft].nratoms+1;
 109   }
 110 }
 111
 112 int main(int argc,char *argv[])
 113 {
 114   static char *desc[] = {
 115     "mk_angndx makes an index file for calculation of",
 116     "angle distributions etc. It uses a run input file ([TT].tpx[tt]) for the",
 117     "definitions of the angles, dihedrals etc."
 118   };
 119   static char *opt[] = { NULL, "angle", "g96-angle", "dihedral", "improper", "ryckaert-bellemans", "phi-psi", NULL };
 120   t_pargs pa[] = {
 121     { "-type", FALSE, etENUM, {opt},
 122       "Type of angle" }
 123   };
 124
 125   FILE       *out;
 126   t_topology *top;
 127   int        i,j,nftype,nang;
 128   int        mult=0,FTYPE=0;
 129   bool       bPP;
 130   int        **index;
 131   int        *ft_ind;
 132   int        *nr;
 133   char       **grpnames;
 134   t_filenm fnm[] = {
 135     { efTPX, NULL, NULL, ffREAD  },
 136     { efNDX, NULL, "angle", ffWRITE }
 137   };
 138 #define NFILE asize(fnm)
 139
 140   CopyRight(stderr,argv[0]);
 141   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,
 142                     asize(desc),desc,0,NULL);
 143
 144   bPP  = FALSE;
 145   mult = 4;
 146   switch(opt[0][0]) {
 147   case 'a':
 148     mult=3;
 149     FTYPE=F_ANGLES;
 150     break;
 151   case 'g':
 152     mult=3;
 153     FTYPE=F_G96ANGLES;
 154     break;
 155   case 'd':
 156     FTYPE=F_PDIHS;
 157     break;
 158   case 'i':
 159     FTYPE=F_IDIHS;
 160     break;
 161   case 'r':
 162     FTYPE=F_RBDIHS;
 163     break;
 164   case 'p':
 165     bPP=TRUE;
 166     break;
 167   default:
 168     break;
 169   }
 170
 171   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 172
 173   if (!bPP) {
 174     nftype=calc_nftype(FTYPE,&(top->idef));
 175     snew(grpnames,nftype);
 176     snew(ft_ind,top->idef.ntypes);
 177     fill_ft_ind(FTYPE,ft_ind,&top->idef,grpnames);
 178     nang=top->idef.il[FTYPE].nr;
 179
 180     snew(nr,nftype);
 181     snew(index,nftype);
 182     for(i=0; (i<nftype); i++)
 183       snew(index[i],nang*mult);
 184
 185     fill_ang(FTYPE,mult,nr,index,ft_ind,&(top->idef));
 186
 187     out=ftp2FILE(efNDX,NFILE,fnm,"w");
 188     for(i=0; (i<nftype); i++) {
 189       fprintf(out,"[ %s ]\n",grpnames[i]);
 190       for(j=0; (j<nr[i]*mult); j++) {
 191         fprintf(out," %5d",index[i][j]+1);
 192         if ((j % 12) == 11)
 193           fprintf(out,"\n");
 194       }
 195       fprintf(out,"\n");
 196     }
 197     fclose(out);
 198   }
 199   else {
 200     fprintf(stderr,"Sorry, maybe later...\n");
 201   }
 202
 203   thanx(stdout);
 204
 205   return 0;
 206 }