src/gromacs/gmxpreprocess/pdb2top.c

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