src/tools/mk_angndx.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 <math.h>
  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                      gmx_bool bNoH,real hq)
 123 {
 124   int     f,ftype,i,j,indg,nr_fac;
 125   gmx_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 (hq) {
 149         for(j=0; j<fac; j++) {
 150           if (atom[ia[1+j]].m < 1.5 && fabs(atom[ia[1+j]].q) < hq)
 151             bUse = FALSE;
 152         }
 153       }
 154       if (bUse) {
 155         if (nr[indg] % 1000 == 0) {
 156           srenew(index[indg],fac*(nr[indg]+1000));
 157         }
 158         nr_fac = fac*nr[indg];
 159         for(j=0; (j<fac); j++)
 160           index[indg][nr_fac+j] = ia[j+1];
 161         nr[indg]++;
 162       }
 163       ia += interaction_function[ftype].nratoms+1;
 164       i  += interaction_function[ftype].nratoms+1;
 165     }
 166   }
 167 }
 168
 169 static int *select_ftype(const char *opt,int *nft,int *mult)
 170 {
 171   int *ft=NULL,ftype;
 172
 173   if (opt[0] == 'a') {
 174     *mult = 3;
 175     for(ftype=0; ftype<F_NRE; ftype++) {
 176       if ((interaction_function[ftype].flags & IF_ATYPE) ||
 177           ftype == F_TABANGLES) {
 178         (*nft)++;
 179         srenew(ft,*nft);
 180         ft[*nft-1] = ftype;
 181       }
 182     }
 183   } else {
 184     *mult = 4;
 185     *nft = 1;
 186     snew(ft,*nft);
 187     switch(opt[0]) {
 188     case 'd':
 189       ft[0] = F_PDIHS;
 190       break;
 191     case 'i':
 192       ft[0] = F_IDIHS;
 193       break;
 194     case 'r':
 195       ft[0] = F_RBDIHS;
 196       break;
 197     default:
 198       break;
 199     }
 200   }
 201
 202   return ft;
 203 }
 204
 205 int main(int argc,char *argv[])
 206 {
 207   static const char *desc[] = {
 208     "[TT]mk_angndx[tt] makes an index file for calculation of",
 209     "angle distributions etc. It uses a run input file ([TT].tpx[tt]) for the",
 210     "definitions of the angles, dihedrals etc."
 211   };
 212   static const char *opt[] = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
 213   static gmx_bool bH=TRUE;
 214   static real hq=-1;
 215   t_pargs pa[] = {
 216     { "-type", FALSE, etENUM, {opt},
 217       "Type of angle" },
 218     { "-hyd", FALSE, etBOOL, {&bH},
 219       "Include angles with atoms with mass < 1.5" },
 220     { "-hq", FALSE, etREAL, {&hq},
 221       "Ignore angles with atoms with mass < 1.5 and |q| < hq" }
 222   };
 223
 224   output_env_t oenv;
 225   FILE       *out;
 226   t_topology *top;
 227   int        i,j,ntype;
 228   int        nft=0,*ft,mult=0;
 229   int        **index;
 230   int        *ft_ind;
 231   int        *nr;
 232   char       **grpnames;
 233   t_filenm fnm[] = {
 234     { efTPX, NULL, NULL, ffREAD  },
 235     { efNDX, NULL, "angle", ffWRITE }
 236   };
 237 #define NFILE asize(fnm)
 238
 239   CopyRight(stderr,argv[0]);
 240   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
 241                     asize(desc),desc,0,NULL,&oenv);
 242
 243
 244   ft = select_ftype(opt[0],&nft,&mult);
 245
 246   top = read_top(ftp2fn(efTPX,NFILE,fnm),NULL);
 247
 248   ntype = calc_ntype(nft,ft,&(top->idef));
 249   snew(grpnames,ntype);
 250   snew(ft_ind,top->idef.ntypes);
 251   fill_ft_ind(nft,ft,&top->idef,ft_ind,grpnames);
 252
 253   snew(nr,ntype);
 254   snew(index,ntype);
 255   fill_ang(nft,ft,mult,nr,index,ft_ind,top,!bH,hq);
 256
 257   out=ftp2FILE(efNDX,NFILE,fnm,"w");
 258   for(i=0; (i<ntype); i++) {
 259     if (nr[i] > 0) {
 260       fprintf(out,"[ %s ]\n",grpnames[i]);
 261       for(j=0; (j<nr[i]*mult); j++) {
 262         fprintf(out," %5d",index[i][j]+1);
 263         if ((j % 12) == 11)
 264           fprintf(out,"\n");
 265       }
 266       fprintf(out,"\n");
 267     }
 268   }
 269   ffclose(out);
 270
 271   thanx(stderr);
 272
 273   return 0;
 274 }