src/gromacs/gmxpreprocess/resall.c

   1 /*  -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   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  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "sysstuff.h"
  41 #include <ctype.h>
  42 #include "string2.h"
  43 #include "strdb.h"
  44 #include "futil.h"
  45 #include "smalloc.h"
  46 #include "gmx_fatal.h"
  47 #include "symtab.h"
  48 #include "macros.h"
  49 #include "resall.h"
  50 #include "pgutil.h"
  51 #include "fflibutil.h"
  52
  53 gpp_atomtype_t read_atype(const char *ffdir,t_symtab *tab)
  54 {
  55     int        nfile,f;
  56     char       **file;
  57     FILE       *in;
  58     char       buf[STRLEN],name[STRLEN];
  59     double     m;
  60     int        nratt=0;
  61     gpp_atomtype_t at;
  62     t_atom     *a;
  63     t_param    *nb;
  64
  65     nfile = fflib_search_file_end(ffdir,".atp",TRUE,&file);
  66     at = init_atomtype();
  67     snew(a,1);
  68     snew(nb,1);
  69
  70     for(f=0; f<nfile; f++)
  71     {
  72         in = fflib_open(file[f]);
  73         while(!feof(in))
  74         {
  75             /* Skip blank or comment-only lines */
  76             do
  77             {
  78                 fgets2(buf,STRLEN,in);
  79                 if (NULL != buf)
  80                 {
  81                     strip_comment(buf);
  82                     trim(buf);
  83                 }
  84             }
  85             while (!feof(in) && NULL!=buf && strlen(buf)==0);
  86
  87             if ((buf != NULL) && (sscanf(buf,"%s%lf",name,&m) == 2))
  88             {
  89                 a->m = m;
  90                 add_atomtype(at,tab,a,name,nb, 0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 );
  91                 fprintf(stderr,"\rAtomtype %d",nratt+1);
  92             }
  93         }
  94         ffclose(in);
  95         sfree(file[f]);
  96     }
  97     fprintf(stderr,"\n");
  98     sfree(file);
  99
 100     return at;
 101 }
 102
 103 static void print_resatoms(FILE *out,gpp_atomtype_t atype,t_restp *rtp)
 104 {
 105   int j,tp;
 106   char *tpnm;
 107
 108   /* fprintf(out,"%5s\n",rtp->resname);
 109      fprintf(out,"%5d\n",rtp->natom); */
 110   fprintf(out,"[ %s ]\n",rtp->resname);
 111   fprintf(out," [ atoms ]\n");
 112
 113   for(j=0; (j<rtp->natom); j++) {
 114     tp = rtp->atom[j].type;
 115     tpnm = get_atomtype_name(tp,atype);
 116     if (tpnm == NULL)
 117       gmx_fatal(FARGS,"Incorrect atomtype (%d)",tp);
 118     fprintf(out,"%6s  %6s  %8.3f  %6d\n",
 119             *(rtp->atomname[j]),tpnm,rtp->atom[j].q,rtp->cgnr[j]);
 120   }
 121 }
 122
 123 static gmx_bool read_atoms(FILE *in,char *line,
 124                        t_restp *r0,t_symtab *tab,gpp_atomtype_t atype)
 125 {
 126   int    i,j,cg,maxentries;
 127   char   buf[256],buf1[256];
 128   double q;
 129
 130   /* Read Atoms */
 131   maxentries=0;
 132   r0->atom=     NULL;
 133   r0->atomname= NULL;
 134   r0->cgnr=     NULL;
 135   i=0;
 136   while (get_a_line(in,line,STRLEN) && (strchr(line,'[')==NULL)) {
 137     if (sscanf(line,"%s%s%lf%d",buf,buf1,&q,&cg) != 4)
 138       return FALSE;
 139     if (i>=maxentries) {
 140       maxentries+=100;
 141       srenew(r0->atom,     maxentries);
 142       srenew(r0->atomname, maxentries);
 143       srenew(r0->cgnr,     maxentries);
 144     }
 145     r0->atomname[i] = put_symtab(tab,buf);
 146     r0->atom[i].q=q;
 147     r0->cgnr[i]=cg;
 148     j = get_atomtype_type(buf1,atype);
 149     if (j == NOTSET)
 150       gmx_fatal(FARGS,"Atom type %s (residue %s) not found in atomtype "
 151                   "database",buf1,r0->resname);
 152     r0->atom[i].type=j;
 153     r0->atom[i].m=get_atomtype_massA(j,atype);
 154     i++;
 155   }
 156   r0->natom=i;
 157   srenew(r0->atom,i);
 158   srenew(r0->atomname,i);
 159   srenew(r0->cgnr,i);
 160
 161   return TRUE;
 162 }
 163
 164 gmx_bool read_bondeds(int bt, FILE *in, char *line, t_restp *rtp)
 165 {
 166   char str[STRLEN];
 167   int  j,n,ni,maxrb;
 168
 169   maxrb = rtp->rb[bt].nb;
 170   while (get_a_line(in,line,STRLEN) && (strchr(line,'[')==NULL)) {
 171     if ( rtp->rb[bt].nb >= maxrb ) {
 172       maxrb+=100;
 173       srenew(rtp->rb[bt].b,maxrb);
 174     }
 175     n=0;
 176     for(j=0; j < btsNiatoms[bt]; j++) {
 177       if ( sscanf(line+n,"%s%n",str,&ni)==1 )
 178         rtp->rb[bt].b[rtp->rb[bt].nb].a[j]=strdup(str);
 179       else
 180         return FALSE;
 181       n+=ni;
 182     }
 183     for(  ; j < MAXATOMLIST; j++)
 184       rtp->rb[bt].b[rtp->rb[bt].nb].a[j]=NULL;
 185     while (isspace(line[n]))
 186       n++;
 187     rtrim(line+n);
 188     rtp->rb[bt].b[rtp->rb[bt].nb].s=strdup(line+n);
 189     rtp->rb[bt].nb++;
 190   }
 191   /* give back unused memory */
 192   srenew(rtp->rb[bt].b,rtp->rb[bt].nb);
 193
 194   return TRUE;
 195 }
 196
 197 static void print_resbondeds(FILE *out, int bt, t_restp *rtp)
 198 {
 199   int i,j;
 200
 201   if (rtp->rb[bt].nb) {
 202     fprintf(out," [ %s ]\n",btsNames[bt]);
 203
 204     for(i=0; i < rtp->rb[bt].nb; i++) {
 205       for(j=0; j<btsNiatoms[bt]; j++)
 206         fprintf(out,"%6s ",rtp->rb[bt].b[i].a[j]);
 207       if (rtp->rb[bt].b[i].s[0])
 208         fprintf(out,"    %s",rtp->rb[bt].b[i].s);
 209       fprintf(out,"\n");
 210     }
 211   }
 212 }
 213
 214 static void check_rtp(int nrtp,t_restp rtp[],char *libfn)
 215 {
 216   int i;
 217
 218   /* check for double entries, assuming list is already sorted */
 219   for(i=1; (i<nrtp); i++) {
 220     if (gmx_strcasecmp(rtp[i-1].resname,rtp[i].resname) == 0)
 221       fprintf(stderr,"WARNING double entry %s in file %s\n",
 222               rtp[i].resname,libfn);
 223   }
 224 }
 225
 226 static int comprtp(const void *a,const void *b)
 227 {
 228   t_restp *ra,*rb;
 229
 230   ra=(t_restp *)a;
 231   rb=(t_restp *)b;
 232
 233   return gmx_strcasecmp(ra->resname,rb->resname);
 234 }
 235
 236 int get_bt(char* header)
 237 {
 238   int i;
 239
 240   for(i=0; i<ebtsNR; i++)
 241     if ( gmx_strcasecmp(btsNames[i],header)==0 )
 242       return i;
 243   return NOTSET;
 244 }
 245
 246 void clear_t_restp(t_restp *rrtp)
 247 {
 248   memset((void *)rrtp, 0, sizeof(t_restp));
 249 }
 250
 251 /* print all the ebtsNR type numbers */
 252 void print_resall_header(FILE *out, t_restp rtp[])
 253 {
 254   fprintf(out,"[ bondedtypes ]\n");
 255   fprintf(out,"; bonds  angles  dihedrals  impropers all_dihedrals nr_exclusions  HH14  remove_dih\n");
 256   fprintf(out," %5d  %6d  %9d  %9d  %14d  %14d %14d %14d\n\n",
 257           rtp[0].rb[0].type,
 258           rtp[0].rb[1].type,
 259           rtp[0].rb[2].type,
 260           rtp[0].rb[3].type,
 261           rtp[0].bKeepAllGeneratedDihedrals,
 262           rtp[0].nrexcl,
 263           rtp[0].bGenerateHH14Interactions,
 264           rtp[0].bRemoveDihedralIfWithImproper);
 265 }
 266
 267 void print_resall(FILE *out, int nrtp, t_restp rtp[],
 268                   gpp_atomtype_t atype)
 269 {
 270   int i,bt;
 271
 272   if (nrtp == 0) {
 273     return;
 274   }
 275
 276   print_resall_header(out, rtp);
 277
 278   for(i=0; i<nrtp; i++) {
 279     if (rtp[i].natom > 0) {
 280       print_resatoms(out,atype,&rtp[i]);
 281       for(bt=0; bt<ebtsNR; bt++)
 282         print_resbondeds(out,bt,&rtp[i]);
 283     }
 284   }
 285 }
 286
 287 void read_resall(char *rrdb, int *nrtpptr, t_restp **rtp,
 288                  gpp_atomtype_t atype, t_symtab *tab,
 289                  gmx_bool bAllowOverrideRTP)
 290 {
 291   FILE      *in;
 292   char      filebase[STRLEN],*ptr,line[STRLEN],header[STRLEN];
 293   int       i,nrtp,maxrtp,bt,nparam;
 294   int       dum1,dum2,dum3;
 295   t_restp   *rrtp, *header_settings;
 296   gmx_bool      bNextResidue,bError;
 297   int       firstrtp;
 298
 299   fflib_filename_base(rrdb,filebase,STRLEN);
 300
 301   in = fflib_open(rrdb);
 302
 303   if (debug) {
 304     fprintf(debug,"%9s %5s", "Residue", "atoms");
 305     for(i=0; i<ebtsNR; i++)
 306       fprintf(debug," %10s",btsNames[i]);
 307     fprintf(debug,"\n");
 308   }
 309   snew(header_settings, 1);
 310
 311   /* these bonded parameters will overwritten be when  *
 312    * there is a [ bondedtypes ] entry in the .rtp file */
 313   header_settings->rb[ebtsBONDS].type  = 1; /* normal bonds     */
 314   header_settings->rb[ebtsANGLES].type = 1; /* normal angles    */
 315   header_settings->rb[ebtsPDIHS].type  = 1; /* normal dihedrals */
 316   header_settings->rb[ebtsIDIHS].type  = 2; /* normal impropers */
 317   header_settings->rb[ebtsEXCLS].type  = 1; /* normal exclusions */
 318   header_settings->rb[ebtsCMAP].type   = 1; /* normal cmap torsions */
 319
 320   header_settings->bKeepAllGeneratedDihedrals    = FALSE;
 321   header_settings->nrexcl                        = 3;
 322   header_settings->bGenerateHH14Interactions     = TRUE;
 323   header_settings->bRemoveDihedralIfWithImproper = TRUE;
 324
 325   /* Column 5 & 6 aren't really bonded types, but we include
 326    * them here to avoid introducing a new section:
 327    * Column 5 : This controls the generation of dihedrals from the bonding.
 328    *            All possible dihedrals are generated automatically. A value of
 329    *            1 here means that all these are retained. A value of
 330    *            0 here requires generated dihedrals be removed if
 331    *              * there are any dihedrals on the same central atoms
 332    *                specified in the residue topology, or
 333    *              * there are other identical generated dihedrals
 334    *                sharing the same central atoms, or
 335    *              * there are other generated dihedrals sharing the
 336    *                same central bond that have fewer hydrogen atoms
 337    * Column 6: Number of bonded neighbors to exclude.
 338    * Column 7: Generate 1,4 interactions between two hydrogen atoms
 339    * Column 8: Remove proper dihedrals if centered on the same bond
 340    *           as an improper dihedral
 341    */
 342   get_a_line(in,line,STRLEN);
 343   if (!get_header(line,header))
 344     gmx_fatal(FARGS,"in .rtp file at line:\n%s\n",line);
 345   if (gmx_strncasecmp("bondedtypes",header,5)==0) {
 346     get_a_line(in,line,STRLEN);
 347     if ((nparam=sscanf(line,"%d %d %d %d %d %d %d %d",
 348                        &header_settings->rb[ebtsBONDS].type,&header_settings->rb[ebtsANGLES].type,
 349                        &header_settings->rb[ebtsPDIHS].type,&header_settings->rb[ebtsIDIHS].type,
 350                        &dum1,&header_settings->nrexcl,&dum2,&dum3)) < 4 )
 351     {
 352       gmx_fatal(FARGS,"need 4 to 8 parameters in the header of .rtp file %s at line:\n%s\n", rrdb, line);
 353     }
 354     header_settings->bKeepAllGeneratedDihedrals    = (dum1 != 0);
 355     header_settings->bGenerateHH14Interactions     = (dum2 != 0);
 356     header_settings->bRemoveDihedralIfWithImproper = (dum3 != 0);
 357     get_a_line(in,line,STRLEN);
 358     if(nparam<5) {
 359       fprintf(stderr,"Using default: not generating all possible dihedrals\n");
 360       header_settings->bKeepAllGeneratedDihedrals = FALSE;
 361     }
 362     if(nparam<6) {
 363       fprintf(stderr,"Using default: excluding 3 bonded neighbors\n");
 364       header_settings->nrexcl = 3;
 365     }
 366     if(nparam<7) {
 367       fprintf(stderr,"Using default: generating 1,4 H--H interactions\n");
 368       header_settings->bGenerateHH14Interactions = TRUE;
 369     }
 370     if(nparam<8) {
 371       fprintf(stderr,"Using default: removing proper dihedrals found on the same bond as a proper dihedral\n");
 372       header_settings->bRemoveDihedralIfWithImproper = TRUE;
 373     }
 374   } else {
 375     fprintf(stderr,
 376             "Reading .rtp file without '[ bondedtypes ]' directive,\n"
 377             "Will proceed as if the entry was:\n");
 378     print_resall_header(stderr, header_settings);
 379   }
 380   /* We don't know the current size of rrtp, but simply realloc immediately */
 381   nrtp = *nrtpptr;
 382   rrtp = *rtp;
 383   maxrtp = nrtp;
 384   while (!feof(in)) {
 385     if (nrtp >= maxrtp) {
 386       maxrtp+=100;
 387       srenew(rrtp,maxrtp);
 388     }
 389
 390     /* Initialise rtp entry structure */
 391     rrtp[nrtp] = *header_settings;
 392     if (!get_header(line,header))
 393       gmx_fatal(FARGS,"in .rtp file at line:\n%s\n",line);
 394     rrtp[nrtp].resname  = strdup(header);
 395     rrtp[nrtp].filebase = strdup(filebase);
 396
 397     get_a_line(in,line,STRLEN);
 398     bError=FALSE;
 399     bNextResidue=FALSE;
 400     do {
 401       if (!get_header(line,header)) {
 402         bError = TRUE;
 403       } else {
 404         bt = get_bt(header);
 405         if (bt != NOTSET) {
 406           /* header is an bonded directive */
 407           bError = !read_bondeds(bt,in,line,&rrtp[nrtp]);
 408         } else if (gmx_strncasecmp("atoms",header,5) == 0) {
 409           /* header is the atoms directive */
 410           bError = !read_atoms(in,line,&(rrtp[nrtp]),tab,atype);
 411         } else {
 412           /* else header must be a residue name */
 413           bNextResidue = TRUE;
 414         }
 415       }
 416       if (bError)
 417         gmx_fatal(FARGS,"in .rtp file in residue %s at line:\n%s\n",
 418                     rrtp[nrtp].resname,line);
 419     } while (!feof(in) && !bNextResidue);
 420
 421     if (rrtp[nrtp].natom == 0)
 422       gmx_fatal(FARGS,"No atoms found in .rtp file in residue %s\n",
 423                   rrtp[nrtp].resname);
 424     if (debug) {
 425       fprintf(debug,"%3d %5s %5d",
 426               nrtp+1,rrtp[nrtp].resname,rrtp[nrtp].natom);
 427       for(i=0; i<ebtsNR; i++)
 428         fprintf(debug," %10d",rrtp[nrtp].rb[i].nb);
 429       fprintf(debug,"\n");
 430     }
 431
 432     firstrtp = -1;
 433     for(i=0; i<nrtp; i++) {
 434         if (gmx_strcasecmp(rrtp[i].resname,rrtp[nrtp].resname) == 0) {
 435             firstrtp = i;
 436         }
 437     }
 438
 439     if (firstrtp == -1) {
 440         nrtp++;
 441         fprintf(stderr,"\rResidue %d",nrtp);
 442     } else {
 443         if (firstrtp >= *nrtpptr)
 444         {
 445             gmx_fatal(FARGS,"Found a second entry for '%s' in '%s'",
 446                       rrtp[nrtp].resname,rrdb);
 447         }
 448         if (bAllowOverrideRTP)
 449         {
 450             fprintf(stderr,"\n");
 451             fprintf(stderr,"Found another rtp entry for '%s' in '%s', ignoring this entry and keeping the one from '%s.rtp'\n",
 452                     rrtp[nrtp].resname,rrdb,rrtp[firstrtp].filebase);
 453             /* We should free all the data for this entry.
 454              * The current code gives a lot of dangling pointers.
 455              */
 456             clear_t_restp(&rrtp[nrtp]);
 457         }
 458         else
 459         {
 460             gmx_fatal(FARGS,"Found rtp entries for '%s' in both '%s' and '%s'. If you want the first definition to override the second one, set the -rtpo option of pdb2gmx.",rrtp[nrtp].resname,rrtp[firstrtp].filebase,rrdb);
 461         }
 462     }
 463   }
 464   ffclose(in);
 465   /* give back unused memory */
 466   srenew(rrtp,nrtp);
 467
 468   fprintf(stderr,"\nSorting it all out...\n");
 469   qsort(rrtp,nrtp,(size_t)sizeof(rrtp[0]),comprtp);
 470
 471   check_rtp(nrtp,rrtp,rrdb);
 472
 473   *nrtpptr = nrtp;
 474   *rtp     = rrtp;
 475 }
 476
 477 /************************************************************
 478  *
 479  *                  SEARCH   ROUTINES
 480  *
 481  ***********************************************************/
 482 static gmx_bool is_sign(char c)
 483 {
 484     return (c == '+' || c == '-');
 485 }
 486
 487 /* Compares if the strins match except for a sign at the end
 488  * of (only) one of the two.
 489  */
 490 static int neq_str_sign(const char *a1,const char *a2)
 491 {
 492     int l1,l2,lm;
 493
 494     l1 = (int)strlen(a1);
 495     l2 = (int)strlen(a2);
 496     lm = min(l1,l2);
 497
 498     if (lm >= 1 &&
 499         ((l1 == l2+1 && is_sign(a1[l1-1])) ||
 500          (l2 == l1+1 && is_sign(a2[l2-1]))) &&
 501         gmx_strncasecmp(a1,a2,lm) == 0)
 502     {
 503         return lm;
 504     }
 505     else
 506     {
 507         return 0;
 508     }
 509 }
 510
 511 char *search_rtp(const char *key,int nrtp,t_restp rtp[])
 512 {
 513     int  i,n,nbest,best,besti;
 514     char bestbuf[STRLEN];
 515
 516     nbest =  0;
 517     besti = -1;
 518     /* We want to match at least one character */
 519     best  =  1;
 520     for(i=0; (i<nrtp); i++)
 521     {
 522         if (gmx_strcasecmp(key,rtp[i].resname) == 0)
 523         {
 524              besti = i;
 525              nbest = 1;
 526              break;
 527         }
 528         else
 529         {
 530             /* Allow a mismatch of at most a sign character (with warning) */
 531             n = neq_str_sign(key,rtp[i].resname);
 532             if (n >= best &&
 533                 n+1 >= (int)strlen(key) &&
 534                 n+1 >= (int)strlen(rtp[i].resname))
 535             {
 536                 if (n == best)
 537                 {
 538                     if (nbest == 1)
 539                     {
 540                         strcpy(bestbuf,rtp[besti].resname);
 541                     }
 542                     if (nbest >= 1)
 543                     {
 544                         strcat(bestbuf," or ");
 545                         strcat(bestbuf,rtp[i].resname);
 546                     }
 547                 }
 548                 else
 549                 {
 550                     nbest = 0;
 551                 }
 552                 besti = i;
 553                 best  = n;
 554                 nbest++;
 555             }
 556         }
 557     }
 558     if (nbest > 1)
 559     {
 560         gmx_fatal(FARGS,"Residue '%s' not found in residue topology database, looks a bit like %s",key,bestbuf);
 561     }
 562     else if (besti == -1)
 563     {
 564         gmx_fatal(FARGS,"Residue '%s' not found in residue topology database",key);
 565     }
 566     if (gmx_strcasecmp(rtp[besti].resname,key) != 0)
 567     {
 568         fprintf(stderr,
 569                 "\nWARNING: '%s' not found in residue topology database, "
 570                 "trying to use '%s'\n\n", key, rtp[besti].resname);
 571     }
 572
 573     return rtp[besti].resname;
 574 }
 575
 576 t_restp *get_restp(const char *rtpname,int nrtp,t_restp rtp[])
 577 {
 578     int  i;
 579
 580     i = 0;
 581     while (i < nrtp && gmx_strcasecmp(rtpname,rtp[i].resname) != 0)
 582     {
 583         i++;
 584     }
 585     if (i >= nrtp)
 586     {
 587         /* This should never happen, since search_rtp should have been called
 588          * before calling get_restp.
 589          */
 590         gmx_fatal(FARGS,"Residue type '%s' not found in residue topology database",rtpname);
 591     }
 592
 593     return &rtp[i];
 594 }