src/gromacs/gmxpreprocess/gen_vsite.cpp

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