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