src/gromacs/gmxpreprocess/gen_vsite.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "string2.h"
  40 #include <stdio.h>
  41 #include <math.h>
  42 #include <string.h>
  43 #include "gen_vsite.h"
  44 #include "smalloc.h"
  45 #include "resall.h"
  46 #include "add_par.h"
  47 #include "vec.h"
  48 #include "toputil.h"
  49 #include "physics.h"
  50 #include "index.h"
  51 #include "names.h"
  52 #include "futil.h"
  53 #include "gpp_atomtype.h"
  54 #include "fflibutil.h"
  55 #include "macros.h"
  56
  57 #define MAXNAME 32
  58 #define OPENDIR     '[' /* starting sign for directive          */
  59 #define CLOSEDIR    ']' /* ending sign for directive            */
  60
  61 typedef struct {
  62     char       atomtype[MAXNAME];  /* Type for the XH3/XH2 atom */
  63     gmx_bool   isplanar;           /* If true, the atomtype above and the three connected
  64                                     * ones are in a planar geometry. The two next entries
  65                                     * are undefined in that case
  66                                     */
  67     int    nhydrogens;             /* number of connected hydrogens */
  68     char   nextheavytype[MAXNAME]; /* Type for the heavy atom bonded to XH2/XH3 */
  69     char   dummymass[MAXNAME];     /* The type of MNH* or MCH3* dummy mass to use */
  70 } t_vsiteconf;
  71
  72
  73 /* Structure to represent average bond and angles values in vsite aromatic
  74  * residues. Note that these are NOT necessarily the bonds and angles from the
  75  * forcefield; many forcefields (like Amber, OPLS) have some inherent strain in
  76  * 5-rings (i.e. the sum of angles is !=540, but impropers keep it planar)
  77  */
  78 typedef struct {
  79     char resname[MAXNAME];
  80     int  nbonds;
  81     int  nangles;
  82     struct vsitetop_bond {
  83         char   atom1[MAXNAME];
  84         char   atom2[MAXNAME];
  85         float  value;
  86     } *bond; /* list of bonds */
  87     struct vsitetop_angle {
  88         char   atom1[MAXNAME];
  89         char   atom2[MAXNAME];
  90         char   atom3[MAXNAME];
  91         float  value;
  92     } *angle; /* list of angles */
  93 } t_vsitetop;
  94
  95
  96 enum {
  97     DDB_CH3, DDB_NH3, DDB_NH2, DDB_PHE, DDB_TYR,
  98     DDB_TRP, DDB_HISA, DDB_HISB, DDB_HISH, DDB_DIR_NR
  99 };
 100
 101 typedef char t_dirname[STRLEN];
 102
 103 static const t_dirname ddb_dirnames[DDB_DIR_NR] = {
 104     "CH3",
 105     "NH3",
 106     "NH2",
 107     "PHE",
 108     "TYR",
 109     "TRP",
 110     "HISA",
 111     "HISB",
 112     "HISH"
 113 };
 114
 115 static int ddb_name2dir(char *name)
 116 {
 117     /* Translate a directive name to the number of the directive.
 118      * HID/HIE/HIP names are translated to the ones we use in Gromacs.
 119      */
 120
 121     int i, index;
 122
 123     index = -1;
 124
 125     for (i = 0; i < DDB_DIR_NR && index < 0; i++)
 126     {
 127         if (!gmx_strcasecmp(name, ddb_dirnames[i]))
 128         {
 129             index = i;
 130         }
 131     }
 132
 133     return index;
 134 }
 135
 136
 137 static void read_vsite_database(const char *ddbname,
 138                                 t_vsiteconf **pvsiteconflist, int *nvsiteconf,
 139                                 t_vsitetop **pvsitetoplist, int *nvsitetop)
 140 {
 141     /* This routine is a quick hack to fix the problem with hardcoded atomtypes
 142      * and aromatic vsite parameters by reading them from a ff???.vsd file.
 143      *
 144      * The file can contain sections [ NH3 ], [ CH3 ], [ NH2 ], and ring residue names.
 145      * For the NH3 and CH3 section each line has three fields. The first is the atomtype
 146      * (nb: not bonded type) of the N/C atom to be replaced, the second field is
 147      * the type of the next heavy atom it is bonded to, and the third field the type
 148      * of dummy mass that will be used for this group.
 149      *
 150      * If the NH2 group planar (sp2 N) a different vsite construct is used, so in this
 151      * case the second field should just be the word planar.
 152      */
 153
 154     FILE        *ddb;
 155     char         dirstr[STRLEN];
 156     char         pline[STRLEN];
 157     int          i, j, n, k, nvsite, ntop, curdir, prevdir;
 158     t_vsiteconf *vsiteconflist;
 159     t_vsitetop  *vsitetoplist;
 160     char        *ch;
 161     char         s1[MAXNAME], s2[MAXNAME], s3[MAXNAME], s4[MAXNAME];
 162
 163     ddb = libopen(ddbname);
 164
 165     nvsite        = *nvsiteconf;
 166     vsiteconflist = *pvsiteconflist;
 167     ntop          = *nvsitetop;
 168     vsitetoplist  = *pvsitetoplist;
 169
 170     curdir = -1;
 171
 172     snew(vsiteconflist, 1);
 173     snew(vsitetoplist, 1);
 174
 175     while (fgets2(pline, STRLEN-2, ddb) != NULL)
 176     {
 177         strip_comment(pline);
 178         trim(pline);
 179         if (strlen(pline) > 0)
 180         {
 181             if (pline[0] == OPENDIR)
 182             {
 183                 strncpy(dirstr, pline+1, STRLEN-2);
 184                 if ((ch = strchr (dirstr, CLOSEDIR)) != NULL)
 185                 {
 186                     (*ch) = 0;
 187                 }
 188                 trim (dirstr);
 189
 190                 if (!gmx_strcasecmp(dirstr, "HID") ||
 191                     !gmx_strcasecmp(dirstr, "HISD"))
 192                 {
 193                     sprintf(dirstr, "HISA");
 194                 }
 195                 else if (!gmx_strcasecmp(dirstr, "HIE") ||
 196                          !gmx_strcasecmp(dirstr, "HISE"))
 197                 {
 198                     sprintf(dirstr, "HISB");
 199                 }
 200                 else if (!gmx_strcasecmp(dirstr, "HIP"))
 201                 {
 202                     sprintf(dirstr, "HISH");
 203                 }
 204
 205                 curdir = ddb_name2dir(dirstr);
 206                 if (curdir < 0)
 207                 {
 208                     gmx_fatal(FARGS, "Invalid directive %s in vsite database %s",
 209                               dirstr, ddbname);
 210                 }
 211             }
 212             else
 213             {
 214                 switch (curdir)
 215                 {
 216                     case -1:
 217                         gmx_fatal(FARGS, "First entry in vsite database must be a directive.\n");
 218                         break;
 219                     case DDB_CH3:
 220                     case DDB_NH3:
 221                     case DDB_NH2:
 222                         n = sscanf(pline, "%s%s%s", s1, s2, s3);
 223                         if (n < 3 && !gmx_strcasecmp(s2, "planar"))
 224                         {
 225                             srenew(vsiteconflist, nvsite+1);
 226                             strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1);
 227                             vsiteconflist[nvsite].isplanar         = TRUE;
 228                             vsiteconflist[nvsite].nextheavytype[0] = 0;
 229                             vsiteconflist[nvsite].dummymass[0]     = 0;
 230                             vsiteconflist[nvsite].nhydrogens       = 2;
 231                             nvsite++;
 232                         }
 233                         else if (n == 3)
 234                         {
 235                             srenew(vsiteconflist, (nvsite+1));
 236                             strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1);
 237                             vsiteconflist[nvsite].isplanar = FALSE;
 238                             strncpy(vsiteconflist[nvsite].nextheavytype, s2, MAXNAME-1);
 239                             strncpy(vsiteconflist[nvsite].dummymass, s3, MAXNAME-1);
 240                             if (curdir == DDB_NH2)
 241                             {
 242                                 vsiteconflist[nvsite].nhydrogens = 2;
 243                             }
 244                             else
 245                             {
 246                                 vsiteconflist[nvsite].nhydrogens = 3;
 247                             }
 248                             nvsite++;
 249                         }
 250                         else
 251                         {
 252                             gmx_fatal(FARGS, "Not enough directives in vsite database line: %s\n", pline);
 253                         }
 254                         break;
 255                     case DDB_PHE:
 256                     case DDB_TYR:
 257                     case DDB_TRP:
 258                     case DDB_HISA:
 259                     case DDB_HISB:
 260                     case DDB_HISH:
 261                         i = 0;
 262                         while ((i < ntop) && gmx_strcasecmp(dirstr, vsitetoplist[i].resname))
 263                         {
 264                             i++;
 265                         }
 266                         /* Allocate a new topology entry if this is a new residue */
 267                         if (i == ntop)
 268                         {
 269                             srenew(vsitetoplist, ntop+1);
 270                             ntop++; /* i still points to current vsite topology entry */
 271                             strncpy(vsitetoplist[i].resname, dirstr, MAXNAME-1);
 272                             vsitetoplist[i].nbonds = vsitetoplist[i].nangles = 0;
 273                             snew(vsitetoplist[i].bond, 1);
 274                             snew(vsitetoplist[i].angle, 1);
 275                         }
 276                         n = sscanf(pline, "%s%s%s%s", s1, s2, s3, s4);
 277                         if (n == 3)
 278                         {
 279                             /* bond */
 280                             k = vsitetoplist[i].nbonds++;
 281                             srenew(vsitetoplist[i].bond, k+1);
 282                             strncpy(vsitetoplist[i].bond[k].atom1, s1, MAXNAME-1);
 283                             strncpy(vsitetoplist[i].bond[k].atom2, s2, MAXNAME-1);
 284                             vsitetoplist[i].bond[k].value = strtod(s3, NULL);
 285                         }
 286                         else if (n == 4)
 287                         {
 288                             /* angle */
 289                             k = vsitetoplist[i].nangles++;
 290                             srenew(vsitetoplist[i].angle, k+1);
 291                             strncpy(vsitetoplist[i].angle[k].atom1, s1, MAXNAME-1);
 292                             strncpy(vsitetoplist[i].angle[k].atom2, s2, MAXNAME-1);
 293                             strncpy(vsitetoplist[i].angle[k].atom3, s3, MAXNAME-1);
 294                             vsitetoplist[i].angle[k].value = strtod(s4, NULL);
 295                         }
 296                         else
 297                         {
 298                             gmx_fatal(FARGS, "Need 3 or 4 values to specify bond/angle values in %s: %s\n", ddbname, pline);
 299                         }
 300                         break;
 301                     default:
 302                         gmx_fatal(FARGS, "Didnt find a case for directive %s in read_vsite_database\n", dirstr);
 303                 }
 304             }
 305         }
 306     }
 307
 308     *pvsiteconflist = vsiteconflist;
 309     *pvsitetoplist  = vsitetoplist;
 310     *nvsiteconf     = nvsite;
 311     *nvsitetop      = ntop;
 312
 313     ffclose(ddb);
 314 }
 315
 316 static int nitrogen_is_planar(t_vsiteconf vsiteconflist[], int nvsiteconf, char atomtype[])
 317 {
 318     /* Return 1 if atomtype exists in database list and is planar, 0 if not,
 319      * and -1 if not found.
 320      */
 321     int      i, res;
 322     gmx_bool found = FALSE;
 323     for (i = 0; i < nvsiteconf && !found; i++)
 324     {
 325         found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atomtype) && (vsiteconflist[i].nhydrogens == 2));
 326     }
 327     if (found)
 328     {
 329         res = (vsiteconflist[i-1].isplanar == TRUE);
 330     }
 331     else
 332     {
 333         res = -1;
 334     }
 335
 336     return res;
 337 }
 338
 339 static char *get_dummymass_name(t_vsiteconf vsiteconflist[], int nvsiteconf, char atom[], char nextheavy[])
 340 {
 341     /* Return the dummy mass name if found, or NULL if not set in ddb database */
 342     int      i;
 343     gmx_bool found = FALSE;
 344     for (i = 0; i < nvsiteconf && !found; i++)
 345     {
 346         found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atom) &&
 347                  !gmx_strcasecmp(vsiteconflist[i].nextheavytype, nextheavy));
 348     }
 349     if (found)
 350     {
 351         return vsiteconflist[i-1].dummymass;
 352     }
 353     else
 354     {
 355         return NULL;
 356     }
 357 }
 358
 359
 360
 361 static real get_ddb_bond(t_vsitetop *vsitetop, int nvsitetop,
 362                          const char res[],
 363                          const char atom1[], const char atom2[])
 364 {
 365     int i, j;
 366
 367     i = 0;
 368     while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
 369     {
 370         i++;
 371     }
 372     if (i == nvsitetop)
 373     {
 374         gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res);
 375     }
 376     j = 0;
 377     while (j < vsitetop[i].nbonds &&
 378            ( strcmp(atom1, vsitetop[i].bond[j].atom1) || strcmp(atom2, vsitetop[i].bond[j].atom2)) &&
 379            ( strcmp(atom2, vsitetop[i].bond[j].atom1) || strcmp(atom1, vsitetop[i].bond[j].atom2)))
 380     {
 381         j++;
 382     }
 383     if (j == vsitetop[i].nbonds)
 384     {
 385         gmx_fatal(FARGS, "Couldnt find bond %s-%s for residue %s in vsite database.\n", atom1, atom2, res);
 386     }
 387
 388     return vsitetop[i].bond[j].value;
 389 }
 390
 391
 392 static real get_ddb_angle(t_vsitetop *vsitetop, int nvsitetop,
 393                           const char res[], const char atom1[],
 394                           const char atom2[], const char atom3[])
 395 {
 396     int i, j;
 397
 398     i = 0;
 399     while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
 400     {
 401         i++;
 402     }
 403     if (i == nvsitetop)
 404     {
 405         gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res);
 406     }
 407     j = 0;
 408     while (j < vsitetop[i].nangles &&
 409            ( strcmp(atom1, vsitetop[i].angle[j].atom1) ||
 410              strcmp(atom2, vsitetop[i].angle[j].atom2) ||
 411              strcmp(atom3, vsitetop[i].angle[j].atom3)) &&
 412            ( strcmp(atom3, vsitetop[i].angle[j].atom1) ||
 413              strcmp(atom2, vsitetop[i].angle[j].atom2) ||
 414              strcmp(atom1, vsitetop[i].angle[j].atom3)))
 415     {
 416         j++;
 417     }
 418     if (j == vsitetop[i].nangles)
 419     {
 420         gmx_fatal(FARGS, "Couldnt find angle %s-%s-%s for residue %s in vsite database.\n", atom1, atom2, atom3, res);
 421     }
 422
 423     return vsitetop[i].angle[j].value;
 424 }
 425
 426
 427 static void count_bonds(int atom, t_params *psb, char ***atomname,
 428                         int *nrbonds, int *nrHatoms, int Hatoms[], int *Heavy,
 429                         int *nrheavies, int heavies[])
 430 {
 431     int i, heavy, other, nrb, nrH, nrhv;
 432
 433     /* find heavy atom bound to this hydrogen */
 434     heavy = NOTSET;
 435     for (i = 0; (i < psb->nr) && (heavy == NOTSET); i++)
 436     {
 437         if (psb->param[i].AI == atom)
 438         {
 439             heavy = psb->param[i].AJ;
 440         }
 441         else if (psb->param[i].AJ == atom)
 442         {
 443             heavy = psb->param[i].AI;
 444         }
 445     }
 446     if (heavy == NOTSET)
 447     {
 448         gmx_fatal(FARGS, "unbound hydrogen atom %d", atom+1);
 449     }
 450     /* find all atoms bound to heavy atom */
 451     other = NOTSET;
 452     nrb   = 0;
 453     nrH   = 0;
 454     nrhv  = 0;
 455     for (i = 0; i < psb->nr; i++)
 456     {
 457         if (psb->param[i].AI == heavy)
 458         {
 459             other = psb->param[i].AJ;
 460         }
 461         else if (psb->param[i].AJ == heavy)
 462         {
 463             other = psb->param[i].AI;
 464         }
 465         if (other != NOTSET)
 466         {
 467             nrb++;
 468             if (is_hydrogen(*(atomname[other])))
 469             {
 470                 Hatoms[nrH] = other;
 471                 nrH++;
 472             }
 473             else
 474             {
 475                 heavies[nrhv] = other;
 476                 nrhv++;
 477             }
 478             other = NOTSET;
 479         }
 480     }
 481     *Heavy     = heavy;
 482     *nrbonds   = nrb;
 483     *nrHatoms  = nrH;
 484     *nrheavies = nrhv;
 485 }
 486
 487 static void print_bonds(FILE *fp, int o2n[],
 488                         int nrHatoms, int Hatoms[], int Heavy,
 489                         int nrheavies, int heavies[])
 490 {
 491     int i;
 492
 493     fprintf(fp, "Found: %d Hatoms: ", nrHatoms);
 494     for (i = 0; i < nrHatoms; i++)
 495     {
 496         fprintf(fp, " %d", o2n[Hatoms[i]]+1);
 497     }
 498     fprintf(fp, "; %d Heavy atoms: %d", nrheavies+1, o2n[Heavy]+1);
 499     for (i = 0; i < nrheavies; i++)
 500     {
 501         fprintf(fp, " %d", o2n[heavies[i]]+1);
 502     }
 503     fprintf(fp, "\n");
 504 }
 505
 506 static int get_atype(int atom, t_atoms *at, int nrtp, t_restp rtp[],
 507                      gmx_residuetype_t rt)
 508 {
 509     int      type;
 510     gmx_bool bNterm;
 511     int      j;
 512     t_restp *rtpp;
 513
 514     if (at->atom[atom].m)
 515     {
 516         type = at->atom[atom].type;
 517     }
 518     else
 519     {
 520         /* get type from rtp */
 521         rtpp   = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
 522         bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
 523             (at->atom[atom].resind == 0);
 524         j    = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
 525         type = rtpp->atom[j].type;
 526     }
 527     return type;
 528 }
 529
 530 static int vsite_nm2type(const char *name, gpp_atomtype_t atype)
 531 {
 532     int tp;
 533
 534     tp = get_atomtype_type(name, atype);
 535     if (tp == NOTSET)
 536     {
 537         gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database",
 538                   name);
 539     }
 540
 541     return tp;
 542 }
 543
 544 static real get_amass(int atom, t_atoms *at, int nrtp, t_restp rtp[],
 545                       gmx_residuetype_t rt)
 546 {
 547     real     mass;
 548     gmx_bool bNterm;
 549     int      j;
 550     t_restp *rtpp;
 551
 552     if (at->atom[atom].m)
 553     {
 554         mass = at->atom[atom].m;
 555     }
 556     else
 557     {
 558         /* get mass from rtp */
 559         rtpp   = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
 560         bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
 561             (at->atom[atom].resind == 0);
 562         j    = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
 563         mass = rtpp->atom[j].m;
 564     }
 565     return mass;
 566 }
 567
 568 static void my_add_param(t_params *plist, int ai, int aj, real b)
 569 {
 570     static real c[MAXFORCEPARAM] =
 571     { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
 572
 573     c[0] = b;
 574     add_param(plist, ai, aj, c, NULL);
 575 }
 576
 577 static void add_vsites(t_params plist[], int vsite_type[],
 578                        int Heavy, int nrHatoms, int Hatoms[],
 579                        int nrheavies, int heavies[])
 580 {
 581     int      i, j, ftype, other, moreheavy, bb;
 582     gmx_bool bSwapParity;
 583
 584     for (i = 0; i < nrHatoms; i++)
 585     {
 586         ftype = vsite_type[Hatoms[i]];
 587         /* Errors in setting the vsite_type should really be caugth earlier,
 588          * because here it's not possible to print any useful error message.
 589          * But it's still better to print a message than to segfault.
 590          */
 591         if (ftype == NOTSET)
 592         {
 593             gmx_incons("Undetected error in setting up virtual sites");
 594         }
 595         bSwapParity           = (ftype < 0);
 596         vsite_type[Hatoms[i]] = ftype = abs(ftype);
 597         if (ftype == F_BONDS)
 598         {
 599             if ( (nrheavies != 1) && (nrHatoms != 1) )
 600             {
 601                 gmx_fatal(FARGS, "cannot make constraint in add_vsites for %d heavy "
 602                           "atoms and %d hydrogen atoms", nrheavies, nrHatoms);
 603             }
 604             my_add_param(&(plist[F_CONSTRNC]), Hatoms[i], heavies[0], NOTSET);
 605         }
 606         else
 607         {
 608             switch (ftype)
 609             {
 610                 case F_VSITE3:
 611                 case F_VSITE3FD:
 612                 case F_VSITE3OUT:
 613                     if (nrheavies < 2)
 614                     {
 615                         gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 3)",
 616                                   nrheavies+1,
 617                                   interaction_function[vsite_type[Hatoms[i]]].name);
 618                     }
 619                     add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1],
 620                                      bSwapParity);
 621                     break;
 622                 case F_VSITE3FAD:
 623                 {
 624                     if (nrheavies > 1)
 625                     {
 626                         moreheavy = heavies[1];
 627                     }
 628                     else
 629                     {
 630                         /* find more heavy atoms */
 631                         other = moreheavy = NOTSET;
 632                         for (j = 0; (j < plist[F_BONDS].nr) && (moreheavy == NOTSET); j++)
 633                         {
 634                             if (plist[F_BONDS].param[j].AI == heavies[0])
 635                             {
 636                                 other = plist[F_BONDS].param[j].AJ;
 637                             }
 638                             else if (plist[F_BONDS].param[j].AJ == heavies[0])
 639                             {
 640                                 other = plist[F_BONDS].param[j].AI;
 641                             }
 642                             if ( (other != NOTSET) && (other != Heavy) )
 643                             {
 644                                 moreheavy = other;
 645                             }
 646                         }
 647                         if (moreheavy == NOTSET)
 648                         {
 649                             gmx_fatal(FARGS, "Unbound molecule part %d-%d", Heavy+1, Hatoms[0]+1);
 650                         }
 651                     }
 652                     add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy,
 653                                      bSwapParity);
 654                     break;
 655                 }
 656                 case F_VSITE4FD:
 657                 case F_VSITE4FDN:
 658                     if (nrheavies < 3)
 659                     {
 660                         gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 4)",
 661                                   nrheavies+1,
 662                                   interaction_function[vsite_type[Hatoms[i]]].name);
 663                     }
 664                     add_vsite4_atoms(&plist[ftype],
 665                                      Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
 666                     break;
 667
 668                 default:
 669                     gmx_fatal(FARGS, "can't use add_vsites for interaction function %s",
 670                               interaction_function[vsite_type[Hatoms[i]]].name);
 671             } /* switch ftype */
 672         }     /* else */
 673     }         /* for i */
 674 }
 675
 676 #define ANGLE_6RING (DEG2RAD*120)
 677
 678 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
 679 /* get a^2 when a, b and alpha are given: */
 680 #define cosrule(b, c, alpha) ( sqr(b) + sqr(c) - 2*b*c*cos(alpha) )
 681 /* get cos(alpha) when a, b and c are given: */
 682 #define acosrule(a, b, c) ( (sqr(b)+sqr(c)-sqr(a))/(2*b*c) )
 683
 684 static int gen_vsites_6ring(t_atoms *at, int *vsite_type[], t_params plist[],
 685                             int nrfound, int *ats, real bond_cc, real bond_ch,
 686                             real xcom, gmx_bool bDoZ)
 687 {
 688     /* these MUST correspond to the atnms array in do_vsite_aromatics! */
 689     enum {
 690         atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
 691         atCZ, atHZ, atNR
 692     };
 693
 694     int  i, nvsite;
 695     real a, b, dCGCE, tmp1, tmp2, mtot, mG, mrest;
 696     real xCG, yCG, xCE1, yCE1, xCE2, yCE2;
 697     /* CG, CE1 and CE2 stay and each get a part of the total mass,
 698      * so the c-o-m stays the same.
 699      */
 700
 701     if (bDoZ)
 702     {
 703         if (atNR != nrfound)
 704         {
 705             gmx_incons("Generating vsites on 6-rings");
 706         }
 707     }
 708
 709     /* constraints between CG, CE1 and CE2: */
 710     dCGCE = sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
 711     my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE);
 712     my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE2], dCGCE);
 713     my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atCE2], dCGCE);
 714
 715     /* rest will be vsite3 */
 716     mtot   = 0;
 717     nvsite = 0;
 718     for (i = 0; i <  (bDoZ ? atNR : atHZ); i++)
 719     {
 720         mtot += at->atom[ats[i]].m;
 721         if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ) ) )
 722         {
 723             at->atom[ats[i]].m    = at->atom[ats[i]].mB = 0;
 724             (*vsite_type)[ats[i]] = F_VSITE3;
 725             nvsite++;
 726         }
 727     }
 728     /* Distribute mass so center-of-mass stays the same.
 729      * The center-of-mass in the call is defined with x=0 at
 730      * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
 731      */
 732     xCG  = -bond_cc+bond_cc*cos(ANGLE_6RING);
 733     yCG  = 0;
 734     xCE1 = 0;
 735     yCE1 = bond_cc*sin(0.5*ANGLE_6RING);
 736     xCE2 = 0;
 737     yCE2 = -bond_cc*sin(0.5*ANGLE_6RING);
 738
 739     mG                             = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom*mtot/xCG;
 740     mrest                          = mtot-mG;
 741     at->atom[ats[atCE1]].m         = at->atom[ats[atCE1]].mB =
 742             at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mrest / 2;
 743
 744     /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
 745     tmp1  = dCGCE*sin(ANGLE_6RING*0.5);
 746     tmp2  = bond_cc*cos(0.5*ANGLE_6RING) + tmp1;
 747     tmp1 *= 2;
 748     a     = b = -bond_ch / tmp1;
 749     /* HE1 and HE2: */
 750     add_vsite3_param(&plist[F_VSITE3],
 751                      ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
 752     add_vsite3_param(&plist[F_VSITE3],
 753                      ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
 754     /* CD1, CD2 and CZ: */
 755     a = b = tmp2 / tmp1;
 756     add_vsite3_param(&plist[F_VSITE3],
 757                      ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
 758     add_vsite3_param(&plist[F_VSITE3],
 759                      ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
 760     if (bDoZ)
 761     {
 762         add_vsite3_param(&plist[F_VSITE3],
 763                          ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
 764     }
 765     /* HD1, HD2 and HZ: */
 766     a = b = ( bond_ch + tmp2 ) / tmp1;
 767     add_vsite3_param(&plist[F_VSITE3],
 768                      ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
 769     add_vsite3_param(&plist[F_VSITE3],
 770                      ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
 771     if (bDoZ)
 772     {
 773         add_vsite3_param(&plist[F_VSITE3],
 774                          ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
 775     }
 776
 777     return nvsite;
 778 }
 779
 780 static int gen_vsites_phe(t_atoms *at, int *vsite_type[], t_params plist[],
 781                           int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
 782 {
 783     real bond_cc, bond_ch;
 784     real xcom, mtot;
 785     int  i;
 786     /* these MUST correspond to the atnms array in do_vsite_aromatics! */
 787     enum {
 788         atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
 789         atCZ, atHZ, atNR
 790     };
 791     real x[atNR], y[atNR];
 792     /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
 793      * (angle is always 120 degrees).
 794      */
 795     bond_cc = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "CE1");
 796     bond_ch = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "HD1");
 797
 798     x[atCG]  = -bond_cc+bond_cc*cos(ANGLE_6RING);
 799     y[atCG]  = 0;
 800     x[atCD1] = -bond_cc;
 801     y[atCD1] = bond_cc*sin(0.5*ANGLE_6RING);
 802     x[atHD1] = x[atCD1]+bond_ch*cos(ANGLE_6RING);
 803     y[atHD1] = y[atCD1]+bond_ch*sin(ANGLE_6RING);
 804     x[atCE1] = 0;
 805     y[atCE1] = y[atCD1];
 806     x[atHE1] = x[atCE1]-bond_ch*cos(ANGLE_6RING);
 807     y[atHE1] = y[atCE1]+bond_ch*sin(ANGLE_6RING);
 808     x[atCD2] = x[atCD1];
 809     y[atCD2] = -y[atCD1];
 810     x[atHD2] = x[atHD1];
 811     y[atHD2] = -y[atHD1];
 812     x[atCE2] = x[atCE1];
 813     y[atCE2] = -y[atCE1];
 814     x[atHE2] = x[atHE1];
 815     y[atHE2] = -y[atHE1];
 816     x[atCZ]  = bond_cc*cos(0.5*ANGLE_6RING);
 817     y[atCZ]  = 0;
 818     x[atHZ]  = x[atCZ]+bond_ch;
 819     y[atHZ]  = 0;
 820
 821     xcom = mtot = 0;
 822     for (i = 0; i < atNR; i++)
 823     {
 824         xcom += x[i]*at->atom[ats[i]].m;
 825         mtot += at->atom[ats[i]].m;
 826     }
 827     xcom /= mtot;
 828
 829     return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, TRUE);
 830 }
 831
 832 static void calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj,
 833                               real xk, real yk, real *a, real *b)
 834 {
 835     /* determine parameters by solving the equation system, since we know the
 836      * virtual site coordinates here.
 837      */
 838     real dx_ij, dx_ik, dy_ij, dy_ik;
 839     real b_ij, b_ik;
 840
 841     dx_ij = xj-xi;
 842     dy_ij = yj-yi;
 843     dx_ik = xk-xi;
 844     dy_ik = yk-yi;
 845     b_ij  = sqrt(dx_ij*dx_ij+dy_ij*dy_ij);
 846     b_ik  = sqrt(dx_ik*dx_ik+dy_ik*dy_ik);
 847
 848     *a = ( (xd-xi)*dy_ik - dx_ik*(yd-yi) ) / (dx_ij*dy_ik - dx_ik*dy_ij);
 849     *b = ( yd - yi - (*a)*dy_ij ) / dy_ik;
 850 }
 851
 852
 853 static int gen_vsites_trp(gpp_atomtype_t atype, rvec *newx[],
 854                           t_atom *newatom[], char ***newatomname[],
 855                           int *o2n[], int *newvsite_type[], int *newcgnr[],
 856                           t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
 857                           t_atoms *at, int *vsite_type[], t_params plist[],
 858                           int nrfound, int *ats, int add_shift,
 859                           t_vsitetop *vsitetop, int nvsitetop)
 860 {
 861 #define NMASS 2
 862     /* these MUST correspond to the atnms array in do_vsite_aromatics! */
 863     enum {
 864         atCB,  atCG,  atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
 865         atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR
 866     };
 867     /* weights for determining the COM's of both rings (M1 and M2): */
 868     real mw[NMASS][atNR] = {
 869         {   0,     1,     1,     1,   0.5,     1,     1,   0.5,     0,     0,
 870             0,     0,     0,     0,     0,     0 },
 871         {   0,     0,     0,     0,   0.5,     0,     0,   0.5,     1,     1,
 872             1,     1,     1,     1,     1,     1 }
 873     };
 874
 875     real xi[atNR], yi[atNR];
 876     real xcom[NMASS], ycom[NMASS], I, alpha;
 877     real lineA, lineB, dist;
 878     real b_CD2_CE2, b_NE1_CE2, b_CG_CD2, b_CH2_HH2, b_CE2_CZ2;
 879     real b_NE1_HE1, b_CD2_CE3, b_CE3_CZ3, b_CB_CG;
 880     real b_CZ2_CH2, b_CZ2_HZ2, b_CD1_HD1, b_CE3_HE3;
 881     real b_CG_CD1, b_CZ3_HZ3;
 882     real a_NE1_CE2_CD2, a_CE2_CD2_CG, a_CB_CG_CD2, a_CE2_CD2_CE3;
 883     real a_CB_CG_CD1, a_CD2_CG_CD1, a_CE2_CZ2_HZ2, a_CZ2_CH2_HH2;
 884     real a_CD2_CE2_CZ2, a_CD2_CE3_CZ3, a_CE3_CZ3_HZ3, a_CG_CD1_HD1;
 885     real a_CE2_CZ2_CH2, a_HE1_NE1_CE2, a_CD2_CE3_HE3;
 886     real xM[NMASS];
 887     int  atM[NMASS], tpM, i, i0, j, nvsite;
 888     real mwtot, mtot, mM[NMASS], dCBM1, dCBM2, dM1M2;
 889     real a, b, c[MAXFORCEPARAM];
 890     rvec r_ij, r_ik, t1, t2;
 891     char name[10];
 892
 893     if (atNR != nrfound)
 894     {
 895         gmx_incons("atom types in gen_vsites_trp");
 896     }
 897     /* Get geometry from database */
 898     b_CD2_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE2");
 899     b_NE1_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "CE2");
 900     b_CG_CD1  = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD1");
 901     b_CG_CD2  = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD2");
 902     b_CB_CG   = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CB", "CG");
 903     b_CE2_CZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE2", "CZ2");
 904     b_CD2_CE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE3");
 905     b_CE3_CZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "CZ3");
 906     b_CZ2_CH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "CH2");
 907
 908     b_CD1_HD1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD1", "HD1");
 909     b_CZ2_HZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "HZ2");
 910     b_NE1_HE1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "HE1");
 911     b_CH2_HH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CH2", "HH2");
 912     b_CE3_HE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "HE3");
 913     b_CZ3_HZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ3", "HZ3");
 914
 915     a_NE1_CE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "NE1", "CE2", "CD2");
 916     a_CE2_CD2_CG  = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CG");
 917     a_CB_CG_CD2   = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD2");
 918     a_CD2_CG_CD1  = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CG", "CD1");
 919     a_CB_CG_CD1   = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD1");
 920
 921     a_CE2_CD2_CE3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CE3");
 922     a_CD2_CE2_CZ2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE2", "CZ2");
 923     a_CD2_CE3_CZ3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "CZ3");
 924     a_CE3_CZ3_HZ3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE3", "CZ3", "HZ3");
 925     a_CZ2_CH2_HH2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CZ2", "CH2", "HH2");
 926     a_CE2_CZ2_HZ2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "HZ2");
 927     a_CE2_CZ2_CH2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "CH2");
 928     a_CG_CD1_HD1  = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CG", "CD1", "HD1");
 929     a_HE1_NE1_CE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "HE1", "NE1", "CE2");
 930     a_CD2_CE3_HE3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "HE3");
 931
 932     /* Calculate local coordinates.
 933      * y-axis (x=0) is the bond CD2-CE2.
 934      * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
 935      * intersects the middle of the bond.
 936      */
 937     xi[atCD2] = 0;
 938     yi[atCD2] = -0.5*b_CD2_CE2;
 939
 940     xi[atCE2] = 0;
 941     yi[atCE2] = 0.5*b_CD2_CE2;
 942
 943     xi[atNE1] = -b_NE1_CE2*sin(a_NE1_CE2_CD2);
 944     yi[atNE1] = yi[atCE2]-b_NE1_CE2*cos(a_NE1_CE2_CD2);
 945
 946     xi[atCG] = -b_CG_CD2*sin(a_CE2_CD2_CG);
 947     yi[atCG] = yi[atCD2]+b_CG_CD2*cos(a_CE2_CD2_CG);
 948
 949     alpha    = a_CE2_CD2_CG + M_PI - a_CB_CG_CD2;
 950     xi[atCB] = xi[atCG]-b_CB_CG*sin(alpha);
 951     yi[atCB] = yi[atCG]+b_CB_CG*cos(alpha);
 952
 953     alpha     = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI;
 954     xi[atCD1] = xi[atCG]-b_CG_CD1*sin(alpha);
 955     yi[atCD1] = yi[atCG]+b_CG_CD1*cos(alpha);
 956
 957     xi[atCE3] = b_CD2_CE3*sin(a_CE2_CD2_CE3);
 958     yi[atCE3] = yi[atCD2]+b_CD2_CE3*cos(a_CE2_CD2_CE3);
 959
 960     xi[atCZ2] = b_CE2_CZ2*sin(a_CD2_CE2_CZ2);
 961     yi[atCZ2] = yi[atCE2]-b_CE2_CZ2*cos(a_CD2_CE2_CZ2);
 962
 963     alpha     = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI;
 964     xi[atCZ3] = xi[atCE3]+b_CE3_CZ3*sin(alpha);
 965     yi[atCZ3] = yi[atCE3]+b_CE3_CZ3*cos(alpha);
 966
 967     alpha     = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI;
 968     xi[atCH2] = xi[atCZ2]+b_CZ2_CH2*sin(alpha);
 969     yi[atCH2] = yi[atCZ2]-b_CZ2_CH2*cos(alpha);
 970
 971     /* hydrogens */
 972     alpha     = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
 973     xi[atHD1] = xi[atCD1]-b_CD1_HD1*sin(alpha);
 974     yi[atHD1] = yi[atCD1]+b_CD1_HD1*cos(alpha);
 975
 976     alpha     = a_NE1_CE2_CD2 + M_PI - a_HE1_NE1_CE2;
 977     xi[atHE1] = xi[atNE1]-b_NE1_HE1*sin(alpha);
 978     yi[atHE1] = yi[atNE1]-b_NE1_HE1*cos(alpha);
 979
 980     alpha     = a_CE2_CD2_CE3 + M_PI - a_CD2_CE3_HE3;
 981     xi[atHE3] = xi[atCE3]+b_CE3_HE3*sin(alpha);
 982     yi[atHE3] = yi[atCE3]+b_CE3_HE3*cos(alpha);
 983
 984     alpha     = a_CD2_CE2_CZ2 + M_PI - a_CE2_CZ2_HZ2;
 985     xi[atHZ2] = xi[atCZ2]+b_CZ2_HZ2*sin(alpha);
 986     yi[atHZ2] = yi[atCZ2]-b_CZ2_HZ2*cos(alpha);
 987
 988     alpha     = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
 989     xi[atHZ3] = xi[atCZ3]+b_CZ3_HZ3*sin(alpha);
 990     yi[atHZ3] = yi[atCZ3]+b_CZ3_HZ3*cos(alpha);
 991
 992     alpha     = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
 993     xi[atHH2] = xi[atCH2]+b_CH2_HH2*sin(alpha);
 994     yi[atHH2] = yi[atCH2]-b_CH2_HH2*cos(alpha);
 995
 996     /* Determine coeff. for the line CB-CG */
 997     lineA = (yi[atCB]-yi[atCG])/(xi[atCB]-xi[atCG]);
 998     lineB = yi[atCG]-lineA*xi[atCG];
 999
1000     /* Calculate masses for each ring and put it on the dummy masses */
1001     for (j = 0; j < NMASS; j++)
1002     {
1003         mM[j] = xcom[j] = ycom[j] = 0;
1004     }
1005     for (i = 0; i < atNR; i++)
1006     {
1007         if (i != atCB)
1008         {
1009             for (j = 0; j < NMASS; j++)
1010             {
1011                 mM[j]   += mw[j][i] * at->atom[ats[i]].m;
1012                 xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
1013                 ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
1014             }
1015         }
1016     }
1017     for (j = 0; j < NMASS; j++)
1018     {
1019         xcom[j] /= mM[j];
1020         ycom[j] /= mM[j];
1021     }
1022
1023     /* get dummy mass type */
1024     tpM = vsite_nm2type("MW", atype);
1025     /* make space for 2 masses: shift all atoms starting with CB */
1026     i0 = ats[atCB];
1027     for (j = 0; j < NMASS; j++)
1028     {
1029         atM[j] = i0+*nadd+j;
1030     }
1031     if (debug)
1032     {
1033         fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0]+1);
1034     }
1035     *nadd += NMASS;
1036     for (j = i0; j < at->nr; j++)
1037     {
1038         (*o2n)[j] = j+*nadd;
1039     }
1040     srenew(*newx, at->nr+*nadd);
1041     srenew(*newatom, at->nr+*nadd);
1042     srenew(*newatomname, at->nr+*nadd);
1043     srenew(*newvsite_type, at->nr+*nadd);
1044     srenew(*newcgnr, at->nr+*nadd);
1045     for (j = 0; j < NMASS; j++)
1046     {
1047         (*newatomname)[at->nr+*nadd-1-j] = NULL;
1048     }
1049
1050     /* Dummy masses will be placed at the center-of-mass in each ring. */
1051
1052     /* calc initial position for dummy masses in real (non-local) coordinates.
1053      * Cheat by using the routine to calculate virtual site parameters. It is
1054      * much easier when we have the coordinates expressed in terms of
1055      * CB, CG, CD2.
1056      */
1057     rvec_sub(x[ats[atCB]], x[ats[atCG]], r_ij);
1058     rvec_sub(x[ats[atCD2]], x[ats[atCG]], r_ik);
1059     calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1060                       xi[atCD2], yi[atCD2], &a, &b);
1061     svmul(a, r_ij, t1);
1062     svmul(b, r_ik, t2);
1063     rvec_add(t1, t2, t1);
1064     rvec_add(t1, x[ats[atCG]], (*newx)[atM[0]]);
1065
1066     calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1067                       xi[atCD2], yi[atCD2], &a, &b);
1068     svmul(a, r_ij, t1);
1069     svmul(b, r_ik, t2);
1070     rvec_add(t1, t2, t1);
1071     rvec_add(t1, x[ats[atCG]], (*newx)[atM[1]]);
1072
1073     /* set parameters for the masses */
1074     for (j = 0; j < NMASS; j++)
1075     {
1076         sprintf(name, "MW%d", j+1);
1077         (*newatomname)  [atM[j]]        = put_symtab(symtab, name);
1078         (*newatom)      [atM[j]].m      = (*newatom)[atM[j]].mB    = mM[j];
1079         (*newatom)      [atM[j]].q      = (*newatom)[atM[j]].qB    = 0.0;
1080         (*newatom)      [atM[j]].type   = (*newatom)[atM[j]].typeB = tpM;
1081         (*newatom)      [atM[j]].ptype  = eptAtom;
1082         (*newatom)      [atM[j]].resind = at->atom[i0].resind;
1083         (*newvsite_type)[atM[j]]        = NOTSET;
1084         (*newcgnr)      [atM[j]]        = (*cgnr)[i0];
1085     }
1086     /* renumber cgnr: */
1087     for (i = i0; i < at->nr; i++)
1088     {
1089         (*cgnr)[i]++;
1090     }
1091
1092     /* constraints between CB, M1 and M2 */
1093     /* 'add_shift' says which atoms won't be renumbered afterwards */
1094     dCBM1 = sqrt( sqr(xcom[0]-xi[atCB]) + sqr(ycom[0]-yi[atCB]) );
1095     dM1M2 = sqrt( sqr(xcom[0]-xcom[1]) + sqr(ycom[0]-ycom[1]) );
1096     dCBM2 = sqrt( sqr(xcom[1]-xi[atCB]) + sqr(ycom[1]-yi[atCB]) );
1097     my_add_param(&(plist[F_CONSTRNC]), ats[atCB],       add_shift+atM[0], dCBM1);
1098     my_add_param(&(plist[F_CONSTRNC]), ats[atCB],       add_shift+atM[1], dCBM2);
1099     my_add_param(&(plist[F_CONSTRNC]), add_shift+atM[0], add_shift+atM[1], dM1M2);
1100
1101     /* rest will be vsite3 */
1102     nvsite = 0;
1103     for (i = 0; i < atNR; i++)
1104     {
1105         if (i != atCB)
1106         {
1107             at->atom[ats[i]].m    = at->atom[ats[i]].mB = 0;
1108             (*vsite_type)[ats[i]] = F_VSITE3;
1109             nvsite++;
1110         }
1111     }
1112
1113     /* now define all vsites from M1, M2, CB, ie:
1114        r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
1115     for (i = 0; i < atNR; i++)
1116     {
1117         if ( (*vsite_type)[ats[i]] == F_VSITE3)
1118         {
1119             calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB], &a, &b);
1120             add_vsite3_param(&plist[F_VSITE3],
1121                              ats[i], add_shift+atM[0], add_shift+atM[1], ats[atCB], a, b);
1122         }
1123     }
1124     return nvsite;
1125 #undef NMASS
1126 }
1127
1128
1129 static int gen_vsites_tyr(gpp_atomtype_t atype, rvec *newx[],
1130                           t_atom *newatom[], char ***newatomname[],
1131                           int *o2n[], int *newvsite_type[], int *newcgnr[],
1132                           t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
1133                           t_atoms *at, int *vsite_type[], t_params plist[],
1134                           int nrfound, int *ats, int add_shift,
1135                           t_vsitetop *vsitetop, int nvsitetop)
1136 {
1137     int  nvsite, i, i0, j, atM, tpM;
1138     real dCGCE, dCEOH, dCGM, tmp1, a, b;
1139     real bond_cc, bond_ch, bond_co, bond_oh, angle_coh;
1140     real xcom, mtot;
1141     real vmass, vdist, mM;
1142     rvec r1;
1143     char name[10];
1144
1145     /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1146     enum {
1147         atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
1148         atCZ, atOH, atHH, atNR
1149     };
1150     real xi[atNR], yi[atNR];
1151     /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
1152        rest gets virtualized.
1153        Now we have two linked triangles with one improper keeping them flat */
1154     if (atNR != nrfound)
1155     {
1156         gmx_incons("Number of atom types in gen_vsites_tyr");
1157     }
1158
1159     /* Aromatic rings have 6-fold symmetry, so we only need one bond length
1160      * for the ring part (angle is always 120 degrees).
1161      */
1162     bond_cc   = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "CE1");
1163     bond_ch   = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "HD1");
1164     bond_co   = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CZ", "OH");
1165     bond_oh   = get_ddb_bond(vsitetop, nvsitetop, "TYR", "OH", "HH");
1166     angle_coh = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TYR", "CZ", "OH", "HH");
1167
1168     xi[atCG]  = -bond_cc+bond_cc*cos(ANGLE_6RING);
1169     yi[atCG]  = 0;
1170     xi[atCD1] = -bond_cc;
1171     yi[atCD1] = bond_cc*sin(0.5*ANGLE_6RING);
1172     xi[atHD1] = xi[atCD1]+bond_ch*cos(ANGLE_6RING);
1173     yi[atHD1] = yi[atCD1]+bond_ch*sin(ANGLE_6RING);
1174     xi[atCE1] = 0;
1175     yi[atCE1] = yi[atCD1];
1176     xi[atHE1] = xi[atCE1]-bond_ch*cos(ANGLE_6RING);
1177     yi[atHE1] = yi[atCE1]+bond_ch*sin(ANGLE_6RING);
1178     xi[atCD2] = xi[atCD1];
1179     yi[atCD2] = -yi[atCD1];
1180     xi[atHD2] = xi[atHD1];
1181     yi[atHD2] = -yi[atHD1];
1182     xi[atCE2] = xi[atCE1];
1183     yi[atCE2] = -yi[atCE1];
1184     xi[atHE2] = xi[atHE1];
1185     yi[atHE2] = -yi[atHE1];
1186     xi[atCZ]  = bond_cc*cos(0.5*ANGLE_6RING);
1187     yi[atCZ]  = 0;
1188     xi[atOH]  = xi[atCZ]+bond_co;
1189     yi[atOH]  = 0;
1190
1191     xcom = mtot = 0;
1192     for (i = 0; i < atOH; i++)
1193     {
1194         xcom += xi[i]*at->atom[ats[i]].m;
1195         mtot += at->atom[ats[i]].m;
1196     }
1197     xcom /= mtot;
1198
1199     /* first do 6 ring as default,
1200        except CZ (we'll do that different) and HZ (we don't have that): */
1201     nvsite = gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, FALSE);
1202
1203     /* then construct CZ from the 2nd triangle */
1204     /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1205     a = b = 0.5 * bond_co / ( bond_co - bond_cc*cos(ANGLE_6RING) );
1206     add_vsite3_param(&plist[F_VSITE3],
1207                      ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
1208     at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
1209
1210     /* constraints between CE1, CE2 and OH */
1211     dCGCE = sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
1212     dCEOH = sqrt( cosrule(bond_cc, bond_co, ANGLE_6RING) );
1213     my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atOH], dCEOH);
1214     my_add_param(&(plist[F_CONSTRNC]), ats[atCE2], ats[atOH], dCEOH);
1215
1216     /* We also want to constrain the angle C-O-H, but since CZ is constructed
1217      * we need to introduce a constraint to CG.
1218      * CG is much further away, so that will lead to instabilities in LINCS
1219      * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1220      * the use of lincs_order=8 we introduce a dummy mass three times further
1221      * away from OH than HH. The mass is accordingly a third, with the remaining
1222      * 2/3 moved to OH. This shouldnt cause any problems since the forces will
1223      * apply to the HH constructed atom and not directly on the virtual mass.
1224      */
1225
1226     vdist                   = 2.0*bond_oh;
1227     mM                      = at->atom[ats[atHH]].m/2.0;
1228     at->atom[ats[atOH]].m  += mM; /* add 1/2 of original H mass */
1229     at->atom[ats[atOH]].mB += mM; /* add 1/2 of original H mass */
1230     at->atom[ats[atHH]].m   = at->atom[ats[atHH]].mB = 0;
1231
1232     /* get dummy mass type */
1233     tpM = vsite_nm2type("MW", atype);
1234     /* make space for 1 mass: shift HH only */
1235     i0  = ats[atHH];
1236     atM = i0+*nadd;
1237     if (debug)
1238     {
1239         fprintf(stderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0]+1);
1240     }
1241     (*nadd)++;
1242     for (j = i0; j < at->nr; j++)
1243     {
1244         (*o2n)[j] = j+*nadd;
1245     }
1246     srenew(*newx, at->nr+*nadd);
1247     srenew(*newatom, at->nr+*nadd);
1248     srenew(*newatomname, at->nr+*nadd);
1249     srenew(*newvsite_type, at->nr+*nadd);
1250     srenew(*newcgnr, at->nr+*nadd);
1251     (*newatomname)[at->nr+*nadd-1] = NULL;
1252
1253     /* Calc the dummy mass initial position */
1254     rvec_sub(x[ats[atHH]], x[ats[atOH]], r1);
1255     svmul(2.0, r1, r1);
1256     rvec_add(r1, x[ats[atHH]], (*newx)[atM]);
1257
1258     strcpy(name, "MW1");
1259     (*newatomname)  [atM]        = put_symtab(symtab, name);
1260     (*newatom)      [atM].m      = (*newatom)[atM].mB    = mM;
1261     (*newatom)      [atM].q      = (*newatom)[atM].qB    = 0.0;
1262     (*newatom)      [atM].type   = (*newatom)[atM].typeB = tpM;
1263     (*newatom)      [atM].ptype  = eptAtom;
1264     (*newatom)      [atM].resind = at->atom[i0].resind;
1265     (*newvsite_type)[atM]        = NOTSET;
1266     (*newcgnr)      [atM]        = (*cgnr)[i0];
1267     /* renumber cgnr: */
1268     for (i = i0; i < at->nr; i++)
1269     {
1270         (*cgnr)[i]++;
1271     }
1272
1273     (*vsite_type)[ats[atHH]] = F_VSITE2;
1274     nvsite++;
1275     /* assume we also want the COH angle constrained: */
1276     tmp1 = bond_cc*cos(0.5*ANGLE_6RING) + dCGCE*sin(ANGLE_6RING*0.5) + bond_co;
1277     dCGM = sqrt( cosrule(tmp1, vdist, angle_coh) );
1278     my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift+atM, dCGM);
1279     my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift+atM, vdist);
1280
1281     add_vsite2_param(&plist[F_VSITE2],
1282                      ats[atHH], ats[atOH], add_shift+atM, 1.0/2.0);
1283     return nvsite;
1284 }
1285
1286 static int gen_vsites_his(t_atoms *at, int *vsite_type[], t_params plist[],
1287                           int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
1288 {
1289     int  nvsite, i;
1290     real a, b, alpha, dCGCE1, dCGNE2;
1291     real sinalpha, cosalpha;
1292     real xcom, ycom, mtot;
1293     real mG, mrest, mCE1, mNE2;
1294     real b_CG_ND1, b_ND1_CE1, b_CE1_NE2, b_CG_CD2, b_CD2_NE2;
1295     real b_ND1_HD1, b_NE2_HE2, b_CE1_HE1, b_CD2_HD2;
1296     real a_CG_ND1_CE1, a_CG_CD2_NE2, a_ND1_CE1_NE2, a_CE1_NE2_CD2;
1297     real a_NE2_CE1_HE1, a_NE2_CD2_HD2, a_CE1_ND1_HD1, a_CE1_NE2_HE2;
1298     char resname[10];
1299
1300     /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1301     enum {
1302         atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR
1303     };
1304     real x[atNR], y[atNR];
1305
1306     /* CG, CE1 and NE2 stay, each gets part of the total mass,
1307        rest gets virtualized */
1308     /* check number of atoms, 3 hydrogens may be missing: */
1309     /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1310      * Don't understand the above logic. Shouldn't it be && rather than || ???
1311      */
1312     if ((nrfound < atNR-3) || (nrfound > atNR))
1313     {
1314         gmx_incons("Generating vsites for HIS");
1315     }
1316
1317     /* avoid warnings about uninitialized variables */
1318     b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 =
1319                         a_NE2_CD2_HD2 = a_CE1_ND1_HD1 = a_CE1_NE2_HE2 = 0;
1320
1321     if (ats[atHD1] != NOTSET)
1322     {
1323         if (ats[atHE2] != NOTSET)
1324         {
1325             sprintf(resname, "HISH");
1326         }
1327         else
1328         {
1329             sprintf(resname, "HISA");
1330         }
1331     }
1332     else
1333     {
1334         sprintf(resname, "HISB");
1335     }
1336
1337     /* Get geometry from database */
1338     b_CG_ND1      = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "ND1");
1339     b_ND1_CE1     = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "CE1");
1340     b_CE1_NE2     = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "NE2");
1341     b_CG_CD2      = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "CD2");
1342     b_CD2_NE2     = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "NE2");
1343     a_CG_ND1_CE1  = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "ND1", "CE1");
1344     a_CG_CD2_NE2  = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "CD2", "NE2");
1345     a_ND1_CE1_NE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "ND1", "CE1", "NE2");
1346     a_CE1_NE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "CD2");
1347
1348     if (ats[atHD1] != NOTSET)
1349     {
1350         b_ND1_HD1     = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "HD1");
1351         a_CE1_ND1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "ND1", "HD1");
1352     }
1353     if (ats[atHE2] != NOTSET)
1354     {
1355         b_NE2_HE2     = get_ddb_bond(vsitetop, nvsitetop, resname, "NE2", "HE2");
1356         a_CE1_NE2_HE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "HE2");
1357     }
1358     if (ats[atHD2] != NOTSET)
1359     {
1360         b_CD2_HD2     = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "HD2");
1361         a_NE2_CD2_HD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CD2", "HD2");
1362     }
1363     if (ats[atHE1] != NOTSET)
1364     {
1365         b_CE1_HE1     = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "HE1");
1366         a_NE2_CE1_HE1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CE1", "HE1");
1367     }
1368
1369     /* constraints between CG, CE1 and NE1 */
1370     dCGCE1   = sqrt( cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1) );
1371     dCGNE2   = sqrt( cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2) );
1372
1373     my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE1);
1374     my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atNE2], dCGNE2);
1375     /* we already have a constraint CE1-NE2, so we don't add it again */
1376
1377     /* calculate the positions in a local frame of reference.
1378      * The x-axis is the line from CG that makes a right angle
1379      * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1380      */
1381     /* First calculate the x-axis intersection with y-axis (=yCE1).
1382      * Get cos(angle CG-CE1-NE2) :
1383      */
1384     cosalpha = acosrule(dCGNE2, dCGCE1, b_CE1_NE2);
1385     x[atCE1] = 0;
1386     y[atCE1] = cosalpha*dCGCE1;
1387     x[atNE2] = 0;
1388     y[atNE2] = y[atCE1]-b_CE1_NE2;
1389     sinalpha = sqrt(1-cosalpha*cosalpha);
1390     x[atCG]  = -sinalpha*dCGCE1;
1391     y[atCG]  = 0;
1392     x[atHE1] = x[atHE2] = x[atHD1] = x[atHD2] = 0;
1393     y[atHE1] = y[atHE2] = y[atHD1] = y[atHD2] = 0;
1394
1395     /* calculate ND1 and CD2 positions from CE1 and NE2 */
1396
1397     x[atND1] = -b_ND1_CE1*sin(a_ND1_CE1_NE2);
1398     y[atND1] = y[atCE1]-b_ND1_CE1*cos(a_ND1_CE1_NE2);
1399
1400     x[atCD2] = -b_CD2_NE2*sin(a_CE1_NE2_CD2);
1401     y[atCD2] = y[atNE2]+b_CD2_NE2*cos(a_CE1_NE2_CD2);
1402
1403     /* And finally the hydrogen positions */
1404     if (ats[atHE1] != NOTSET)
1405     {
1406         x[atHE1] = x[atCE1] + b_CE1_HE1*sin(a_NE2_CE1_HE1);
1407         y[atHE1] = y[atCE1] - b_CE1_HE1*cos(a_NE2_CE1_HE1);
1408     }
1409     /* HD2 - first get (ccw) angle from (positive) y-axis */
1410     if (ats[atHD2] != NOTSET)
1411     {
1412         alpha    = a_CE1_NE2_CD2 + M_PI - a_NE2_CD2_HD2;
1413         x[atHD2] = x[atCD2] - b_CD2_HD2*sin(alpha);
1414         y[atHD2] = y[atCD2] + b_CD2_HD2*cos(alpha);
1415     }
1416     if (ats[atHD1] != NOTSET)
1417     {
1418         /* HD1 - first get (cw) angle from (positive) y-axis */
1419         alpha    = a_ND1_CE1_NE2 + M_PI - a_CE1_ND1_HD1;
1420         x[atHD1] = x[atND1] - b_ND1_HD1*sin(alpha);
1421         y[atHD1] = y[atND1] - b_ND1_HD1*cos(alpha);
1422     }
1423     if (ats[atHE2] != NOTSET)
1424     {
1425         x[atHE2] = x[atNE2] + b_NE2_HE2*sin(a_CE1_NE2_HE2);
1426         y[atHE2] = y[atNE2] + b_NE2_HE2*cos(a_CE1_NE2_HE2);
1427     }
1428     /* Have all coordinates now */
1429
1430     /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1431      * set the rest to vsite3
1432      */
1433     mtot   = xcom = ycom = 0;
1434     nvsite = 0;
1435     for (i = 0; i < atNR; i++)
1436     {
1437         if (ats[i] != NOTSET)
1438         {
1439             mtot += at->atom[ats[i]].m;
1440             xcom += x[i]*at->atom[ats[i]].m;
1441             ycom += y[i]*at->atom[ats[i]].m;
1442             if (i != atCG && i != atCE1 && i != atNE2)
1443             {
1444                 at->atom[ats[i]].m    = at->atom[ats[i]].mB = 0;
1445                 (*vsite_type)[ats[i]] = F_VSITE3;
1446                 nvsite++;
1447             }
1448         }
1449     }
1450     if (nvsite+3 != nrfound)
1451     {
1452         gmx_incons("Generating vsites for HIS");
1453     }
1454
1455     xcom /= mtot;
1456     ycom /= mtot;
1457
1458     /* distribute mass so that com stays the same */
1459     mG    = xcom*mtot/x[atCG];
1460     mrest = mtot-mG;
1461     mCE1  = (ycom-y[atNE2])*mrest/(y[atCE1]-y[atNE2]);
1462     mNE2  = mrest-mCE1;
1463
1464     at->atom[ats[atCG]].m  = at->atom[ats[atCG]].mB = mG;
1465     at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
1466     at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
1467
1468     /* HE1 */
1469     if (ats[atHE1] != NOTSET)
1470     {
1471         calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1472                           x[atCG], y[atCG], &a, &b);
1473         add_vsite3_param(&plist[F_VSITE3],
1474                          ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1475     }
1476     /* HE2 */
1477     if (ats[atHE2] != NOTSET)
1478     {
1479         calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1480                           x[atCG], y[atCG], &a, &b);
1481         add_vsite3_param(&plist[F_VSITE3],
1482                          ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1483     }
1484
1485     /* ND1 */
1486     calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1487                       x[atCG], y[atCG], &a, &b);
1488     add_vsite3_param(&plist[F_VSITE3],
1489                      ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1490
1491     /* CD2 */
1492     calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1493                       x[atCG], y[atCG], &a, &b);
1494     add_vsite3_param(&plist[F_VSITE3],
1495                      ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1496
1497     /* HD1 */
1498     if (ats[atHD1] != NOTSET)
1499     {
1500         calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1501                           x[atCG], y[atCG], &a, &b);
1502         add_vsite3_param(&plist[F_VSITE3],
1503                          ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1504     }
1505     /* HD2 */
1506     if (ats[atHD2] != NOTSET)
1507     {
1508         calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1509                           x[atCG], y[atCG], &a, &b);
1510         add_vsite3_param(&plist[F_VSITE3],
1511                          ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1512     }
1513     return nvsite;
1514 }
1515
1516 static gmx_bool is_vsite(int vsite_type)
1517 {
1518     if (vsite_type == NOTSET)
1519     {
1520         return FALSE;
1521     }
1522     switch (abs(vsite_type) )
1523     {
1524         case F_VSITE3:
1525         case F_VSITE3FD:
1526         case F_VSITE3OUT:
1527         case F_VSITE3FAD:
1528         case F_VSITE4FD:
1529         case F_VSITE4FDN:
1530             return TRUE;
1531         default:
1532             return FALSE;
1533     }
1534 }
1535
1536 static char atomnamesuffix[] = "1234";
1537
1538 void do_vsites(int nrtp, t_restp rtp[], gpp_atomtype_t atype,
1539                t_atoms *at, t_symtab *symtab, rvec *x[],
1540                t_params plist[], int *vsite_type[], int *cgnr[],
1541                real mHmult, gmx_bool bVsiteAromatics,
1542                const char *ffdir)
1543 {
1544 #define MAXATOMSPERRESIDUE 16
1545     int               i, j, k, m, i0, ni0, whatres, resind, add_shift, ftype, nvsite, nadd;
1546     int               ai, aj, ak, al;
1547     int               nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
1548     int               Hatoms[4], heavies[4], bb;
1549     gmx_bool          bWARNING, bAddVsiteParam, bFirstWater;
1550     matrix            tmpmat;
1551     gmx_bool         *bResProcessed;
1552     real              mHtot, mtot, fact, fact2;
1553     rvec              rpar, rperp, temp;
1554     char              name[10], tpname[32], nexttpname[32], *ch;
1555     rvec             *newx;
1556     int              *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
1557     t_atom           *newatom;
1558     t_params         *params;
1559     char           ***newatomname;
1560     char             *resnm = NULL;
1561     int               ndb, f;
1562     char            **db;
1563     int               nvsiteconf, nvsitetop, cmplength;
1564     gmx_bool          isN, planarN, bFound;
1565     gmx_residuetype_t rt;
1566
1567     t_vsiteconf      *vsiteconflist;
1568     /* pointer to a list of CH3/NH3/NH2 configuration entries.
1569      * See comments in read_vsite_database. It isnt beautiful,
1570      * but it had to be fixed, and I dont even want to try to
1571      * maintain this part of the code...
1572      */
1573     t_vsitetop *vsitetop;
1574     /* Pointer to a list of geometry (bond/angle) entries for
1575      * residues like PHE, TRP, TYR, HIS, etc., where we need
1576      * to know the geometry to construct vsite aromatics.
1577      * Note that equilibrium geometry isnt necessarily the same
1578      * as the individual bond and angle values given in the
1579      * force field (rings can be strained).
1580      */
1581
1582     /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1583        PHE, TRP, TYR and HIS to a construction of virtual sites */
1584     enum                    {
1585         resPHE, resTRP, resTYR, resHIS, resNR
1586     };
1587     const char *resnms[resNR]   = {   "PHE",  "TRP",  "TYR",  "HIS" };
1588     /* Amber03 alternative names for termini */
1589     const char *resnmsN[resNR]  = {  "NPHE", "NTRP", "NTYR", "NHIS" };
1590     const char *resnmsC[resNR]  = {  "CPHE", "CTRP", "CTYR", "CHIS" };
1591     /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1592     gmx_bool    bPartial[resNR]  = {  FALSE,  FALSE,  FALSE,   TRUE  };
1593     /* the atnms for every residue MUST correspond to the enums in the
1594        gen_vsites_* (one for each residue) routines! */
1595     /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1596     const char *atnms[resNR][MAXATOMSPERRESIDUE+1] = {
1597         { "CG", /* PHE */
1598           "CD1", "HD1", "CD2", "HD2",
1599           "CE1", "HE1", "CE2", "HE2",
1600           "CZ", "HZ", NULL },
1601         { "CB", /* TRP */
1602           "CG",
1603           "CD1", "HD1", "CD2",
1604           "NE1", "HE1", "CE2", "CE3", "HE3",
1605           "CZ2", "HZ2", "CZ3", "HZ3",
1606           "CH2", "HH2", NULL },
1607         { "CG", /* TYR */
1608           "CD1", "HD1", "CD2", "HD2",
1609           "CE1", "HE1", "CE2", "HE2",
1610           "CZ", "OH", "HH", NULL },
1611         { "CG", /* HIS */
1612           "ND1", "HD1", "CD2", "HD2",
1613           "CE1", "HE1", "NE2", "HE2", NULL }
1614     };
1615
1616     if (debug)
1617     {
1618         printf("Searching for atoms to make virtual sites ...\n");
1619         fprintf(debug, "# # # VSITES # # #\n");
1620     }
1621
1622     ndb           = fflib_search_file_end(ffdir, ".vsd", FALSE, &db);
1623     nvsiteconf    = 0;
1624     vsiteconflist = NULL;
1625     nvsitetop     = 0;
1626     vsitetop      = NULL;
1627     for (f = 0; f < ndb; f++)
1628     {
1629         read_vsite_database(db[f], &vsiteconflist, &nvsiteconf, &vsitetop, &nvsitetop);
1630         sfree(db[f]);
1631     }
1632     sfree(db);
1633
1634     bFirstWater = TRUE;
1635     nvsite      = 0;
1636     nadd        = 0;
1637     /* we need a marker for which atoms should *not* be renumbered afterwards */
1638     add_shift = 10*at->nr;
1639     /* make arrays where masses can be inserted into */
1640     snew(newx, at->nr);
1641     snew(newatom, at->nr);
1642     snew(newatomname, at->nr);
1643     snew(newvsite_type, at->nr);
1644     snew(newcgnr, at->nr);
1645     /* make index array to tell where the atoms go to when masses are inserted */
1646     snew(o2n, at->nr);
1647     for (i = 0; i < at->nr; i++)
1648     {
1649         o2n[i] = i;
1650     }
1651     /* make index to tell which residues were already processed */
1652     snew(bResProcessed, at->nres);
1653
1654     gmx_residuetype_init(&rt);
1655
1656     /* generate vsite constructions */
1657     /* loop over all atoms */
1658     resind = -1;
1659     for (i = 0; (i < at->nr); i++)
1660     {
1661         if (at->atom[i].resind != resind)
1662         {
1663             resind = at->atom[i].resind;
1664             resnm  = *(at->resinfo[resind].name);
1665         }
1666         /* first check for aromatics to virtualize */
1667         /* don't waste our effort on DNA, water etc. */
1668         /* Only do the vsite aromatic stuff when we reach the
1669          * CA atom, since there might be an X2/X3 group on the
1670          * N-terminus that must be treated first.
1671          */
1672         if (bVsiteAromatics &&
1673             !strcmp(*(at->atomname[i]), "CA") &&
1674             !bResProcessed[resind] &&
1675             gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name)) )
1676         {
1677             /* mark this residue */
1678             bResProcessed[resind] = TRUE;
1679             /* find out if this residue needs converting */
1680             whatres = NOTSET;
1681             for (j = 0; j < resNR && whatres == NOTSET; j++)
1682             {
1683
1684                 cmplength = bPartial[j] ? strlen(resnm)-1 : strlen(resnm);
1685
1686                 bFound = ((gmx_strncasecmp(resnm, resnms[j], cmplength) == 0) ||
1687                           (gmx_strncasecmp(resnm, resnmsN[j], cmplength) == 0) ||
1688                           (gmx_strncasecmp(resnm, resnmsC[j], cmplength) == 0));
1689
1690                 if (bFound)
1691                 {
1692                     whatres = j;
1693                     /* get atoms we will be needing for the conversion */
1694                     nrfound = 0;
1695                     for (k = 0; atnms[j][k]; k++)
1696                     {
1697                         ats[k] = NOTSET;
1698                         for (m = i; m < at->nr && at->atom[m].resind == resind && ats[k] == NOTSET; m++)
1699                         {
1700                             if (gmx_strcasecmp(*(at->atomname[m]), atnms[j][k]) == 0)
1701                             {
1702                                 ats[k] = m;
1703                                 nrfound++;
1704                             }
1705                         }
1706                     }
1707
1708                     /* now k is number of atom names in atnms[j] */
1709                     if (j == resHIS)
1710                     {
1711                         needed = k-3;
1712                     }
1713                     else
1714                     {
1715                         needed = k;
1716                     }
1717                     if (nrfound < needed)
1718                     {
1719                         gmx_fatal(FARGS, "not enough atoms found (%d, need %d) in "
1720                                   "residue %s %d while\n             "
1721                                   "generating aromatics virtual site construction",
1722                                   nrfound, needed, resnm, at->resinfo[resind].nr);
1723                     }
1724                     /* Advance overall atom counter */
1725                     i++;
1726                 }
1727             }
1728             /* the enums for every residue MUST correspond to atnms[residue] */
1729             switch (whatres)
1730             {
1731                 case resPHE:
1732                     if (debug)
1733                     {
1734                         fprintf(stderr, "PHE at %d\n", o2n[ats[0]]+1);
1735                     }
1736                     nvsite += gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1737                     break;
1738                 case resTRP:
1739                     if (debug)
1740                     {
1741                         fprintf(stderr, "TRP at %d\n", o2n[ats[0]]+1);
1742                     }
1743                     nvsite += gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
1744                                              &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1745                                              at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1746                     break;
1747                 case resTYR:
1748                     if (debug)
1749                     {
1750                         fprintf(stderr, "TYR at %d\n", o2n[ats[0]]+1);
1751                     }
1752                     nvsite += gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
1753                                              &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1754                                              at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1755                     break;
1756                 case resHIS:
1757                     if (debug)
1758                     {
1759                         fprintf(stderr, "HIS at %d\n", o2n[ats[0]]+1);
1760                     }
1761                     nvsite += gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1762                     break;
1763                 case NOTSET:
1764                     /* this means this residue won't be processed */
1765                     break;
1766                 default:
1767                     gmx_fatal(FARGS, "DEATH HORROR in do_vsites (%s:%d)",
1768                               __FILE__, __LINE__);
1769             } /* switch whatres */
1770               /* skip back to beginning of residue */
1771             while (i > 0 && at->atom[i-1].resind == resind)
1772             {
1773                 i--;
1774             }
1775         } /* if bVsiteAromatics & is protein */
1776
1777         /* now process the rest of the hydrogens */
1778         /* only process hydrogen atoms which are not already set */
1779         if ( ((*vsite_type)[i] == NOTSET) && is_hydrogen(*(at->atomname[i])))
1780         {
1781             /* find heavy atom, count #bonds from it and #H atoms bound to it
1782                and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1783             count_bonds(i, &plist[F_BONDS], at->atomname,
1784                         &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
1785             /* get Heavy atom type */
1786             tpHeavy = get_atype(Heavy, at, nrtp, rtp, rt);
1787             strcpy(tpname, get_atomtype_name(tpHeavy, atype));
1788
1789             bWARNING       = FALSE;
1790             bAddVsiteParam = TRUE;
1791             /* nested if's which check nrHatoms, nrbonds and atomname */
1792             if (nrHatoms == 1)
1793             {
1794                 switch (nrbonds)
1795                 {
1796                     case 2: /* -O-H */
1797                         (*vsite_type)[i] = F_BONDS;
1798                         break;
1799                     case 3: /* =CH-, -NH- or =NH+- */
1800                         (*vsite_type)[i] = F_VSITE3FD;
1801                         break;
1802                     case 4: /* --CH- (tert) */
1803                         /* The old type 4FD had stability issues, so
1804                          * all new constructs should use 4FDN
1805                          */
1806                         (*vsite_type)[i] = F_VSITE4FDN;
1807
1808                         /* Check parity of heavy atoms from coordinates */
1809                         ai = Heavy;
1810                         aj = heavies[0];
1811                         ak = heavies[1];
1812                         al = heavies[2];
1813                         rvec_sub((*x)[aj], (*x)[ai], tmpmat[0]);
1814                         rvec_sub((*x)[ak], (*x)[ai], tmpmat[1]);
1815                         rvec_sub((*x)[al], (*x)[ai], tmpmat[2]);
1816
1817                         if (det(tmpmat) > 0)
1818                         {
1819                             /* swap parity */
1820                             heavies[1] = aj;
1821                             heavies[0] = ak;
1822                         }
1823
1824                         break;
1825                     default: /* nrbonds != 2, 3 or 4 */
1826                         bWARNING = TRUE;
1827                 }
1828
1829             }
1830             else if ( /*(nrHatoms == 2) && (nrbonds == 2) && REMOVED this test
1831                          DvdS 19-01-04 */
1832                 (gmx_strncasecmp(*at->atomname[Heavy], "OW", 2) == 0) )
1833             {
1834                 bAddVsiteParam = FALSE; /* this is water: skip these hydrogens */
1835                 if (bFirstWater)
1836                 {
1837                     bFirstWater = FALSE;
1838                     if (debug)
1839                     {
1840                         fprintf(debug,
1841                                 "Not converting hydrogens in water to virtual sites\n");
1842                     }
1843                 }
1844             }
1845             else if ( (nrHatoms == 2) && (nrbonds == 4) )
1846             {
1847                 /* -CH2- , -NH2+- */
1848                 (*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
1849                 (*vsite_type)[Hatoms[1]] = -F_VSITE3OUT;
1850             }
1851             else
1852             {
1853                 /* 2 or 3 hydrogen atom, with 3 or 4 bonds in total to the heavy atom.
1854                  * If it is a nitrogen, first check if it is planar.
1855                  */
1856                 isN = planarN = FALSE;
1857                 if ((nrHatoms == 2) && ((*at->atomname[Heavy])[0] == 'N'))
1858                 {
1859                     isN = TRUE;
1860                     j   = nitrogen_is_planar(vsiteconflist, nvsiteconf, tpname);
1861                     if (j < 0)
1862                     {
1863                         gmx_fatal(FARGS, "No vsite database NH2 entry for type %s\n", tpname);
1864                     }
1865                     planarN = (j == 1);
1866                 }
1867                 if ( (nrHatoms == 2) && (nrbonds == 3) && ( !isN || planarN ) )
1868                 {
1869                     /* =CH2 or, if it is a nitrogen NH2, it is a planar one */
1870                     (*vsite_type)[Hatoms[0]] = F_VSITE3FAD;
1871                     (*vsite_type)[Hatoms[1]] = -F_VSITE3FAD;
1872                 }
1873                 else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
1874                             ( isN && !planarN ) ) ||
1875                           ( (nrHatoms == 3) && (nrbonds == 4) ) )
1876                 {
1877                     /* CH3, NH3 or non-planar NH2 group */
1878                     int      Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
1879                     gmx_bool Hat_SwapParity[3] = { FALSE,    TRUE,        FALSE };
1880
1881                     if (debug)
1882                     {
1883                         fprintf(stderr, "-XH3 or nonplanar NH2 group at %d\n", i+1);
1884                     }
1885                     bAddVsiteParam = FALSE; /* we'll do this ourselves! */
1886                     /* -NH2 (umbrella), -NH3+ or -CH3 */
1887                     (*vsite_type)[Heavy]       = F_VSITE3;
1888                     for (j = 0; j < nrHatoms; j++)
1889                     {
1890                         (*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
1891                     }
1892                     /* get dummy mass type from first char of heavy atom type (N or C) */
1893
1894                     strcpy(nexttpname, get_atomtype_name(get_atype(heavies[0], at, nrtp, rtp, rt), atype));
1895                     ch = get_dummymass_name(vsiteconflist, nvsiteconf, tpname, nexttpname);
1896
1897                     if (ch == NULL)
1898                     {
1899                         if (ndb > 0)
1900                         {
1901                             gmx_fatal(FARGS, "Can't find dummy mass for type %s bonded to type %s in the virtual site database (.vsd files). Add it to the database!\n", tpname, nexttpname);
1902                         }
1903                         else
1904                         {
1905                             gmx_fatal(FARGS, "A dummy mass for type %s bonded to type %s is required, but no virtual site database (.vsd) files where found.\n", tpname, nexttpname);
1906                         }
1907                     }
1908                     else
1909                     {
1910                         strcpy(name, ch);
1911                     }
1912
1913                     tpM = vsite_nm2type(name, atype);
1914                     /* make space for 2 masses: shift all atoms starting with 'Heavy' */
1915 #define NMASS 2
1916                     i0  = Heavy;
1917                     ni0 = i0+nadd;
1918                     if (debug)
1919                     {
1920                         fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0]+1);
1921                     }
1922                     nadd += NMASS;
1923                     for (j = i0; j < at->nr; j++)
1924                     {
1925                         o2n[j] = j+nadd;
1926                     }
1927
1928                     srenew(newx, at->nr+nadd);
1929                     srenew(newatom, at->nr+nadd);
1930                     srenew(newatomname, at->nr+nadd);
1931                     srenew(newvsite_type, at->nr+nadd);
1932                     srenew(newcgnr, at->nr+nadd);
1933
1934                     for (j = 0; j < NMASS; j++)
1935                     {
1936                         newatomname[at->nr+nadd-1-j] = NULL;
1937                     }
1938
1939                     /* calculate starting position for the masses */
1940                     mHtot = 0;
1941                     /* get atom masses, and set Heavy and Hatoms mass to zero */
1942                     for (j = 0; j < nrHatoms; j++)
1943                     {
1944                         mHtot                += get_amass(Hatoms[j], at, nrtp, rtp, rt);
1945                         at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
1946                     }
1947                     mtot              = mHtot + get_amass(Heavy, at, nrtp, rtp, rt);
1948                     at->atom[Heavy].m = at->atom[Heavy].mB = 0;
1949                     if (mHmult != 1.0)
1950                     {
1951                         mHtot *= mHmult;
1952                     }
1953                     fact2 = mHtot/mtot;
1954                     fact  = sqrt(fact2);
1955                     /* generate vectors parallel and perpendicular to rotational axis:
1956                      * rpar  = Heavy -> Hcom
1957                      * rperp = Hcom  -> H1   */
1958                     clear_rvec(rpar);
1959                     for (j = 0; j < nrHatoms; j++)
1960                     {
1961                         rvec_inc(rpar, (*x)[Hatoms[j]]);
1962                     }
1963                     svmul(1.0/nrHatoms, rpar, rpar); /* rpar = ( H1+H2+H3 ) / 3 */
1964                     rvec_dec(rpar, (*x)[Heavy]);     /*        - Heavy          */
1965                     rvec_sub((*x)[Hatoms[0]], (*x)[Heavy], rperp);
1966                     rvec_dec(rperp, rpar);           /* rperp = H1 - Heavy - rpar */
1967                     /* calc mass positions */
1968                     svmul(fact2, rpar, temp);
1969                     for (j = 0; (j < NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
1970                     {
1971                         rvec_add((*x)[Heavy], temp, newx[ni0+j]);
1972                     }
1973                     svmul(fact, rperp, temp);
1974                     rvec_inc(newx[ni0  ], temp);
1975                     rvec_dec(newx[ni0+1], temp);
1976                     /* set atom parameters for the masses */
1977                     for (j = 0; (j < NMASS); j++)
1978                     {
1979                         /* make name: "M??#" or "M?#" (? is atomname, # is number) */
1980                         name[0] = 'M';
1981                         for (k = 0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++)
1982                         {
1983                             name[k+1] = (*at->atomname[Heavy])[k];
1984                         }
1985                         name[k+1]             = atomnamesuffix[j];
1986                         name[k+2]             = '\0';
1987                         newatomname[ni0+j]    = put_symtab(symtab, name);
1988                         newatom[ni0+j].m      = newatom[ni0+j].mB    = mtot/NMASS;
1989                         newatom[ni0+j].q      = newatom[ni0+j].qB    = 0.0;
1990                         newatom[ni0+j].type   = newatom[ni0+j].typeB = tpM;
1991                         newatom[ni0+j].ptype  = eptAtom;
1992                         newatom[ni0+j].resind = at->atom[i0].resind;
1993                         newvsite_type[ni0+j]  = NOTSET;
1994                         newcgnr[ni0+j]        = (*cgnr)[i0];
1995                     }
1996                     /* add constraints between dummy masses and to heavies[0] */
1997                     /* 'add_shift' says which atoms won't be renumbered afterwards */
1998                     my_add_param(&(plist[F_CONSTRNC]), heavies[0],  add_shift+ni0,  NOTSET);
1999                     my_add_param(&(plist[F_CONSTRNC]), heavies[0],  add_shift+ni0+1, NOTSET);
2000                     my_add_param(&(plist[F_CONSTRNC]), add_shift+ni0, add_shift+ni0+1, NOTSET);
2001
2002                     /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
2003                     /* note that vsite_type cannot be NOTSET, because we just set it */
2004                     add_vsite3_atoms  (&plist[(*vsite_type)[Heavy]],
2005                                        Heavy,     heavies[0], add_shift+ni0, add_shift+ni0+1,
2006                                        FALSE);
2007                     for (j = 0; j < nrHatoms; j++)
2008                     {
2009                         add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]],
2010                                          Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
2011                                          Hat_SwapParity[j]);
2012                     }
2013 #undef NMASS
2014                 }
2015                 else
2016                 {
2017                     bWARNING = TRUE;
2018                 }
2019
2020             }
2021             if (bWARNING)
2022             {
2023                 fprintf(stderr,
2024                         "Warning: cannot convert atom %d %s (bound to a heavy atom "
2025                         "%s with \n"
2026                         "         %d bonds and %d bound hydrogens atoms) to virtual site\n",
2027                         i+1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
2028             }
2029             if (bAddVsiteParam)
2030             {
2031                 /* add vsite parameters to topology,
2032                    also get rid of negative vsite_types */
2033                 add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms,
2034                            nrheavies, heavies);
2035                 /* transfer mass of virtual site to Heavy atom */
2036                 for (j = 0; j < nrHatoms; j++)
2037                 {
2038                     if (is_vsite((*vsite_type)[Hatoms[j]]))
2039                     {
2040                         at->atom[Heavy].m    += at->atom[Hatoms[j]].m;
2041                         at->atom[Heavy].mB    = at->atom[Heavy].m;
2042                         at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
2043                     }
2044                 }
2045             }
2046             nvsite += nrHatoms;
2047             if (debug)
2048             {
2049                 fprintf(debug, "atom %d: ", o2n[i]+1);
2050                 print_bonds(debug, o2n, nrHatoms, Hatoms, Heavy, nrheavies, heavies);
2051             }
2052         } /* if vsite NOTSET & is hydrogen */
2053
2054     }     /* for i < at->nr */
2055
2056     gmx_residuetype_destroy(rt);
2057
2058     if (debug)
2059     {
2060         fprintf(debug, "Before inserting new atoms:\n");
2061         for (i = 0; i < at->nr; i++)
2062         {
2063             fprintf(debug, "%4d %4d %4s %4d %4s %6d %-10s\n", i+1, o2n[i]+1,
2064                     at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2065                     at->resinfo[at->atom[i].resind].nr,
2066                     at->resinfo[at->atom[i].resind].name ?
2067                     *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2068                     (*cgnr)[i],
2069                     ((*vsite_type)[i] == NOTSET) ?
2070                     "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2071         }
2072         fprintf(debug, "new atoms to be inserted:\n");
2073         for (i = 0; i < at->nr+nadd; i++)
2074         {
2075             if (newatomname[i])
2076             {
2077                 fprintf(debug, "%4d %4s %4d %6d %-10s\n", i+1,
2078                         newatomname[i] ? *(newatomname[i]) : "(NULL)",
2079                         newatom[i].resind, newcgnr[i],
2080                         (newvsite_type[i] == NOTSET) ?
2081                         "NOTSET" : interaction_function[newvsite_type[i]].name);
2082             }
2083         }
2084     }
2085
2086     /* add all original atoms to the new arrays, using o2n index array */
2087     for (i = 0; i < at->nr; i++)
2088     {
2089         newatomname  [o2n[i]] = at->atomname [i];
2090         newatom      [o2n[i]] = at->atom     [i];
2091         newvsite_type[o2n[i]] = (*vsite_type)[i];
2092         newcgnr      [o2n[i]] = (*cgnr)      [i];
2093         copy_rvec((*x)[i], newx[o2n[i]]);
2094     }
2095     /* throw away old atoms */
2096     sfree(at->atom);
2097     sfree(at->atomname);
2098     sfree(*vsite_type);
2099     sfree(*cgnr);
2100     sfree(*x);
2101     /* put in the new ones */
2102     at->nr      += nadd;
2103     at->atom     = newatom;
2104     at->atomname = newatomname;
2105     *vsite_type  = newvsite_type;
2106     *cgnr        = newcgnr;
2107     *x           = newx;
2108     if (at->nr > add_shift)
2109     {
2110         gmx_fatal(FARGS, "Added impossible amount of dummy masses "
2111                   "(%d on a total of %d atoms)\n", nadd, at->nr-nadd);
2112     }
2113
2114     if (debug)
2115     {
2116         fprintf(debug, "After inserting new atoms:\n");
2117         for (i = 0; i < at->nr; i++)
2118         {
2119             fprintf(debug, "%4d %4s %4d %4s %6d %-10s\n", i+1,
2120                     at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2121                     at->resinfo[at->atom[i].resind].nr,
2122                     at->resinfo[at->atom[i].resind].name ?
2123                     *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2124                     (*cgnr)[i],
2125                     ((*vsite_type)[i] == NOTSET) ?
2126                     "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2127         }
2128     }
2129
2130     /* now renumber all the interactions because of the added atoms */
2131     for (ftype = 0; ftype < F_NRE; ftype++)
2132     {
2133         params = &(plist[ftype]);
2134         if (debug)
2135         {
2136             fprintf(debug, "Renumbering %d %s\n", params->nr,
2137                     interaction_function[ftype].longname);
2138         }
2139         for (i = 0; i < params->nr; i++)
2140         {
2141             for (j = 0; j < NRAL(ftype); j++)
2142             {
2143                 if (params->param[i].a[j] >= add_shift)
2144                 {
2145                     if (debug)
2146                     {
2147                         fprintf(debug, " [%u -> %u]", params->param[i].a[j],
2148                                 params->param[i].a[j]-add_shift);
2149                     }
2150                     params->param[i].a[j] = params->param[i].a[j]-add_shift;
2151                 }
2152                 else
2153                 {
2154                     if (debug)
2155                     {
2156                         fprintf(debug, " [%u -> %d]", params->param[i].a[j],
2157                                 o2n[params->param[i].a[j]]);
2158                     }
2159                     params->param[i].a[j] = o2n[params->param[i].a[j]];
2160                 }
2161             }
2162             if (debug)
2163             {
2164                 fprintf(debug, "\n");
2165             }
2166         }
2167     }
2168     /* now check if atoms in the added constraints are in increasing order */
2169     params = &(plist[F_CONSTRNC]);
2170     for (i = 0; i < params->nr; i++)
2171     {
2172         if (params->param[i].AI > params->param[i].AJ)
2173         {
2174             j                   = params->param[i].AJ;
2175             params->param[i].AJ = params->param[i].AI;
2176             params->param[i].AI = j;
2177         }
2178     }
2179
2180     /* clean up */
2181     sfree(o2n);
2182
2183     /* tell the user what we did */
2184     fprintf(stderr, "Marked %d virtual sites\n", nvsite);
2185     fprintf(stderr, "Added %d dummy masses\n", nadd);
2186     fprintf(stderr, "Added %d new constraints\n", plist[F_CONSTRNC].nr);
2187 }
2188
2189 void do_h_mass(t_params *psb, int vsite_type[], t_atoms *at, real mHmult,
2190                gmx_bool bDeuterate)
2191 {
2192     int i, j, a;
2193
2194     /* loop over all atoms */
2195     for (i = 0; i < at->nr; i++)
2196     {
2197         /* adjust masses if i is hydrogen and not a virtual site */
2198         if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])) )
2199         {
2200             /* find bonded heavy atom */
2201             a = NOTSET;
2202             for (j = 0; (j < psb->nr) && (a == NOTSET); j++)
2203             {
2204                 /* if other atom is not a virtual site, it is the one we want */
2205                 if ( (psb->param[j].AI == i) &&
2206                      !is_vsite(vsite_type[psb->param[j].AJ]) )
2207                 {
2208                     a = psb->param[j].AJ;
2209                 }
2210                 else if ( (psb->param[j].AJ == i) &&
2211                           !is_vsite(vsite_type[psb->param[j].AI]) )
2212                 {
2213                     a = psb->param[j].AI;
2214                 }
2215             }
2216             if (a == NOTSET)
2217             {
2218                 gmx_fatal(FARGS, "Unbound hydrogen atom (%d) found while adjusting mass",
2219                           i+1);
2220             }
2221
2222             /* adjust mass of i (hydrogen) with mHmult
2223                and correct mass of a (bonded atom) with same amount */
2224             if (!bDeuterate)
2225             {
2226                 at->atom[a].m  -= (mHmult-1.0)*at->atom[i].m;
2227                 at->atom[a].mB -= (mHmult-1.0)*at->atom[i].m;
2228             }
2229             at->atom[i].m  *= mHmult;
2230             at->atom[i].mB *= mHmult;
2231         }
2232     }
2233 }