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