src/kernel/pdb2top.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,2013, 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 <stdio.h>
  43 #include <math.h>
  44 #include <ctype.h>
  45
  46 #ifdef GMX_NATIVE_WINDOWS
  47 #include <direct.h>
  48 #include <io.h>
  49 #endif
  50
  51
  52 #include "vec.h"
  53 #include "copyrite.h"
  54 #include "smalloc.h"
  55 #include "macros.h"
  56 #include "symtab.h"
  57 #include "futil.h"
  58 #include "statutil.h"
  59 #include "gmx_fatal.h"
  60 #include "pdb2top.h"
  61 #include "gpp_nextnb.h"
  62 #include "topdirs.h"
  63 #include "toputil.h"
  64 #include "h_db.h"
  65 #include "pgutil.h"
  66 #include "resall.h"
  67 #include "topio.h"
  68 #include "string2.h"
  69 #include "physics.h"
  70 #include "pdbio.h"
  71 #include "gen_ad.h"
  72 #include "filenm.h"
  73 #include "index.h"
  74 #include "gen_vsite.h"
  75 #include "add_par.h"
  76 #include "toputil.h"
  77 #include "fflibutil.h"
  78 #include "strdb.h"
  79
  80 /* this must correspond to enum in pdb2top.h */
  81 const char *hh[ehisNR]   = { "HISD", "HISE", "HISH", "HIS1" };
  82
  83 static int missing_atoms(t_restp *rp, int resind,t_atoms *at, int i0, int i)
  84 {
  85     int  j,k,nmiss;
  86     char *name;
  87     gmx_bool bFound, bRet;
  88
  89     nmiss = 0;
  90     for (j=0; j<rp->natom; j++)
  91     {
  92         name=*(rp->atomname[j]);
  93         bFound=FALSE;
  94         for (k=i0; k<i; k++)
  95         {
  96             bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]),name));
  97         }
  98         if (!bFound)
  99         {
 100             nmiss++;
 101             fprintf(stderr,"\nWARNING: "
 102                     "atom %s is missing in residue %s %d in the pdb file\n",
 103                     name,*(at->resinfo[resind].name),at->resinfo[resind].nr);
 104             if (name[0]=='H' || name[0]=='h')
 105             {
 106                 fprintf(stderr,"         You might need to add atom %s to the hydrogen database of building block %s\n"
 107                         "         in the file %s.hdb (see the manual)\n",
 108                         name,*(at->resinfo[resind].rtp),rp->filebase);
 109             }
 110             fprintf(stderr,"\n");
 111         }
 112     }
 113
 114     return nmiss;
 115 }
 116
 117 gmx_bool is_int(double x)
 118 {
 119   const double tol = 1e-4;
 120   int   ix;
 121
 122   if (x < 0)
 123     x=-x;
 124   ix=gmx_nint(x);
 125
 126   return (fabs(x-ix) < tol);
 127 }
 128
 129 static void swap_strings(char **s,int i,int j)
 130 {
 131     char *tmp;
 132
 133     tmp  = s[i];
 134     s[i] = s[j];
 135     s[j] = tmp;
 136 }
 137
 138 void
 139 choose_ff(const char *ffsel,
 140           char *forcefield, int ff_maxlen,
 141           char *ffdir, int ffdir_maxlen)
 142 {
 143     int  nff;
 144     char **ffdirs,**ffs,**ffs_dir,*ptr;
 145     int  i,j,sel,cwdsel,nfound;
 146     char buf[STRLEN],**desc;
 147     FILE *fp;
 148     char *pret;
 149
 150     nff = fflib_search_file_in_dirend(fflib_forcefield_itp(),
 151                                       fflib_forcefield_dir_ext(),
 152                                       &ffdirs);
 153
 154     if (nff == 0)
 155     {
 156         gmx_fatal(FARGS,"No force fields found (files with name '%s' in subdirectories ending on '%s')",
 157                   fflib_forcefield_itp(),fflib_forcefield_dir_ext());
 158     }
 159
 160     /* Replace with unix path separators */
 161     if(DIR_SEPARATOR!='/')
 162     {
 163         for(i=0;i<nff;i++)
 164         {
 165             while( (ptr=strchr(ffdirs[i],DIR_SEPARATOR))!=NULL )
 166             {
 167                 *ptr='/';
 168             }
 169         }
 170     }
 171
 172     /* Store the force field names in ffs */
 173     snew(ffs,nff);
 174     snew(ffs_dir,nff);
 175     for(i=0; i<nff; i++)
 176     {
 177         /* Remove the path from the ffdir name - use our unix standard here! */
 178         ptr = strrchr(ffdirs[i],'/');
 179         if (ptr == NULL)
 180         {
 181             ffs[i] = strdup(ffdirs[i]);
 182             ffs_dir[i] = low_gmxlibfn(ffdirs[i],FALSE,FALSE);
 183             if (ffs_dir[i] == NULL)
 184             {
 185                 gmx_fatal(FARGS,"Can no longer find file '%s'",ffdirs[i]);
 186             }
 187         }
 188         else
 189         {
 190             ffs[i] = strdup(ptr+1);
 191             ffs_dir[i] = strdup(ffdirs[i]);
 192         }
 193         ffs_dir[i][strlen(ffs_dir[i])-strlen(ffs[i])-1] = '\0';
 194         /* Remove the extension from the ffdir name */
 195         ffs[i][strlen(ffs[i])-strlen(fflib_forcefield_dir_ext())] = '\0';
 196     }
 197
 198     if (ffsel != NULL)
 199     {
 200         sel     = -1;
 201         cwdsel  = -1;
 202         nfound  = 0;
 203         for(i=0; i<nff; i++)
 204         {
 205             if ( strcmp(ffs[i],ffsel)==0 )
 206             {
 207                 /* Matching ff name */
 208                 sel = i;
 209                 nfound++;
 210
 211                 if( strncmp(ffs_dir[i],".",1)==0 )
 212                 {
 213                     cwdsel = i;
 214                 }
 215             }
 216         }
 217
 218         if(cwdsel != -1)
 219         {
 220             sel = cwdsel;
 221         }
 222
 223         if(nfound>1)
 224         {
 225             if(cwdsel!=-1)
 226             {
 227                 fprintf(stderr,
 228                         "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
 229                         "current directory. Use interactive selection (not the -ff option) if\n"
 230                         "you would prefer a different one.\n",ffsel,nfound);
 231             }
 232             else
 233             {
 234                 gmx_fatal(FARGS,
 235                           "Force field '%s' occurs in %d places, but not in the current directory.\n"
 236                           "Run without the -ff switch and select the force field interactively.",ffsel,nfound);
 237             }
 238         }
 239         else if (nfound==0)
 240         {
 241             gmx_fatal(FARGS,"Could not find force field '%s' in current directory, install tree or GMXDATA path.",ffsel);
 242         }
 243     }
 244     else if (nff > 1)
 245     {
 246         snew(desc,nff);
 247         for(i=0; (i<nff); i++)
 248         {
 249             sprintf(buf,"%s%c%s%s%c%s",
 250                     ffs_dir[i],DIR_SEPARATOR,
 251                     ffs[i],fflib_forcefield_dir_ext(),DIR_SEPARATOR,
 252                     fflib_forcefield_doc());
 253             if (gmx_fexist(buf))
 254             {
 255                 /* We don't use fflib_open, because we don't want printf's */
 256                 fp = ffopen(buf,"r");
 257                 snew(desc[i],STRLEN);
 258                 get_a_line(fp,desc[i],STRLEN);
 259                 ffclose(fp);
 260             }
 261             else
 262             {
 263                 desc[i] = strdup(ffs[i]);
 264             }
 265         }
 266         /* Order force fields from the same dir alphabetically
 267          * and put deprecated force fields at the end.
 268          */
 269         for(i=0; (i<nff); i++)
 270         {
 271             for(j=i+1; (j<nff); j++)
 272             {
 273                 if (strcmp(ffs_dir[i],ffs_dir[j]) == 0 &&
 274                     ((desc[i][0] == '[' && desc[j][0] != '[') ||
 275                      ((desc[i][0] == '[' || desc[j][0] != '[') &&
 276                       gmx_strcasecmp(desc[i],desc[j]) > 0)))
 277                 {
 278                     swap_strings(ffdirs,i,j);
 279                     swap_strings(ffs   ,i,j);
 280                     swap_strings(desc  ,i,j);
 281                 }
 282             }
 283         }
 284
 285         printf("\nSelect the Force Field:\n");
 286         for(i=0; (i<nff); i++)
 287         {
 288             if (i == 0 || strcmp(ffs_dir[i-1],ffs_dir[i]) != 0)
 289             {
 290                 if( strcmp(ffs_dir[i],".")==0 )
 291                 {
 292                     printf("From current directory:\n");
 293                 }
 294                 else
 295                 {
 296                     printf("From '%s':\n",ffs_dir[i]);
 297                 }
 298             }
 299             printf("%2d: %s\n",i+1,desc[i]);
 300             sfree(desc[i]);
 301         }
 302         sfree(desc);
 303
 304         do
 305         {
 306             pret = fgets(buf,STRLEN,stdin);
 307
 308             if (pret != NULL)
 309             {
 310                 sscanf(buf,"%d",&sel);
 311                 sel--;
 312             }
 313         }
 314         while ( pret==NULL || (sel < 0) || (sel >= nff));
 315
 316         /* Check for a current limitation of the fflib code.
 317          * It will always read from the first ff directory in the list.
 318          * This check assumes that the order of ffs matches the order
 319          * in which fflib_open searches ff library files.
 320          */
 321         for(i=0; i<sel; i++)
 322         {
 323             if (strcmp(ffs[i],ffs[sel]) == 0)
 324             {
 325                 gmx_fatal(FARGS,"Can only select the first of multiple force field entries with directory name '%s%s' in the list. If you want to use the next entry, run pdb2gmx in a different directory or rename or move the force field directory present in the current working directory.",
 326                           ffs[sel],fflib_forcefield_dir_ext());
 327             }
 328         }
 329     }
 330     else
 331     {
 332         sel = 0;
 333     }
 334
 335     if (strlen(ffs[sel]) >= (size_t)ff_maxlen)
 336     {
 337         gmx_fatal(FARGS,"Length of force field name (%d) >= maxlen (%d)",
 338                   strlen(ffs[sel]),ff_maxlen);
 339     }
 340     strcpy(forcefield,ffs[sel]);
 341
 342     if (strlen(ffdirs[sel]) >= (size_t)ffdir_maxlen)
 343     {
 344         gmx_fatal(FARGS,"Length of force field dir (%d) >= maxlen (%d)",
 345                   strlen(ffdirs[sel]),ffdir_maxlen);
 346     }
 347     strcpy(ffdir,ffdirs[sel]);
 348
 349     for(i=0; (i<nff); i++)
 350     {
 351         sfree(ffdirs[i]);
 352         sfree(ffs[i]);
 353         sfree(ffs_dir[i]);
 354     }
 355     sfree(ffdirs);
 356     sfree(ffs);
 357     sfree(ffs_dir);
 358 }
 359
 360 void choose_watermodel(const char *wmsel,const char *ffdir,
 361                        char **watermodel)
 362 {
 363     const char *fn_watermodels="watermodels.dat";
 364     char fn_list[STRLEN];
 365     FILE *fp;
 366     char buf[STRLEN];
 367     int  nwm,sel,i;
 368     char **model;
 369     char *pret;
 370
 371     if (strcmp(wmsel,"none") == 0)
 372     {
 373         *watermodel = NULL;
 374
 375         return;
 376     }
 377     else if (strcmp(wmsel,"select") != 0)
 378     {
 379         *watermodel = strdup(wmsel);
 380
 381         return;
 382     }
 383
 384     sprintf(fn_list,"%s%c%s",ffdir,DIR_SEPARATOR,fn_watermodels);
 385
 386     if (!fflib_fexist(fn_list))
 387     {
 388         fprintf(stderr,"No file '%s' found, will not include a water model\n",
 389                 fn_watermodels);
 390         *watermodel = NULL;
 391
 392         return;
 393     }
 394
 395     fp = fflib_open(fn_list);
 396     printf("\nSelect the Water Model:\n");
 397     nwm = 0;
 398     model = NULL;
 399     while (get_a_line(fp,buf,STRLEN))
 400     {
 401         srenew(model,nwm+1);
 402         snew(model[nwm],STRLEN);
 403         sscanf(buf,"%s%n",model[nwm],&i);
 404         if (i > 0)
 405         {
 406             ltrim(buf+i);
 407             fprintf(stderr,"%2d: %s\n",nwm+1,buf+i);
 408             nwm++;
 409         }
 410         else
 411         {
 412             sfree(model[nwm]);
 413         }
 414     }
 415     ffclose(fp);
 416     fprintf(stderr,"%2d: %s\n",nwm+1,"None");
 417
 418     do
 419     {
 420         pret = fgets(buf,STRLEN,stdin);
 421
 422         if (pret != NULL)
 423         {
 424             sscanf(buf,"%d",&sel);
 425             sel--;
 426         }
 427     }
 428     while (pret == NULL || sel < 0 || sel > nwm);
 429
 430     if (sel == nwm)
 431     {
 432         *watermodel = NULL;
 433     }
 434     else
 435     {
 436         *watermodel = strdup(model[sel]);
 437     }
 438
 439     for(i=0; i<nwm; i++)
 440     {
 441         sfree(model[i]);
 442     }
 443     sfree(model);
 444 }
 445
 446 static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
 447                      t_restp restp[], gmx_residuetype_t rt)
 448 {
 449   int     i,j,prevresind,resind,i0,prevcg,cg,curcg;
 450   char    *name;
 451   gmx_bool    bProt, bNterm;
 452   double  qt;
 453   int     nmissat;
 454
 455   nmissat = 0;
 456
 457   resind=-1;
 458   bProt=FALSE;
 459   bNterm=FALSE;
 460   i0=0;
 461   snew(*cgnr,at->nr);
 462   qt=0;
 463   prevcg=NOTSET;
 464   curcg=0;
 465   cg=-1;
 466   j=NOTSET;
 467
 468   for(i=0; (i<at->nr); i++) {
 469     prevresind=resind;
 470     if (at->atom[i].resind != resind) {
 471       resind = at->atom[i].resind;
 472       bProt = gmx_residuetype_is_protein(rt,*(at->resinfo[resind].name));
 473       bNterm=bProt && (resind == 0);
 474       if (resind > 0) {
 475           nmissat += missing_atoms(&restp[prevresind],prevresind,at,i0,i);
 476       }
 477       i0=i;
 478     }
 479     if (at->atom[i].m == 0) {
 480       if (debug)
 481         fprintf(debug,"atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
 482                 i+1,*(at->atomname[i]),curcg,prevcg,
 483                 j==NOTSET ? NOTSET : restp[resind].cgnr[j]);
 484       qt=0;
 485       prevcg=cg;
 486       name=*(at->atomname[i]);
 487       j=search_jtype(&restp[resind],name,bNterm);
 488       at->atom[i].type = restp[resind].atom[j].type;
 489       at->atom[i].q    = restp[resind].atom[j].q;
 490       at->atom[i].m    = get_atomtype_massA(restp[resind].atom[j].type,
 491                                             atype);
 492       cg = restp[resind].cgnr[j];
 493       /* A charge group number -1 signals a separate charge group
 494        * for this atom.
 495        */
 496       if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) ) {
 497           curcg++;
 498       }
 499     } else {
 500       if (debug)
 501         fprintf(debug,"atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
 502                 i+1,*(at->atomname[i]),curcg,qt,is_int(qt));
 503       cg=-1;
 504       if (is_int(qt)) {
 505         qt=0;
 506         curcg++;
 507       }
 508       qt+=at->atom[i].q;
 509     }
 510     (*cgnr)[i]=curcg;
 511     at->atom[i].typeB = at->atom[i].type;
 512     at->atom[i].qB    = at->atom[i].q;
 513     at->atom[i].mB    = at->atom[i].m;
 514   }
 515   nmissat += missing_atoms(&restp[resind],resind,at,i0,i);
 516
 517   return nmissat;
 518 }
 519
 520 static void print_top_heavy_H(FILE *out, real mHmult)
 521 {
 522   if (mHmult == 2.0)
 523     fprintf(out,"; Using deuterium instead of hydrogen\n\n");
 524   else if (mHmult == 4.0)
 525     fprintf(out,"#define HEAVY_H\n\n");
 526   else if (mHmult != 1.0)
 527     fprintf(stderr,"WARNING: unsupported proton mass multiplier (%g) "
 528             "in pdb2top\n",mHmult);
 529 }
 530
 531 void print_top_comment(FILE *out,
 532                        const char *filename,
 533                        const char *generator,
 534                        const char *ffdir,
 535                        gmx_bool bITP)
 536 {
 537   char tmp[256];
 538   char ffdir_parent[STRLEN];
 539   char *p;
 540
 541   nice_header(out,filename);
 542   fprintf(out,";\tThis is a %s topology file\n;\n",bITP ? "include" : "standalone");
 543   fprintf(out,";\tIt was generated using program:\n;\t%s\n;\n",
 544           (NULL == generator) ? "unknown" : generator);
 545   fprintf(out,";\tCommand line was:\n;\t%s\n;\n",command_line());
 546
 547   if(strchr(ffdir,'/')==NULL)
 548   {
 549       fprintf(out,";\tForce field was read from the standard Gromacs share directory.\n;\n\n");
 550   }
 551   else if(ffdir[0]=='.')
 552   {
 553       fprintf(out,";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
 554   }
 555   else
 556   {
 557       strncpy(ffdir_parent,ffdir,STRLEN-1);
 558       ffdir_parent[STRLEN-1]='\0'; /*make sure it is 0-terminated even for long string*/
 559       p=strrchr(ffdir_parent,'/');
 560
 561       *p='\0';
 562
 563       fprintf(out,
 564               ";\tForce field data was read from:\n"
 565               ";\t%s\n"
 566               ";\n"
 567               ";\tNote:\n"
 568               ";\tThis might be a non-standard force field location. When you use this topology, the\n"
 569               ";\tforce field must either be present in the current directory, or the location\n"
 570               ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
 571               ffdir_parent);
 572   }
 573 }
 574
 575 void print_top_header(FILE *out,const char *filename,
 576                       const char *title,gmx_bool bITP,const char *ffdir,real mHmult)
 577 {
 578     const char *p;
 579
 580     print_top_comment(out,filename,title,ffdir,bITP);
 581
 582     print_top_heavy_H(out, mHmult);
 583     fprintf(out,"; Include forcefield parameters\n");
 584
 585     p=strrchr(ffdir,'/');
 586     p = (ffdir[0]=='.' || p==NULL) ? ffdir : p+1;
 587
 588     fprintf(out,"#include \"%s/%s\"\n\n",p,fflib_forcefield_itp());
 589 }
 590
 591 static void print_top_posre(FILE *out,const char *pr)
 592 {
 593   fprintf(out,"; Include Position restraint file\n");
 594   fprintf(out,"#ifdef POSRES\n");
 595   fprintf(out,"#include \"%s\"\n",pr);
 596   fprintf(out,"#endif\n\n");
 597 }
 598
 599 static void print_top_water(FILE *out,const char *ffdir,const char *water)
 600 {
 601   const char *p;
 602   char  buf[STRLEN];
 603
 604   fprintf(out,"; Include water topology\n");
 605
 606   p=strrchr(ffdir,'/');
 607   p = (ffdir[0]=='.' || p==NULL) ? ffdir : p+1;
 608   fprintf(out,"#include \"%s/%s.itp\"\n",p,water);
 609
 610   fprintf(out,"\n");
 611   fprintf(out,"#ifdef POSRES_WATER\n");
 612   fprintf(out,"; Position restraint for each water oxygen\n");
 613   fprintf(out,"[ position_restraints ]\n");
 614   fprintf(out,";%3s %5s %9s %10s %10s\n","i","funct","fcx","fcy","fcz");
 615   fprintf(out,"%4d %4d %10g %10g %10g\n",1,1,1000.0,1000.0,1000.0);
 616   fprintf(out,"#endif\n");
 617   fprintf(out,"\n");
 618
 619   sprintf(buf,"%s/ions.itp",p);
 620
 621   if (fflib_fexist(buf))
 622   {
 623     fprintf(out,"; Include topology for ions\n");
 624     fprintf(out,"#include \"%s\"\n",buf);
 625     fprintf(out,"\n");
 626   }
 627 }
 628
 629 static void print_top_system(FILE *out, const char *title)
 630 {
 631   fprintf(out,"[ %s ]\n",dir2str(d_system));
 632   fprintf(out,"; Name\n");
 633   fprintf(out,"%s\n\n",title[0]?title:"Protein");
 634 }
 635
 636 void print_top_mols(FILE *out,
 637                     const char *title, const char *ffdir, const char *water,
 638                     int nincl, char **incls, int nmol, t_mols *mols)
 639 {
 640   int  i;
 641   char *incl;
 642
 643   if (nincl>0) {
 644     fprintf(out,"; Include chain topologies\n");
 645     for (i=0; (i<nincl); i++) {
 646         incl = strrchr(incls[i],DIR_SEPARATOR);
 647         if (incl == NULL) {
 648             incl = incls[i];
 649         } else {
 650             /* Remove the path from the include name */
 651             incl = incl + 1;
 652         }
 653       fprintf(out,"#include \"%s\"\n",incl);
 654     }
 655     fprintf(out,"\n");
 656   }
 657
 658     if (water)
 659     {
 660       print_top_water(out,ffdir,water);
 661     }
 662     print_top_system(out, title);
 663
 664   if (nmol) {
 665     fprintf(out,"[ %s ]\n",dir2str(d_molecules));
 666     fprintf(out,"; %-15s %5s\n","Compound","#mols");
 667     for (i=0; (i<nmol); i++)
 668       fprintf(out,"%-15s %5d\n",mols[i].name,mols[i].nr);
 669   }
 670 }
 671
 672 void write_top(FILE *out, char *pr,char *molname,
 673                t_atoms *at,gmx_bool bRTPresname,
 674                int bts[],t_params plist[],t_excls excls[],
 675                gpp_atomtype_t atype,int *cgnr, int nrexcl)
 676      /* NOTE: nrexcl is not the size of *excl! */
 677 {
 678   if (at && atype && cgnr) {
 679     fprintf(out,"[ %s ]\n",dir2str(d_moleculetype));
 680     fprintf(out,"; %-15s %5s\n","Name","nrexcl");
 681     fprintf(out,"%-15s %5d\n\n",molname?molname:"Protein",nrexcl);
 682
 683     print_atoms(out, atype, at, cgnr, bRTPresname);
 684     print_bondeds(out,at->nr,d_bonds,      F_BONDS,    bts[ebtsBONDS], plist);
 685     print_bondeds(out,at->nr,d_constraints,F_CONSTR,   0,              plist);
 686     print_bondeds(out,at->nr,d_constraints,F_CONSTRNC, 0,              plist);
 687     print_bondeds(out,at->nr,d_pairs,      F_LJ14,     0,              plist);
 688     print_excl(out,at->nr,excls);
 689     print_bondeds(out,at->nr,d_angles,     F_ANGLES,   bts[ebtsANGLES],plist);
 690     print_bondeds(out,at->nr,d_dihedrals,  F_PDIHS,    bts[ebtsPDIHS], plist);
 691     print_bondeds(out,at->nr,d_dihedrals,  F_IDIHS,    bts[ebtsIDIHS], plist);
 692     print_bondeds(out,at->nr,d_cmap,       F_CMAP,     bts[ebtsCMAP],  plist);
 693     print_bondeds(out,at->nr,d_polarization,F_POLARIZATION,   0,       plist);
 694     print_bondeds(out,at->nr,d_thole_polarization,F_THOLE_POL,0,       plist);
 695     print_bondeds(out,at->nr,d_vsites2,    F_VSITE2,   0,              plist);
 696     print_bondeds(out,at->nr,d_vsites3,    F_VSITE3,   0,              plist);
 697     print_bondeds(out,at->nr,d_vsites3,    F_VSITE3FD, 0,              plist);
 698     print_bondeds(out,at->nr,d_vsites3,    F_VSITE3FAD,0,              plist);
 699     print_bondeds(out,at->nr,d_vsites3,    F_VSITE3OUT,0,              plist);
 700     print_bondeds(out,at->nr,d_vsites4,    F_VSITE4FD, 0,              plist);
 701     print_bondeds(out,at->nr,d_vsites4,    F_VSITE4FDN, 0,             plist);
 702
 703     if (pr)
 704       print_top_posre(out,pr);
 705   }
 706 }
 707
 708 static atom_id search_res_atom(const char *type,int resind,
 709                    t_atoms *atoms,
 710                                const char *bondtype,gmx_bool bAllowMissing)
 711 {
 712   int i;
 713
 714   for(i=0; (i<atoms->nr); i++)
 715   {
 716     if (atoms->atom[i].resind == resind)
 717     {
 718       return search_atom(type,i,atoms,bondtype,bAllowMissing);
 719     }
 720   }
 721
 722   return NO_ATID;
 723 }
 724
 725 static void do_ssbonds(t_params *ps,t_atoms *atoms,
 726                        int nssbonds,t_ssbond *ssbonds,gmx_bool bAllowMissing)
 727 {
 728   int     i,ri,rj;
 729   atom_id ai,aj;
 730
 731   for(i=0; (i<nssbonds); i++) {
 732     ri = ssbonds[i].res1;
 733     rj = ssbonds[i].res2;
 734     ai = search_res_atom(ssbonds[i].a1,ri,atoms,
 735                          "special bond",bAllowMissing);
 736     aj = search_res_atom(ssbonds[i].a2,rj,atoms,
 737                          "special bond",bAllowMissing);
 738     if ((ai == NO_ATID) || (aj == NO_ATID))
 739       gmx_fatal(FARGS,"Trying to make impossible special bond (%s-%s)!",
 740                   ssbonds[i].a1,ssbonds[i].a2);
 741     add_param(ps,ai,aj,NULL,NULL);
 742   }
 743 }
 744
 745 static void at2bonds(t_params *psb, t_hackblock *hb,
 746                      t_atoms *atoms,
 747                      rvec x[],
 748                      real long_bond_dist, real short_bond_dist,
 749                      gmx_bool bAllowMissing)
 750 {
 751   int     resind,i,j,k;
 752   atom_id ai,aj;
 753   real    dist2, long_bond_dist2, short_bond_dist2;
 754   const char *ptr;
 755
 756   long_bond_dist2  = sqr(long_bond_dist);
 757   short_bond_dist2 = sqr(short_bond_dist);
 758
 759   if (debug)
 760     ptr = "bond";
 761   else
 762     ptr = "check";
 763
 764   fprintf(stderr,"Making bonds...\n");
 765   i=0;
 766   for(resind=0; (resind < atoms->nres) && (i<atoms->nr); resind++) {
 767     /* add bonds from list of bonded interactions */
 768     for(j=0; j < hb[resind].rb[ebtsBONDS].nb; j++) {
 769       /* Unfortunately we can not issue errors or warnings
 770        * for missing atoms in bonds, as the hydrogens and terminal atoms
 771        * have not been added yet.
 772        */
 773       ai=search_atom(hb[resind].rb[ebtsBONDS].b[j].AI,i,atoms,
 774                      ptr,TRUE);
 775       aj=search_atom(hb[resind].rb[ebtsBONDS].b[j].AJ,i,atoms,
 776                      ptr,TRUE);
 777       if (ai != NO_ATID && aj != NO_ATID) {
 778           dist2 = distance2(x[ai],x[aj]);
 779           if (dist2 > long_bond_dist2 )
 780           {
 781               fprintf(stderr,"Warning: Long Bond (%d-%d = %g nm)\n",
 782                       ai+1,aj+1,sqrt(dist2));
 783           }
 784           else if (dist2 < short_bond_dist2 )
 785           {
 786               fprintf(stderr,"Warning: Short Bond (%d-%d = %g nm)\n",
 787                       ai+1,aj+1,sqrt(dist2));
 788           }
 789           add_param(psb,ai,aj,NULL,hb[resind].rb[ebtsBONDS].b[j].s);
 790       }
 791     }
 792     /* add bonds from list of hacks (each added atom gets a bond) */
 793     while( (i<atoms->nr) && (atoms->atom[i].resind == resind) ) {
 794       for(j=0; j < hb[resind].nhack; j++)
 795         if ( ( hb[resind].hack[j].tp > 0 ||
 796                hb[resind].hack[j].oname==NULL ) &&
 797              strcmp(hb[resind].hack[j].AI,*(atoms->atomname[i])) == 0 ) {
 798           switch(hb[resind].hack[j].tp) {
 799           case 9:          /* COOH terminus */
 800             add_param(psb,i,i+1,NULL,NULL);     /* C-O  */
 801             add_param(psb,i,i+2,NULL,NULL);     /* C-OA */
 802             add_param(psb,i+2,i+3,NULL,NULL);   /* OA-H */
 803             break;
 804           default:
 805             for(k=0; (k<hb[resind].hack[j].nr); k++)
 806               add_param(psb,i,i+k+1,NULL,NULL);
 807           }
 808         }
 809       i++;
 810     }
 811     /* we're now at the start of the next residue */
 812   }
 813 }
 814
 815 static int pcompar(const void *a, const void *b)
 816 {
 817   t_param *pa,*pb;
 818   int     d;
 819   pa=(t_param *)a;
 820   pb=(t_param *)b;
 821
 822   d = pa->AI - pb->AI;
 823   if (d == 0)
 824     d = pa->AJ - pb->AJ;
 825   if (d == 0)
 826     return strlen(pb->s) - strlen(pa->s);
 827   else
 828     return d;
 829 }
 830
 831 static void clean_bonds(t_params *ps)
 832 {
 833   int     i,j;
 834   atom_id a;
 835
 836   if (ps->nr > 0) {
 837     /* swap atomnumbers in bond if first larger than second: */
 838     for(i=0; (i<ps->nr); i++)
 839       if ( ps->param[i].AJ < ps->param[i].AI ) {
 840         a = ps->param[i].AI;
 841         ps->param[i].AI = ps->param[i].AJ;
 842         ps->param[i].AJ = a;
 843       }
 844
 845     /* Sort bonds */
 846     qsort(ps->param,ps->nr,(size_t)sizeof(ps->param[0]),pcompar);
 847
 848     /* remove doubles, keep the first one always. */
 849     j = 1;
 850     for(i=1; (i<ps->nr); i++) {
 851       if ((ps->param[i].AI != ps->param[j-1].AI) ||
 852           (ps->param[i].AJ != ps->param[j-1].AJ) ) {
 853         if (j != i) {
 854           cp_param(&(ps->param[j]),&(ps->param[i]));
 855         }
 856         j++;
 857       }
 858     }
 859     fprintf(stderr,"Number of bonds was %d, now %d\n",ps->nr,j);
 860     ps->nr=j;
 861   }
 862   else
 863     fprintf(stderr,"No bonds\n");
 864 }
 865
 866 void print_sums(t_atoms *atoms, gmx_bool bSystem)
 867 {
 868   double m,qtot;
 869   int    i;
 870   const char *where;
 871
 872   if (bSystem)
 873     where=" in system";
 874   else
 875     where="";
 876
 877   m=0;
 878   qtot=0;
 879   for(i=0; (i<atoms->nr); i++) {
 880     m+=atoms->atom[i].m;
 881     qtot+=atoms->atom[i].q;
 882   }
 883
 884   fprintf(stderr,"Total mass%s %.3f a.m.u.\n",where,m);
 885   fprintf(stderr,"Total charge%s %.3f e\n",where,qtot);
 886 }
 887
 888 static void check_restp_type(const char *name,int t1,int t2)
 889 {
 890     if (t1 != t2)
 891     {
 892         gmx_fatal(FARGS,"Residues in one molecule have a different '%s' type: %d and %d",name,t1,t2);
 893     }
 894 }
 895
 896 static void check_restp_types(t_restp *r0,t_restp *r1)
 897 {
 898     int i;
 899
 900     check_restp_type("all dihedrals",r0->bKeepAllGeneratedDihedrals,r1->bKeepAllGeneratedDihedrals);
 901     check_restp_type("nrexcl",r0->nrexcl,r1->nrexcl);
 902     check_restp_type("HH14",r0->bGenerateHH14Interactions,r1->bGenerateHH14Interactions);
 903     check_restp_type("remove dihedrals",r0->bRemoveDihedralIfWithImproper,r1->bRemoveDihedralIfWithImproper);
 904
 905     for(i=0; i<ebtsNR; i++)
 906     {
 907         check_restp_type(btsNames[i],r0->rb[i].type,r1->rb[i].type);
 908     }
 909 }
 910
 911 void add_atom_to_restp(t_restp *restp,int resnr,int at_start,const t_hack *hack)
 912 {
 913     char buf[STRLEN];
 914     int  k;
 915     const char *Hnum="123456";
 916
 917     /*if (debug)
 918     {
 919         fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
 920                 hack->nname,
 921                 *restp->atomname[at_start], resnr, restp->resname);
 922                 }*/
 923     strcpy(buf, hack->nname);
 924     buf[strlen(buf)+1]='\0';
 925     if ( hack->nr > 1 )
 926     {
 927         buf[strlen(buf)]='-';
 928     }
 929     /* make space */
 930     restp->natom += hack->nr;
 931     srenew(restp->atom,     restp->natom);
 932     srenew(restp->atomname, restp->natom);
 933     srenew(restp->cgnr,     restp->natom);
 934     /* shift the rest */
 935     for(k=restp->natom-1; k > at_start+hack->nr; k--)
 936     {
 937         restp->atom[k] =
 938             restp->atom    [k - hack->nr];
 939         restp->atomname[k] =
 940             restp->atomname[k - hack->nr];
 941         restp->cgnr[k] =
 942             restp->cgnr    [k - hack->nr];
 943     }
 944     /* now add them */
 945     for(k=0; k < hack->nr; k++)
 946     {
 947         /* set counter in atomname */
 948         if ( hack->nr > 1 )
 949         {
 950             buf[strlen(buf)-1] = Hnum[k];
 951         }
 952         snew( restp->atomname[at_start+1+k], 1);
 953         restp->atom     [at_start+1+k] = *hack->atom;
 954         *restp->atomname[at_start+1+k] = strdup(buf);
 955         if ( hack->cgnr != NOTSET )
 956         {
 957             restp->cgnr   [at_start+1+k] = hack->cgnr;
 958         }
 959         else
 960         {
 961             restp->cgnr   [at_start+1+k] = restp->cgnr[at_start];
 962         }
 963     }
 964 }
 965
 966 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
 967                                int nrtp, t_restp rtp[],
 968                                int nres, t_resinfo *resinfo,
 969                                int nterpairs,
 970                                t_hackblock **ntdb, t_hackblock **ctdb,
 971                                int *rn, int *rc)
 972 {
 973   int i, j, k, l;
 974   char *key;
 975   t_restp *res;
 976   char buf[STRLEN];
 977   const char *Hnum="123456";
 978   int tern,terc;
 979   gmx_bool bN,bC,bRM;
 980
 981   snew(*hb,nres);
 982   snew(*restp,nres);
 983   /* first the termini */
 984   for(i=0; i<nterpairs; i++) {
 985       if (rn[i] >= 0 && ntdb[i] != NULL) {
 986           copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
 987       }
 988       if (rc[i] >= 0 && ctdb[i] != NULL) {
 989           merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
 990       }
 991   }
 992
 993   /* then the whole rtp */
 994   for(i=0; i < nres; i++) {
 995     /* Here we allow a mismatch of one character when looking for the rtp entry.
 996      * For such a mismatch there should be only one mismatching name.
 997      * This is mainly useful for small molecules such as ions.
 998      * Note that this will usually not work for protein, DNA and RNA,
 999      * since there the residue names should be listed in residuetypes.dat
1000      * and an error will have been generated earlier in the process.
1001      */
1002     key = *resinfo[i].rtp;
1003     snew(resinfo[i].rtp,1);
1004     *resinfo[i].rtp = search_rtp(key,nrtp,rtp);
1005     res = get_restp(*resinfo[i].rtp,nrtp,rtp);
1006     copy_t_restp(res, &(*restp)[i]);
1007
1008     /* Check that we do not have different bonded types in one molecule */
1009     check_restp_types(&(*restp)[0],&(*restp)[i]);
1010
1011     tern = -1;
1012     for(j=0; j<nterpairs && tern==-1; j++) {
1013         if (i == rn[j]) {
1014             tern = j;
1015         }
1016     }
1017     terc = -1;
1018     for(j=0; j<nterpairs && terc == -1; j++) {
1019         if (i == rc[j]) {
1020             terc = j;
1021         }
1022     }
1023     bRM = merge_t_bondeds(res->rb, (*hb)[i].rb,tern>=0,terc>=0);
1024
1025     if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1026                 (terc >= 0 && ctdb[terc] == NULL))) {
1027         gmx_fatal(FARGS,"There is a dangling bond at at least one of the terminal ends and the force field does not provide terminal entries or files. Fix your terminal residues so that they match the residue database (.rtp) entries, or provide terminal database entries (.tdb).");
1028     }
1029     if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1030                 (terc >= 0 && ctdb[terc]->nhack == 0))) {
1031         gmx_fatal(FARGS,"There is a dangling bond at at least one of the terminal ends. Fix your coordinate file, add a new terminal database entry (.tdb), or select the proper existing terminal entry.");
1032     }
1033   }
1034
1035   /* now perform t_hack's on t_restp's,
1036      i.e. add's and deletes from termini database will be
1037      added to/removed from residue topology
1038      we'll do this on one big dirty loop, so it won't make easy reading! */
1039     for(i=0; i < nres; i++)
1040     {
1041         for(j=0; j < (*hb)[i].nhack; j++)
1042         {
1043             if ( (*hb)[i].hack[j].nr )
1044             {
1045                 /* find atom in restp */
1046                 for(l=0; l < (*restp)[i].natom; l++)
1047                     if ( ( (*hb)[i].hack[j].oname==NULL &&
1048                            strcmp((*hb)[i].hack[j].AI, *(*restp)[i].atomname[l])==0 ) ||
1049                          ( (*hb)[i].hack[j].oname!=NULL &&
1050                            strcmp((*hb)[i].hack[j].oname,*(*restp)[i].atomname[l])==0 ) )
1051                         break;
1052                 if (l == (*restp)[i].natom)
1053                 {
1054                     /* If we are doing an atom rename only, we don't need
1055                      * to generate a fatal error if the old name is not found
1056                      * in the rtp.
1057                      */
1058                     /* Deleting can happen also only on the input atoms,
1059                      * not necessarily always on the rtp entry.
1060                      */
1061                     if (!((*hb)[i].hack[j].oname != NULL &&
1062                           (*hb)[i].hack[j].nname != NULL) &&
1063                         !((*hb)[i].hack[j].oname != NULL &&
1064                           (*hb)[i].hack[j].nname == NULL))
1065                     {
1066                         gmx_fatal(FARGS,
1067                                   "atom %s not found in buiding block %d%s "
1068                                   "while combining tdb and rtp",
1069                                   (*hb)[i].hack[j].oname!=NULL ?
1070                                   (*hb)[i].hack[j].oname : (*hb)[i].hack[j].AI,
1071                                   i+1,*resinfo[i].rtp);
1072                     }
1073                 }
1074                 else
1075                 {
1076                     if ( (*hb)[i].hack[j].oname==NULL ) {
1077                         /* we're adding: */
1078                         add_atom_to_restp(&(*restp)[i],resinfo[i].nr,l,
1079                                           &(*hb)[i].hack[j]);
1080                     }
1081                     else
1082                     {
1083                         /* oname != NULL */
1084                         if ( (*hb)[i].hack[j].nname==NULL ) {
1085                             /* we're deleting */
1086                             if (debug)
1087                                 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1088                                         *(*restp)[i].atomname[l],
1089                                         i+1,(*restp)[i].resname);
1090                             /* shift the rest */
1091                             (*restp)[i].natom--;
1092                             for(k=l; k < (*restp)[i].natom; k++) {
1093                                 (*restp)[i].atom    [k] = (*restp)[i].atom    [k+1];
1094                                 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1095                                 (*restp)[i].cgnr    [k] = (*restp)[i].cgnr    [k+1];
1096                             }
1097                             /* give back space */
1098                             srenew((*restp)[i].atom,     (*restp)[i].natom);
1099                             srenew((*restp)[i].atomname, (*restp)[i].natom);
1100                             srenew((*restp)[i].cgnr,     (*restp)[i].natom);
1101                         } else { /* nname != NULL */
1102                             /* we're replacing */
1103                             if (debug)
1104                                 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1105                                         *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1106                                         i+1,(*restp)[i].resname);
1107                             snew( (*restp)[i].atomname[l], 1);
1108                             (*restp)[i].atom[l]      =       *(*hb)[i].hack[j].atom;
1109                             *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1110                             if ( (*hb)[i].hack[j].cgnr != NOTSET )
1111                                 (*restp)[i].cgnr   [l] = (*hb)[i].hack[j].cgnr;
1112                         }
1113                     }
1114                 }
1115             }
1116         }
1117     }
1118 }
1119
1120 static gmx_bool atomname_cmp_nr(const char *anm,t_hack *hack,int *nr)
1121 {
1122
1123     if (hack->nr == 1)
1124     {
1125         *nr = 0;
1126
1127         return (gmx_strcasecmp(anm,hack->nname) == 0);
1128     }
1129     else
1130     {
1131         if (isdigit(anm[strlen(anm)-1]))
1132         {
1133             *nr = anm[strlen(anm)-1] - '0';
1134         }
1135         else
1136         {
1137             *nr = 0;
1138         }
1139         if (*nr <= 0 || *nr > hack->nr)
1140         {
1141             return FALSE;
1142         }
1143         else
1144         {
1145             return (strlen(anm) == strlen(hack->nname) + 1 &&
1146                     gmx_strncasecmp(anm,hack->nname,strlen(hack->nname)) == 0);
1147         }
1148     }
1149 }
1150
1151 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba,rvec *x,int atind,
1152                                           t_restp *rptr,t_hackblock *hbr,
1153                                           gmx_bool bVerbose)
1154 {
1155     int  resnr;
1156     int  i,j,k;
1157     char *oldnm,*newnm;
1158     int  anmnr;
1159     char *start_at,buf[STRLEN];
1160     int  start_nr;
1161     gmx_bool bReplaceReplace,bFoundInAdd;
1162     gmx_bool bDeleted;
1163
1164     oldnm = *pdba->atomname[atind];
1165     resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1166
1167     bDeleted = FALSE;
1168     for(j=0; j<hbr->nhack; j++)
1169     {
1170         if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1171             gmx_strcasecmp(oldnm,hbr->hack[j].oname) == 0)
1172         {
1173             /* This is a replace entry. */
1174             /* Check if we are not replacing a replaced atom. */
1175             bReplaceReplace = FALSE;
1176             for(k=0; k<hbr->nhack; k++) {
1177                 if (k != j &&
1178                     hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1179                     gmx_strcasecmp(hbr->hack[k].nname,hbr->hack[j].oname) == 0)
1180                 {
1181                     /* The replace in hack[j] replaces an atom that
1182                      * was already replaced in hack[k], we do not want
1183                      * second or higher level replaces at this stage.
1184                      */
1185                     bReplaceReplace = TRUE;
1186                 }
1187             }
1188             if (bReplaceReplace)
1189             {
1190                 /* Skip this replace. */
1191                 continue;
1192             }
1193
1194             /* This atom still has the old name, rename it */
1195             newnm = hbr->hack[j].nname;
1196             for(k=0; k<rptr->natom; k++)
1197             {
1198                 if (gmx_strcasecmp(newnm,*rptr->atomname[k]) == 0)
1199                 {
1200                     break;
1201                 }
1202             }
1203             if (k == rptr->natom)
1204             {
1205                 /* The new name is not present in the rtp.
1206                  * We need to apply the replace also to the rtp entry.
1207                  */
1208
1209                 /* We need to find the add hack that can add this atom
1210                  * to find out after which atom it should be added.
1211                  */
1212                 bFoundInAdd = FALSE;
1213                 for(k=0; k<hbr->nhack; k++)
1214                 {
1215                     if (hbr->hack[k].oname == NULL &&
1216                         hbr->hack[k].nname != NULL &&
1217                         atomname_cmp_nr(newnm,&hbr->hack[k],&anmnr))
1218                     {
1219                         if (anmnr <= 1)
1220                         {
1221                             start_at = hbr->hack[k].a[0];
1222                         }
1223                         else
1224                         {
1225                             sprintf(buf,"%s%d",hbr->hack[k].nname,anmnr-1);
1226                             start_at = buf;
1227                         }
1228                         for(start_nr=0; start_nr<rptr->natom; start_nr++)
1229                         {
1230                             if (gmx_strcasecmp(start_at,(*rptr->atomname[start_nr])) == 0)
1231                             {
1232                                 break;
1233                             }
1234                         }
1235                         if (start_nr == rptr->natom)
1236                         {
1237                             gmx_fatal(FARGS,"Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1238                                       start_at,rptr->resname,newnm);
1239                         }
1240                         /* We can add the atom after atom start_nr */
1241                         add_atom_to_restp(rptr,resnr,start_nr,
1242                                           &hbr->hack[j]);
1243
1244                         bFoundInAdd = TRUE;
1245                     }
1246                 }
1247
1248                 if (!bFoundInAdd)
1249                 {
1250                     gmx_fatal(FARGS,"Could not find an 'add' entry for atom named '%s' corresponding to the 'replace' entry from atom name '%s' to '%s' for tdb or hdb database of residue type '%s'",
1251                               newnm,
1252                               hbr->hack[j].oname,hbr->hack[j].nname,
1253                               rptr->resname);
1254                 }
1255             }
1256
1257             if (bVerbose)
1258             {
1259                 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1260                        oldnm,rptr->resname,resnr,newnm);
1261             }
1262             /* Rename the atom in pdba */
1263             snew(pdba->atomname[atind],1);
1264             *pdba->atomname[atind] = strdup(newnm);
1265         }
1266         else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1267                  gmx_strcasecmp(oldnm,hbr->hack[j].oname) == 0)
1268         {
1269             /* This is a delete entry, check if this atom is present
1270              * in the rtp entry of this residue.
1271              */
1272             for(k=0; k<rptr->natom; k++)
1273             {
1274                 if (gmx_strcasecmp(oldnm,*rptr->atomname[k]) == 0)
1275                 {
1276                     break;
1277                 }
1278             }
1279             if (k == rptr->natom)
1280             {
1281                 /* This atom is not present in the rtp entry,
1282                  * delete is now from the input pdba.
1283                  */
1284                 if (bVerbose)
1285                 {
1286                     printf("Deleting atom '%s' in residue '%s' %d\n",
1287                            oldnm,rptr->resname,resnr);
1288                 }
1289                 /* We should free the atom name,
1290                  * but it might be used multiple times in the symtab.
1291                  * sfree(pdba->atomname[atind]);
1292                  */
1293                 for(k=atind+1; k<pdba->nr; k++)
1294                 {
1295                     pdba->atom[k-1]     = pdba->atom[k];
1296                     pdba->atomname[k-1] = pdba->atomname[k];
1297                     copy_rvec(x[k],x[k-1]);
1298                 }
1299                 pdba->nr--;
1300                 bDeleted = TRUE;
1301             }
1302         }
1303     }
1304
1305     return bDeleted;
1306 }
1307
1308 void match_atomnames_with_rtp(t_restp restp[],t_hackblock hb[],
1309                               t_atoms *pdba,rvec *x,
1310                               gmx_bool bVerbose)
1311 {
1312     int  i,j,k;
1313     char *oldnm,*newnm;
1314     int  resnr;
1315     t_restp *rptr;
1316     t_hackblock *hbr;
1317     int  anmnr;
1318     char *start_at,buf[STRLEN];
1319     int  start_nr;
1320     gmx_bool bFoundInAdd;
1321
1322     for(i=0; i<pdba->nr; i++)
1323     {
1324         oldnm = *pdba->atomname[i];
1325         resnr = pdba->resinfo[pdba->atom[i].resind].nr;
1326         rptr  = &restp[pdba->atom[i].resind];
1327         for(j=0; (j<rptr->natom); j++)
1328         {
1329             if (gmx_strcasecmp(oldnm,*(rptr->atomname[j])) == 0)
1330             {
1331                 break;
1332             }
1333         }
1334         if (j == rptr->natom)
1335         {
1336             /* Not found yet, check if we have to rename this atom */
1337             if (match_atomnames_with_rtp_atom(pdba,x,i,
1338                                               rptr,&(hb[pdba->atom[i].resind]),
1339                                               bVerbose))
1340             {
1341                 /* We deleted this atom, decrease the atom counter by 1. */
1342                 i--;
1343             }
1344         }
1345     }
1346 }
1347
1348 #define NUM_CMAP_ATOMS 5
1349 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t rt)
1350 {
1351     int residx,i,j,k;
1352     const char *ptr;
1353     t_resinfo *resinfo = atoms->resinfo;
1354     int nres = atoms->nres;
1355     gmx_bool bAddCMAP;
1356     atom_id cmap_atomid[NUM_CMAP_ATOMS];
1357     int cmap_chainnum=-1, this_residue_index;
1358
1359         if (debug)
1360                 ptr = "cmap";
1361         else
1362                 ptr = "check";
1363
1364         fprintf(stderr,"Making cmap torsions...");
1365         i=0;
1366         /* End loop at nres-1, since the very last residue does not have a +N atom, and
1367          * therefore we get a valgrind invalid 4 byte read error with atom am */
1368         for(residx=0; residx<nres-1; residx++)
1369         {
1370                 /* Add CMAP terms from the list of CMAP interactions */
1371                 for(j=0;j<restp[residx].rb[ebtsCMAP].nb; j++)
1372                 {
1373             bAddCMAP = TRUE;
1374             /* Loop over atoms in a candidate CMAP interaction and
1375              * check that they exist, are from the same chain and are
1376              * from residues labelled as protein. */
1377             for(k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1378             {
1379                 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
1380                                              i,atoms,ptr,TRUE);
1381                 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1382                 if (!bAddCMAP)
1383                 {
1384                     /* This break is necessary, because cmap_atomid[k]
1385                      * == NO_ATID cannot be safely used as an index
1386                      * into the atom array. */
1387                     break;
1388                 }
1389                 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1390                 if (0 == k)
1391                 {
1392                     cmap_chainnum = resinfo[this_residue_index].chainnum;
1393                 }
1394                 else
1395                 {
1396                     /* Does the residue for this atom have the same
1397                      * chain number as the residues for previous
1398                      * atoms? */
1399                     bAddCMAP = bAddCMAP &&
1400                         cmap_chainnum == resinfo[this_residue_index].chainnum;
1401                 }
1402                 bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt,*(resinfo[this_residue_index].name));
1403             }
1404
1405             if(bAddCMAP)
1406             {
1407                 add_cmap_param(psb,cmap_atomid[0],cmap_atomid[1],cmap_atomid[2],cmap_atomid[3],cmap_atomid[4],restp[residx].rb[ebtsCMAP].b[j].s);
1408                         }
1409                 }
1410
1411                 if(residx<nres-1)
1412                 {
1413                         while(atoms->atom[i].resind<residx+1)
1414                         {
1415                                 i++;
1416                         }
1417                 }
1418         }
1419
1420         /* Start the next residue */
1421 }
1422
1423 static void
1424 scrub_charge_groups(int *cgnr, int natoms)
1425 {
1426         int i;
1427
1428         for(i=0;i<natoms;i++)
1429         {
1430                 cgnr[i]=i+1;
1431         }
1432 }
1433
1434
1435 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1436              t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1437              int nrtp, t_restp rtp[],
1438              t_restp *restp, t_hackblock *hb,
1439              int nterpairs,t_hackblock **ntdb, t_hackblock **ctdb,
1440              gmx_bool bAllowMissing,
1441              gmx_bool bVsites, gmx_bool bVsiteAromatics,
1442              const char *ff, const char *ffdir,
1443              real mHmult,
1444              int nssbonds, t_ssbond *ssbonds,
1445              real long_bond_dist, real short_bond_dist,
1446              gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1447              gmx_bool bRenumRes,gmx_bool bRTPresname)
1448 {
1449     /*
1450   t_hackblock *hb;
1451   t_restp  *restp;
1452     */
1453   t_params plist[F_NRE];
1454   t_excls  *excls;
1455   t_nextnb nnb;
1456   int      *cgnr;
1457   int      *vsite_type;
1458   int      i,nmissat;
1459   int      bts[ebtsNR];
1460   gmx_residuetype_t rt;
1461
1462   init_plist(plist);
1463   gmx_residuetype_init(&rt);
1464
1465   if (debug) {
1466     print_resall(debug, atoms->nres, restp, atype);
1467     dump_hb(debug, atoms->nres, hb);
1468   }
1469
1470   /* Make bonds */
1471   at2bonds(&(plist[F_BONDS]), hb,
1472            atoms, *x,
1473            long_bond_dist, short_bond_dist, bAllowMissing);
1474
1475   /* specbonds: disulphide bonds & heme-his */
1476   do_ssbonds(&(plist[F_BONDS]),
1477              atoms, nssbonds, ssbonds,
1478              bAllowMissing);
1479
1480   nmissat = name2type(atoms, &cgnr, atype, restp, rt);
1481   if (nmissat) {
1482     if (bAllowMissing)
1483       fprintf(stderr,"There were %d missing atoms in molecule %s\n",
1484               nmissat,molname);
1485     else
1486       gmx_fatal(FARGS,"There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1487                   nmissat,molname);
1488   }
1489
1490   /* Cleanup bonds (sort and rm doubles) */
1491   clean_bonds(&(plist[F_BONDS]));
1492
1493   snew(vsite_type,atoms->nr);
1494   for(i=0; i<atoms->nr; i++)
1495     vsite_type[i]=NOTSET;
1496   if (bVsites) {
1497     /* determine which atoms will be vsites and add dummy masses
1498        also renumber atom numbers in plist[0..F_NRE]! */
1499     do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1500               &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1501   }
1502
1503   /* Make Angles and Dihedrals */
1504   fprintf(stderr,"Generating angles, dihedrals and pairs...\n");
1505   snew(excls,atoms->nr);
1506   init_nnb(&nnb,atoms->nr,4);
1507   gen_nnb(&nnb,plist);
1508   print_nnb(&nnb,"NNB");
1509   gen_pad(&nnb,atoms,restp,plist,excls,hb,bAllowMissing);
1510   done_nnb(&nnb);
1511
1512     /* Make CMAP */
1513     if (TRUE == bCmap)
1514     {
1515         gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1516         if (plist[F_CMAP].nr > 0)
1517         {
1518             fprintf(stderr, "There are %4d cmap torsion pairs\n",
1519                     plist[F_CMAP].nr);
1520         }
1521     }
1522
1523   /* set mass of all remaining hydrogen atoms */
1524   if (mHmult != 1.0)
1525     do_h_mass(&(plist[F_BONDS]),vsite_type,atoms,mHmult,bDeuterate);
1526   sfree(vsite_type);
1527
1528   /* Cleanup bonds (sort and rm doubles) */
1529   /* clean_bonds(&(plist[F_BONDS]));*/
1530
1531   fprintf(stderr,
1532           "There are %4d dihedrals, %4d impropers, %4d angles\n"
1533           "          %4d pairs,     %4d bonds and  %4d virtual sites\n",
1534           plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1535           plist[F_LJ14].nr, plist[F_BONDS].nr,
1536           plist[F_VSITE2].nr +
1537           plist[F_VSITE3].nr +
1538           plist[F_VSITE3FD].nr +
1539           plist[F_VSITE3FAD].nr +
1540           plist[F_VSITE3OUT].nr +
1541       plist[F_VSITE4FD].nr +
1542       plist[F_VSITE4FDN].nr );
1543
1544   print_sums(atoms, FALSE);
1545
1546   if (FALSE == bChargeGroups)
1547   {
1548           scrub_charge_groups(cgnr, atoms->nr);
1549   }
1550
1551     if (bRenumRes)
1552     {
1553         for(i=0; i<atoms->nres; i++)
1554         {
1555             atoms->resinfo[i].nr = i + 1;
1556             atoms->resinfo[i].ic = ' ';
1557         }
1558     }
1559
1560   if (top_file) {
1561     fprintf(stderr,"Writing topology\n");
1562     /* We can copy the bonded types from the first restp,
1563      * since the types have to be identical for all residues in one molecule.
1564      */
1565     for(i=0; i<ebtsNR; i++) {
1566         bts[i] = restp[0].rb[i].type;
1567     }
1568     write_top(top_file, posre_fn, molname,
1569               atoms, bRTPresname,
1570               bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1571   }
1572
1573   /* cleaning up */
1574   free_t_hackblock(atoms->nres, &hb);
1575   free_t_restp(atoms->nres, &restp);
1576   gmx_residuetype_destroy(rt);
1577
1578   /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1579   sfree(cgnr);
1580   for (i=0; i<F_NRE; i++)
1581     sfree(plist[i].param);
1582   for (i=0; i<atoms->nr; i++)
1583     sfree(excls[i].e);
1584   sfree(excls);
1585 }