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