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