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