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