src/gromacs/gmxpreprocess/gen_vsite.cpp

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