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 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 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "typedefs.h"
  41 #include "smalloc.h"
  42 #include "copyrite.h"
  43 #include "statutil.h"
  44 #include "macros.h"
  45 #include "string2.h"
  46 #include "futil.h"
  47 #include "gmx_fatal.h"
  48
  49 static int calc_ntype(int nft,int *ft,t_idef *idef)
  50 {
  51   int  i,f,nf=0;
  52
  53   for(i=0; (i<idef->ntypes); i++) {
  54     for(f=0; f<nft; f++) {
  55       if (idef->functype[i] == ft[f])
  56         nf++;
  57     }
  58   }
  59
  60   return nf;
  61 }
  62
  63 static void fill_ft_ind(int nft,int *ft,t_idef *idef,
  64                         int ft_ind[],char *grpnames[])
  65 {
  66   char buf[125];
  67   int  i,f,ftype,ind=0;
  68
  69   /* Loop over all the function types in the topology */
  70   for(i=0; (i<idef->ntypes); i++) {
  71     ft_ind[i] = -1;
  72     /* Check all the selected function types */
  73     for(f=0; f<nft; f++) {
  74       ftype = ft[f];
  75       if (idef->functype[i] == ftype) {
  76         ft_ind[i] = ind;
  77         switch (ftype) {
  78         case F_ANGLES:
  79           sprintf(buf,"Theta=%.1f_%.2f",idef->iparams[i].harmonic.rA,
  80                   idef->iparams[i].harmonic.krA);
  81           break;
  82         case F_G96ANGLES:
  83           sprintf(buf,"Cos_th=%.1f_%.2f",idef->iparams[i].harmonic.rA,
  84                   idef->iparams[i].harmonic.krA);
  85           break;
  86         case F_UREY_BRADLEY:
  87           sprintf(buf,"UB_th=%.1f_%.2f2f",idef->iparams[i].u_b.theta,
  88                   idef->iparams[i].u_b.ktheta);
  89           break;
  90         case F_QUARTIC_ANGLES:
  91           sprintf(buf,"Q_th=%.1f_%.2f_%.2f",idef->iparams[i].qangle.theta,
  92                   idef->iparams[i].qangle.c[0],idef->iparams[i].qangle.c[1]);
  93           break;
  94         case F_TABANGLES:
  95           sprintf(buf,"Table=%d_%.2f",idef->iparams[i].tab.table,
  96                   idef->iparams[i].tab.kA);
  97           break;
  98         case F_PDIHS:
  99           sprintf(buf,"Phi=%.1f_%d_%.2f",idef->iparams[i].pdihs.phiA,
 100                   idef->iparams[i].pdihs.mult,idef->iparams[i].pdihs.cpA);
 101           break;
 102         case F_IDIHS:
 103           sprintf(buf,"Xi=%.1f_%.2f",idef->iparams[i].harmonic.rA,
 104                   idef->iparams[i].harmonic.krA);
 105           break;
 106         case F_RBDIHS:
 107           sprintf(buf,"RB-A1=%.2f",idef->iparams[i].rbdihs.rbcA[1]);
 108           break;
 109         default:
 110           gmx_fatal(FARGS,"Unsupported function type '%s' selected",
 111                     interaction_function[ftype].longname);
 112         }
 113         grpnames[ind]=strdup(buf);
 114         ind++;
 115       }
 116     }
 117   }
 118 }
 119
 120 static void fill_ang(int nft,int *ft,int fac,
 121                      int nr[],int *index[],int ft_ind[],t_topology *top,
 122                      bool bNoH)
 123 {
 124   int     f,ftype,i,j,indg,nr_fac;
 125   bool    bUse;
 126   t_idef  *idef;
 127   t_atom  *atom;
 128   t_iatom *ia;
 129
 130
 131   idef = &top->idef;
 132   atom = top->atoms.atom;
 133
 134   for(f=0; f<nft; f++) {
 135     ftype = ft[f];
 136     ia = idef->il[ftype].iatoms;
 137     for(i=0; (i<idef->il[ftype].nr); ) {
 138       indg = ft_ind[ia[0]];
 139       if (indg == -1)
 140         gmx_incons("Routine fill_ang");
 141       bUse = TRUE;
 142       if (bNoH) {
 143         for(j=0; j<fac; j++) {
 144           if (atom[ia[1+j]].m < 1.5)
 145             bUse = FALSE;
 146         }
 147       }
 148       if (bUse) {
 149         if (nr[indg] % 1000 == 0) {
 150           srenew(index[indg],fac*(nr[indg]+1000));
 151         }
 152         nr_fac = fac*nr[indg];
 153         for(j=0; (j<fac); j++)
 154           index[indg][nr_fac+j] = ia[j+1];
 155         nr[indg]++;
 156       }
 157       ia += interaction_function[ftype].nratoms+1;
 158       i  += interaction_function[ftype].nratoms+1;
 159     }
 160   }
 161 }
 162
 163 static int *select_ftype(char *opt,int *nft,int *mult)
 164 {
 165   int *ft=NULL,ftype;
 166
 167   if (opt[0] == 'a') {
 168     *mult = 3;
 169     for(ftype=0; ftype<F_NRE; ftype++) {
 170       if (interaction_function[ftype].flags & IF_ATYPE ||
 171           ftype == F_TABANGLES) {
 172         (*nft)++;
 173         srenew(ft,*nft);
 174         ft[*nft-1] = ftype;
 175       }
 176     }
 177   } else {
 178     *mult = 4;
 179     *nft = 1;
 180     snew(ft,*nft);
 181     switch(opt[0]) {
 182     case 'd':
 183       ft[0] = F_PDIHS;
 184       break;
 185     case 'i':
 186       ft[0] = F_IDIHS;
 187       break;
 188     case 'r':
 189       ft[0] = F_RBDIHS;
 190       break;
 191     default:
 192       break;
 193     }
 194   }
 195
 196   return ft;
 197 }
 198
 199 int main(int argc,char *argv[])
 200 {
 201   static char *desc[] = {
 202     "mk_angndx makes an index file for calculation of",
 203     "angle distributions etc. It uses a run input file ([TT].tpx[tt]) for the",
 204     "definitions of the angles, dihedrals etc."
 205   };
 206   static char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
 207   static bool bH=TRUE;
 208   t_pargs pa[] = {
 209     { "-type", FALSE, etENUM, {opt},
 210       "Type of angle" },
 211     { "-hyd", FALSE, etBOOL, {&bH},
 212       "Include angles with atoms with mass < 1.5" }
 213   };
 214
 215   FILE       *out;
 216   t_topology *top;
 217   int        i,j,ntype;
 218   int        nft=0,*ft,mult=0;
 219   int        **index;
 220   int        *ft_ind;
 221   int        *nr;
 222   char       **grpnames;
 223   t_filenm fnm[] = {
 224     { efTPX, NULL, NULL, ffREAD  },
 225     { efNDX, NULL, "angle", ffWRITE }
 226   };
 227 #define NFILE asize(fnm)
 228
 229   CopyRight(stderr,argv[0]);
 230   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 231                     asize(desc),desc,0,NULL);
 232
 233
 234   ft = select_ftype(opt[0],&nft,&mult);
 235
 236   top = read_top(ftp2fn(efTPX,NFILE,fnm),NULL);
 237
 238   ntype = calc_ntype(nft,ft,&(top->idef));
 239   snew(grpnames,ntype);
 240   snew(ft_ind,top->idef.ntypes);
 241   fill_ft_ind(nft,ft,&top->idef,ft_ind,grpnames);
 242
 243   snew(nr,ntype);
 244   snew(index,ntype);
 245   fill_ang(nft,ft,mult,nr,index,ft_ind,top,!bH);
 246
 247   out=ftp2FILE(efNDX,NFILE,fnm,"w");
 248   for(i=0; (i<ntype); i++) {
 249     if (nr[i] > 0) {
 250       fprintf(out,"[ %s ]\n",grpnames[i]);
 251       for(j=0; (j<nr[i]*mult); j++) {
 252         fprintf(out," %5d",index[i][j]+1);
 253         if ((j % 12) == 11)
 254           fprintf(out,"\n");
 255       }
 256       fprintf(out,"\n");
 257     }
 258   }
 259   fclose(out);
 260
 261   thanx(stderr);
 262
 263   return 0;
 264 }