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