src/kernel/gen_vsite.c

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