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