src/gromacs/gmxpreprocess/gen_vsite.cpp

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