src/kernel/resall.c

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