src/gromacs/gmxpreprocess/gen_vsite.c

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