src/gromacs/gmxpreprocess/gen_vsite.cpp

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