src/kernel/resall.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_resall_c = "$Id$";
  30
  31 #include "sysstuff.h"
  32 #include "assert.h"
  33 #include "string2.h"
  34 #include "strdb.h"
  35 #include "futil.h"
  36 #include "smalloc.h"
  37 #include "fatal.h"
  38 #include "symtab.h"
  39 #include "macros.h"
  40 #include "resall.h"
  41 #include "pgutil.h"
  42
  43 t_atomtype *read_atype(char *adb,t_symtab *tab)
  44 {
  45 #define MAXAT 1000
  46   FILE       *in;
  47   char       aadb[STRLEN];
  48   char       buf[STRLEN],name[STRLEN];
  49   double     m;
  50   int        nratt;
  51   t_atomtype *at;
  52
  53   sprintf(aadb,"%s.atp",adb);
  54   in=libopen(aadb);
  55
  56   snew(at,1);
  57   snew(at->atom,MAXAT);
  58   snew(at->atomname,MAXAT);
  59
  60   for(nratt=0; ; nratt++) {
  61     if (nratt >= MAXAT)
  62       fatal_error(0,"nratt >= MAXAT(%d). Increase the latter",MAXAT);
  63     if (feof(in))
  64       break;
  65     if (fgets2(buf,STRLEN,in) == NULL)
  66       break;
  67     if (sscanf(buf,"%s%lf",name,&m) != 2)
  68       break;
  69     set_at(&(at->atom[nratt]),m,0.0,nratt,0);
  70     at->atomname[nratt]=put_symtab(tab,name);
  71     fprintf(stderr,"\rAtomtype %d",nratt+1);
  72   }
  73   fclose(in);
  74   fprintf(stderr,"\n");
  75   at->nr=nratt;
  76
  77   return at;
  78 }
  79
  80 static void print_resatoms(FILE *out,t_atomtype *atype,t_restp *rtp)
  81 {
  82   int j,tp;
  83
  84   /* fprintf(out,"%5s\n",rtp->resname);
  85      fprintf(out,"%5d\n",rtp->natom); */
  86   fprintf(out,"[ %s ]\n",rtp->resname);
  87   fprintf(out," [ atoms ]\n");
  88
  89   for(j=0; (j<rtp->natom); j++) {
  90     tp=rtp->atom[j].type;
  91     assert (tp >= 0);
  92     assert (tp < atype->nr);
  93     fprintf(out,"%6s%6s%8.3f%6d\n",
  94             *(rtp->atomname[j]),*(atype->atomname[tp]),
  95             rtp->atom[j].q,rtp->cgnr[j]);
  96   }
  97 }
  98
  99 static bool read_atoms(FILE *in,char *line,
 100                        t_restp *r0,t_symtab *tab,t_atomtype *atype)
 101 {
 102   int    i,j,cg,maxentries;
 103   char   buf[256],buf1[256];
 104   double q;
 105
 106   /* Read Atoms */
 107   maxentries=0;
 108   r0->atom=     NULL;
 109   r0->atomname= NULL;
 110   r0->cgnr=     NULL;
 111   i=0;
 112   while (get_a_line(in,line,STRLEN) && (strchr(line,'[')==NULL)) {
 113     if (sscanf(line,"%s%s%lf%d",buf,buf1,&q,&cg) != 4)
 114       return FALSE;
 115     if (i>=maxentries) {
 116       maxentries+=100;
 117       srenew(r0->atom,     maxentries);
 118       srenew(r0->atomname, maxentries);
 119       srenew(r0->cgnr,     maxentries);
 120     }
 121     r0->atomname[i]=put_symtab(tab,buf);
 122     r0->atom[i].q=q;
 123     r0->cgnr[i]=cg;
 124     for(j=0; (j<atype->nr); j++)
 125       if (strcasecmp(buf1,*(atype->atomname[j])) == 0)
 126         break;
 127     if (j == atype->nr)
 128       fatal_error(0,"Atom type %s (residue %s) not found in atomtype "
 129                   "database",buf1,r0->resname);
 130     r0->atom[i].type=j;
 131     r0->atom[i].m=atype->atom[j].m;
 132     i++;
 133   }
 134   r0->natom=i;
 135   srenew(r0->atom,i);
 136   srenew(r0->atomname,i);
 137   srenew(r0->cgnr,i);
 138
 139   return TRUE;
 140 }
 141
 142 bool read_bondeds(int bt, FILE *in, char *line, t_restp *rtp)
 143 {
 144   char str[STRLEN];
 145   int  j,n,ni,maxrb;
 146
 147   maxrb = rtp->rb[bt].nb;
 148   while (get_a_line(in,line,STRLEN) && (strchr(line,'[')==NULL)) {
 149     if ( rtp->rb[bt].nb >= maxrb ) {
 150       maxrb+=100;
 151       srenew(rtp->rb[bt].b,maxrb);
 152     }
 153     n=0;
 154     for(j=0; j < btsNiatoms[bt]; j++) {
 155       if ( sscanf(line+n,"%s%n",str,&ni)==1 )
 156         rtp->rb[bt].b[rtp->rb[bt].nb].a[j]=strdup(str);
 157       else
 158         return FALSE;
 159       n+=ni;
 160     }
 161     for(  ; j < MAXATOMLIST; j++)
 162       rtp->rb[bt].b[rtp->rb[bt].nb].a[j]=NULL;
 163     strcpy(str, "");
 164     sscanf(line+n,"%s",str);
 165     rtp->rb[bt].b[rtp->rb[bt].nb].s=strdup(str);
 166     rtp->rb[bt].nb++;
 167   }
 168   /* give back unused memory */
 169   srenew(rtp->rb[bt].b,rtp->rb[bt].nb);
 170
 171   return TRUE;
 172 }
 173
 174 static void print_resbondeds(FILE *out, int bt, t_restp *rtp)
 175 {
 176   int i,j;
 177
 178   if (rtp->rb[bt].nb) {
 179     fprintf(out," [ %s ]\n",btsNames[bt]);
 180
 181     for(i=0; i < rtp->rb[bt].nb; i++) {
 182       for(j=0; j<btsNiatoms[bt]; j++)
 183         fprintf(out,"%6s ",rtp->rb[bt].b[i].a[j]);
 184       if (rtp->rb[bt].b[i].s[0])
 185         fprintf(out,"    %s",rtp->rb[bt].b[i].s);
 186       fprintf(out,"\n");
 187     }
 188   }
 189 }
 190
 191 static void check_rtp(int nrtp,t_restp rtp[],char *libfn)
 192 {
 193   int i;
 194
 195   /* check for double entries, assuming list is already sorted */
 196   for(i=1; (i<nrtp); i++) {
 197     if (strcasecmp(rtp[i-1].resname,rtp[i].resname) == 0)
 198       fprintf(stderr,"WARNING double entry %s in file %s\n",
 199               rtp[i].resname,libfn);
 200   }
 201 }
 202
 203 static int comprtp(const void *a,const void *b)
 204 {
 205   t_restp *ra,*rb;
 206
 207   ra=(t_restp *)a;
 208   rb=(t_restp *)b;
 209
 210   return strcasecmp(ra->resname,rb->resname);
 211 }
 212
 213 int get_bt(char* header)
 214 {
 215   int i;
 216
 217   for(i=0; i<ebtsNR; i++)
 218     if ( strcasecmp(btsNames[i],header)==0 )
 219       return i;
 220   return NOTSET;
 221 }
 222
 223 void clear_t_restp(t_restp *rrtp)
 224 {
 225   memset((void *)rrtp, 0, sizeof(t_restp));
 226 }
 227
 228 int read_resall(char *ff, int bts[], t_restp **rtp,
 229                 t_atomtype *atype, t_symtab *tab)
 230 {
 231   FILE      *in;
 232   char      rrdb[STRLEN],line[STRLEN],header[STRLEN];
 233   int       i,nrtp,maxrtp,bt;
 234   t_restp   *rrtp;
 235   bool      bNextResidue,bError;
 236
 237   sprintf(rrdb,"%s.rtp",ff);
 238   in=libopen(rrdb);
 239   rrtp=NULL;
 240   if (debug) {
 241     fprintf(debug,"%9s %5s", "Residue", "atoms");
 242     for(i=0; i<ebtsNR; i++)
 243       fprintf(debug," %10s",btsNames[i]);
 244     fprintf(debug,"\n");
 245   }
 246
 247   /* these bonded parameters will overwritten be when  *
 248    * there is a [ bondedtypes ] entry in the .rtp file */
 249   bts[0] = 1; /* normal bonds     */
 250   bts[1] = 1; /* normal angles    */
 251   bts[2] = 1; /* normal dihedrals */
 252   bts[3] = 2; /* normal impropers */
 253
 254   nrtp=0;
 255   maxrtp=0;
 256   get_a_line(in,line,STRLEN);
 257   if (!get_header(line,header))
 258     fatal_error(0,"in .rtp file at line:\n%s\n",line);
 259   if (strncasecmp("bondedtypes",header,5)==0) {
 260     get_a_line(in,line,STRLEN);
 261     if (sscanf(line,"%d %d %d %d",&bts[0],&bts[1],&bts[2],&bts[3]) != ebtsNR)
 262       fatal_error(0,"need %d parameters in .rtp file at line:\n%s\n",
 263                   ebtsNR,line);
 264     get_a_line(in,line,STRLEN);
 265   } else {
 266     fprintf(stderr,
 267             "Reading .rtp file without '[ bondedtypes ]' directive,\n"
 268             "Will proceed as if the entry\n"
 269             "\n"
 270             "\n[ bondedtypes ]"
 271             "\n; bonds  angles  dihedrals  impropers"
 272             "\n   %3d     %3d        %3d        %3d"
 273             "\n"
 274             "was present at the beginning of %s",
 275             bts[0],bts[1],bts[2],bts[3],rrdb);
 276   }
 277   nrtp=0;
 278   while (!feof(in)) {
 279     if (nrtp >= maxrtp) {
 280       maxrtp+=100;
 281       srenew(rrtp,maxrtp);
 282     }
 283     clear_t_restp(&rrtp[nrtp]);
 284     if (!get_header(line,header))
 285       fatal_error(0,"in .rtp file at line:\n%s\n",line);
 286     rrtp[nrtp].resname=strdup(header);
 287
 288     get_a_line(in,line,STRLEN);
 289     bError=FALSE;
 290     bNextResidue=FALSE;
 291     while (!bNextResidue) {
 292       if (!get_header(line,header))
 293         if (feof(in))
 294           bNextResidue=TRUE;
 295         else
 296           bError=TRUE;
 297       else {
 298         bt=get_bt(header);
 299         if ( bt!=NOTSET )
 300           bError=!read_bondeds(bt,in,line,&rrtp[nrtp]);
 301         else
 302           if (strncasecmp("atoms",header,5)==0)
 303             bError=!read_atoms(in,line,&(rrtp[nrtp]),tab,atype);
 304           else {
 305             if (!feof(in) && !get_header(line,header))
 306               bError=TRUE;
 307             bNextResidue=TRUE;
 308           }
 309       }
 310       if (bError)
 311         fatal_error(0,"in .rtp file in residue %s at line:\n%s\n",
 312                     rrtp[nrtp].resname,line);
 313     }
 314     if (rrtp[nrtp].natom == 0)
 315       fatal_error(0,"No atoms found in .rtp file in residue %s\n",
 316                   rrtp[nrtp].resname);
 317     if (debug) {
 318       fprintf(debug,"%3d %5s %5d",
 319               nrtp+1,rrtp[nrtp].resname,rrtp[nrtp].natom);
 320       for(i=0; i<ebtsNR; i++)
 321         fprintf(debug," %10d",rrtp[nrtp].rb[i].nb);
 322       fprintf(debug,"\n");
 323     }
 324     nrtp++;
 325     fprintf(stderr,"\rResidue %d",nrtp);
 326   }
 327   fclose(in);
 328   /* give back unused memory */
 329   srenew(rrtp,nrtp);
 330
 331   fprintf(stderr,"\nSorting it all out...\n");
 332   qsort(rrtp,nrtp,(size_t)sizeof(rrtp[0]),comprtp);
 333
 334   check_rtp(nrtp,rrtp,rrdb);
 335
 336   *rtp  = rrtp;
 337
 338   return nrtp;
 339 }
 340
 341 void print_resall(FILE *out, int bts[], int nrtp, t_restp rtp[],
 342                   t_atomtype *atype)
 343 {
 344   int i,bt;
 345
 346   /* print all the ebtsNR type numbers */
 347   fprintf(out,"[ bondedtypes ]\n");
 348   fprintf(out,"; bonds  angles  dihedrals  impropers\n");
 349   fprintf(out," %5d  %6d  %9d  %9d\n\n",bts[0],bts[1],bts[2],bts[3]);
 350
 351   for(i=0; i<nrtp; i++) {
 352     if (rtp[i].natom > 0) {
 353       print_resatoms(out,atype,&rtp[i]);
 354       for(bt=0; bt<ebtsNR; bt++)
 355         print_resbondeds(out,bt,&rtp[i]);
 356     }
 357   }
 358 }
 359
 360 /************************************************************
 361  *
 362  *                  SEARCH   ROUTINES
 363  *
 364  ***********************************************************/
 365 int neq_str(char *a1,char *a2)
 366 {
 367   int j,l;
 368
 369   l=min((int)strlen(a1),(int)strlen(a2));
 370   j=0;
 371   while ( (j<l) && (toupper(a1[j]) == toupper(a2[j])) )
 372     j++;
 373
 374   return j;
 375 }
 376
 377 t_restp *search_rtp(char *key,int nrtp,t_restp rtp[])
 378 {
 379   int i,n,best,besti;
 380
 381   besti=-1;
 382   best=1;
 383   for(i=0; (i<nrtp); i++) {
 384     n=neq_str(key,rtp[i].resname);
 385     if (n > best) {
 386       besti=i;
 387       best=n;
 388     }
 389   }
 390   if (besti == -1)
 391     fatal_error(0,"Residue '%s' not found in residue topology database\n",key);
 392   if (strlen(rtp[besti].resname) != strlen(key))
 393     fprintf(stderr,"Warning: '%s' not found in residue topology database, "
 394             "trying to use '%s'\n", key, rtp[besti].resname);
 395
 396   return &rtp[besti];
 397 }