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