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