src/gromacs/gmxpreprocess/pdb2top.c

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