src/kernel/ter_db.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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_ter_db_c = "$Id$";
  30
  31 #include "sysstuff.h"
  32 #include "smalloc.h"
  33 #include "typedefs.h"
  34 #include "symtab.h"
  35 #include "futil.h"
  36 #include "resall.h"
  37 #include "h_db.h"
  38 #include "string2.h"
  39 #include "strdb.h"
  40 #include "fatal.h"
  41 #include "ter_db.h"
  42 #include "toputil.h"
  43
  44 /* use bonded types definitions in hackblock.h */
  45 #define ekwRepl ebtsNR+1
  46 #define ekwAdd  ebtsNR+2
  47 #define ekwDel  ebtsNR+3
  48 #define ekwNR   3
  49 char *kw_names[ekwNR] = {
  50   "replace", "add", "delete"
  51 };
  52
  53 int find_kw(char *keyw)
  54 {
  55   int i;
  56
  57   for(i=0; i<ebtsNR; i++)
  58     if (strcasecmp(btsNames[i],keyw) == 0)
  59       return i;
  60   for(i=0; i<ekwNR; i++)
  61     if (strcasecmp(kw_names[i],keyw) == 0)
  62       return ebtsNR + 1 + i;
  63
  64   return NOTSET;
  65 }
  66
  67 #define FATAL() fatal_error(0,"Reading Termini Database: not enough items on line\n%s",line)
  68
  69 static void read_atom(char *line, t_atom *a, t_atomtype *atype,
  70                       char nnew[], int *cgnr)
  71 {
  72   int    i, n;
  73   char   type[12];
  74   double m, q;
  75
  76   if ( (i=sscanf(line,"%s%s%lf%lf%n", nnew, type, &m, &q, &n)) != 4 )
  77     fatal_error(0,"Reading Termini Database: expected %d items of atom data in stead of %d on line\n%s", 4, i, line);
  78   a->m=m;
  79   a->q=q;
  80   a->type=at2type(type,atype);
  81   if ( sscanf(line+n,"%d", cgnr) != 1 )
  82     *cgnr = NOTSET;
  83 }
  84
  85 int read_ter_db(char *inf,t_hackblock **tbptr,t_atomtype *atype)
  86 {
  87   FILE       *in;
  88   char       header[STRLEN],buf[STRLEN],line[STRLEN];
  89  t_hackblock *tb;
  90   int        i,j,n,ni,kwnr,nb,maxnb,nh;
  91
  92   in=libopen(inf);
  93   if (debug)
  94     fprintf(debug,"Opened %s\n",inf);
  95
  96   tb=NULL;
  97   nb=-1;
  98   maxnb=0;
  99   kwnr=NOTSET;
 100   get_a_line(in,line,STRLEN);
 101   while (!feof(in)) {
 102     if (get_header(line,header)) {
 103       /* this is a new block, or a new keyword */
 104       kwnr=find_kw(header);
 105
 106       if (kwnr == NOTSET) {
 107         nb++;
 108         /* here starts a new block */
 109         if ( nb >= maxnb ) {
 110           maxnb+=100;
 111           srenew(tb,maxnb);
 112         }
 113         clear_t_hackblock(&tb[nb]);
 114         tb[nb].name=strdup(header);
 115       }
 116     } else {
 117       if (nb < 0)
 118         fatal_error(0,"reading termini database: "
 119                     "directive expected before line:\n%s\n"
 120                     "This might be a file in an old format.",line);
 121       /* this is not a header, so it must be data */
 122       if (kwnr >= ebtsNR) {
 123         /* this is a hack: add/rename/delete atoms */
 124         /* make space for hacks */
 125         if (tb[nb].nhack >= tb[nb].maxhack) {
 126           tb[nb].maxhack+=10;
 127           srenew(tb[nb].hack, tb[nb].maxhack);
 128         }
 129         nh=tb[nb].nhack;
 130         clear_t_hack(&(tb[nb].hack[nh]));
 131         for(i=0; i<4; i++)
 132           tb[nb].hack[nh].a[i]=NULL;
 133         tb[nb].nhack++;
 134
 135         /* get data */
 136         n=0;
 137         if ( kwnr==ekwRepl || kwnr==ekwDel ) {
 138           if (sscanf(line, "%s%n", buf, &n) != 1)
 139             fatal_error(0,"Reading Termini Database: "
 140                         "expected atom name on line\n%s",line);
 141           tb[nb].hack[nh].oname = strdup(buf);
 142           /* we only replace or delete one atom at a time */
 143           tb[nb].hack[nh].nr = 1;
 144         } else if ( kwnr==ekwAdd ) {
 145           read_ab(line, inf, &(tb[nb].hack[nh]));
 146           get_a_line(in, line, STRLEN);
 147         } else
 148           fatal_error(0,"unimplemented keyword number %d (%s:%d)",
 149                       kwnr,__FILE__,__LINE__);
 150         if ( kwnr==ekwRepl || kwnr==ekwAdd ) {
 151           snew(tb[nb].hack[nh].atom, 1);
 152           read_atom(line+n, tb[nb].hack[nh].atom, atype,
 153                     buf, &tb[nb].hack[nh].cgnr);
 154           tb[nb].hack[nh].nname = strdup(buf);
 155         }
 156       } else if (kwnr >= 0 && kwnr < ebtsNR) {
 157         /* this is bonded data: bonds, angles, dihedrals or impropers */
 158         srenew(tb[nb].rb[kwnr].b,tb[nb].rb[kwnr].nb+1);
 159         n=0;
 160         for(j=0; j<btsNiatoms[kwnr]; j++) {
 161           if ( sscanf(line+n, "%s%n", buf, &ni) == 1 )
 162             tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].a[j] = strdup(buf);
 163           else
 164             fatal_error(0,"Reading Termini Database: expected %d atom names (found %d) on line\n%s", btsNiatoms[kwnr], j-1, line);
 165           n+=ni;
 166         }
 167         for(   ; j<MAXATOMLIST; j++)
 168           tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].a[j] = NULL;
 169         strcpy(buf, "");
 170         sscanf(line+n, "%s", buf);
 171         tb[nb].rb[kwnr].b[tb[nb].rb[kwnr].nb].s = strdup(buf);
 172         tb[nb].rb[kwnr].nb++;
 173       } else
 174         fatal_error(0,"Reading Termini Database: Expecting a header at line\n"
 175                     "%s",line);
 176     }
 177     get_a_line(in,line,STRLEN);
 178   }
 179   nb++;
 180   srenew(tb,nb);
 181
 182   fclose(in);
 183
 184   if (debug) print_ter_db(debug, nb, tb, atype);
 185
 186   *tbptr=tb;
 187   return nb;
 188 }
 189
 190 t_hackblock *choose_ter(int nb,t_hackblock tb[],char *title)
 191 {
 192   int i,sel;
 193
 194   if (nb == 1)
 195     return &(tb[0]);
 196
 197   printf("%s\n",title);
 198   for(i=0; (i<nb); i++)
 199     printf("%2d: %s\n",i,tb[i].name);
 200   do {
 201     fscanf(stdin,"%d",&sel);
 202   } while ((sel < 0) || (sel >= nb));
 203
 204   return &(tb[sel]);
 205 }
 206
 207 static void print_atom(FILE *out,t_atom *a,t_atomtype *atype,char *newnm)
 208 {
 209   fprintf(out,"%s\t%s\t%g\t%g\n",
 210           newnm,type2nm(a->type,atype),a->m,a->q);
 211 }
 212
 213 void print_ter_db(FILE *out,int nb,t_hackblock tb[],t_atomtype *atype)
 214 {
 215   int i,j,k,bt,nrepl,nadd,ndel;
 216
 217   for(i=0; (i<nb); i++) {
 218     fprintf(out,"[ %s ]\n",tb[i].name);
 219
 220     /* first count: */
 221     nrepl=0;
 222     nadd=0;
 223     ndel=0;
 224     for(j=0; j<tb[i].nhack; j++)
 225       if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname!=NULL )
 226         nrepl++;
 227       else if ( tb[i].hack[j].oname==NULL && tb[i].hack[j].nname!=NULL )
 228         nadd++;
 229       else if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname==NULL )
 230         ndel++;
 231       else if ( tb[i].hack[j].oname==NULL && tb[i].hack[j].nname==NULL )
 232         fatal_error(0,"invalid hack (%s) in termini database",tb[i].name);
 233     if (nrepl) {
 234       fprintf(out,"[ %s ]\n",kw_names[ekwRepl-ebtsNR-1]);
 235       for(j=0; j<tb[i].nhack; j++)
 236         if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname!=NULL ) {
 237           fprintf(out,"%s\t",tb[i].hack[j].oname);
 238           print_atom(out,tb[i].hack[j].atom,atype,tb[i].hack[j].nname);
 239         }
 240     }
 241     if (nadd) {
 242       fprintf(out,"[ %s ]\n",kw_names[ekwAdd-ebtsNR-1]);
 243       for(j=0; j<tb[i].nhack; j++)
 244         if ( tb[i].hack[j].oname==NULL && tb[i].hack[j].nname!=NULL ) {
 245           print_ab(out,&(tb[i].hack[j]));
 246           fprintf(out,"\t");
 247           print_atom(out,tb[i].hack[j].atom,atype,tb[i].hack[j].nname);
 248         }
 249     }
 250     if (ndel) {
 251       fprintf(out,"[ %s ]\n",kw_names[ekwDel-ebtsNR-1]);
 252       for(j=0; j<tb[i].nhack; j++)
 253         if ( tb[i].hack[j].oname!=NULL && tb[i].hack[j].nname==NULL )
 254           fprintf(out,"%s\n",tb[i].hack[j].oname);
 255     }
 256     for(bt=0; bt<ebtsNR; bt++)
 257       if (tb[i].rb[bt].nb) {
 258         fprintf(out,"[ %s ]\n", btsNames[bt]);
 259         for(j=0; j<tb[i].rb[bt].nb; j++) {
 260           for(k=0; k<btsNiatoms[bt]; k++)
 261             fprintf(out,"%s%s",k?"\t":"",tb[i].rb[bt].b[j].a[k]);
 262           if ( tb[i].rb[bt].b[j].s )
 263             fprintf(out,"\t%s",tb[i].rb[bt].b[j].s);
 264           fprintf(out,"\n");
 265         }
 266       }
 267     fprintf(out,"\n");
 268   }
 269 }
 270
 271 #ifdef DBTDB
 272 int main(int argc,char *argv[])
 273 {
 274   t_symtab   tab;
 275   t_atomtype *atype;
 276   t_hackblock *tb,*seltb;
 277   int nb;
 278
 279   open_symtab(&tab);
 280   atype=read_atype("ffgmx",&tab);
 281   nb=read_ter_db("ffgmx-c.tdb",&tb,atype);
 282   seltb=choose_ter(nb,tb,"What do you want ?");
 283   print_ter_db(stdout,1,seltb,atype);
 284
 285   return 0;
 286 }
 287 #endif