src/gromacs/mdlib/vsite.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <stdio.h>
  42
  43 #include "typedefs.h"
  44 #include "types/commrec.h"
  45 #include "vsite.h"
  46 #include "macros.h"
  47 #include "gromacs/utility/smalloc.h"
  48 #include "nrnb.h"
  49 #include "gromacs/math/vec.h"
  50 #include "network.h"
  51 #include "gromacs/pbcutil/mshift.h"
  52 #include "pbc.h"
  53 #include "domdec.h"
  54 #include "mtop_util.h"
  55 #include "gmx_omp_nthreads.h"
  56
  57 #include "gromacs/utility/gmxomp.h"
  58
  59 /* Routines to send/recieve coordinates and force
  60  * of constructing atoms.
  61  */
  62
  63 static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx)
  64 {
  65     if (pbc)
  66     {
  67         return pbc_dx_aiuc(pbc, xi, xj, dx);
  68     }
  69     else
  70     {
  71         rvec_sub(xi, xj, dx);
  72         return CENTRAL;
  73     }
  74 }
  75
  76 /* Vsite construction routines */
  77
  78 static void constr_vsite2(rvec xi, rvec xj, rvec x, real a, t_pbc *pbc)
  79 {
  80     real b;
  81     rvec dx;
  82
  83     b = 1.0-a;
  84     /* 1 flop */
  85
  86     if (pbc)
  87     {
  88         pbc_dx_aiuc(pbc, xj, xi, dx);
  89         x[XX] = xi[XX] + a*dx[XX];
  90         x[YY] = xi[YY] + a*dx[YY];
  91         x[ZZ] = xi[ZZ] + a*dx[ZZ];
  92     }
  93     else
  94     {
  95         x[XX] = b*xi[XX] + a*xj[XX];
  96         x[YY] = b*xi[YY] + a*xj[YY];
  97         x[ZZ] = b*xi[ZZ] + a*xj[ZZ];
  98         /* 9 Flops */
  99     }
 100
 101     /* TOTAL: 10 flops */
 102 }
 103
 104 static void constr_vsite3(rvec xi, rvec xj, rvec xk, rvec x, real a, real b,
 105                           t_pbc *pbc)
 106 {
 107     real c;
 108     rvec dxj, dxk;
 109
 110     c = 1.0-a-b;
 111     /* 2 flops */
 112
 113     if (pbc)
 114     {
 115         pbc_dx_aiuc(pbc, xj, xi, dxj);
 116         pbc_dx_aiuc(pbc, xk, xi, dxk);
 117         x[XX] = xi[XX] + a*dxj[XX] + b*dxk[XX];
 118         x[YY] = xi[YY] + a*dxj[YY] + b*dxk[YY];
 119         x[ZZ] = xi[ZZ] + a*dxj[ZZ] + b*dxk[ZZ];
 120     }
 121     else
 122     {
 123         x[XX] = c*xi[XX] + a*xj[XX] + b*xk[XX];
 124         x[YY] = c*xi[YY] + a*xj[YY] + b*xk[YY];
 125         x[ZZ] = c*xi[ZZ] + a*xj[ZZ] + b*xk[ZZ];
 126         /* 15 Flops */
 127     }
 128
 129     /* TOTAL: 17 flops */
 130 }
 131
 132 static void constr_vsite3FD(rvec xi, rvec xj, rvec xk, rvec x, real a, real b,
 133                             t_pbc *pbc)
 134 {
 135     rvec xij, xjk, temp;
 136     real c;
 137
 138     pbc_rvec_sub(pbc, xj, xi, xij);
 139     pbc_rvec_sub(pbc, xk, xj, xjk);
 140     /* 6 flops */
 141
 142     /* temp goes from i to a point on the line jk */
 143     temp[XX] = xij[XX] + a*xjk[XX];
 144     temp[YY] = xij[YY] + a*xjk[YY];
 145     temp[ZZ] = xij[ZZ] + a*xjk[ZZ];
 146     /* 6 flops */
 147
 148     c = b*gmx_invsqrt(iprod(temp, temp));
 149     /* 6 + 10 flops */
 150
 151     x[XX] = xi[XX] + c*temp[XX];
 152     x[YY] = xi[YY] + c*temp[YY];
 153     x[ZZ] = xi[ZZ] + c*temp[ZZ];
 154     /* 6 Flops */
 155
 156     /* TOTAL: 34 flops */
 157 }
 158
 159 static void constr_vsite3FAD(rvec xi, rvec xj, rvec xk, rvec x, real a, real b, t_pbc *pbc)
 160 {
 161     rvec xij, xjk, xp;
 162     real a1, b1, c1, invdij;
 163
 164     pbc_rvec_sub(pbc, xj, xi, xij);
 165     pbc_rvec_sub(pbc, xk, xj, xjk);
 166     /* 6 flops */
 167
 168     invdij = gmx_invsqrt(iprod(xij, xij));
 169     c1     = invdij * invdij * iprod(xij, xjk);
 170     xp[XX] = xjk[XX] - c1*xij[XX];
 171     xp[YY] = xjk[YY] - c1*xij[YY];
 172     xp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
 173     a1     = a*invdij;
 174     b1     = b*gmx_invsqrt(iprod(xp, xp));
 175     /* 45 */
 176
 177     x[XX] = xi[XX] + a1*xij[XX] + b1*xp[XX];
 178     x[YY] = xi[YY] + a1*xij[YY] + b1*xp[YY];
 179     x[ZZ] = xi[ZZ] + a1*xij[ZZ] + b1*xp[ZZ];
 180     /* 12 Flops */
 181
 182     /* TOTAL: 63 flops */
 183 }
 184
 185 static void constr_vsite3OUT(rvec xi, rvec xj, rvec xk, rvec x,
 186                              real a, real b, real c, t_pbc *pbc)
 187 {
 188     rvec xij, xik, temp;
 189
 190     pbc_rvec_sub(pbc, xj, xi, xij);
 191     pbc_rvec_sub(pbc, xk, xi, xik);
 192     cprod(xij, xik, temp);
 193     /* 15 Flops */
 194
 195     x[XX] = xi[XX] + a*xij[XX] + b*xik[XX] + c*temp[XX];
 196     x[YY] = xi[YY] + a*xij[YY] + b*xik[YY] + c*temp[YY];
 197     x[ZZ] = xi[ZZ] + a*xij[ZZ] + b*xik[ZZ] + c*temp[ZZ];
 198     /* 18 Flops */
 199
 200     /* TOTAL: 33 flops */
 201 }
 202
 203 static void constr_vsite4FD(rvec xi, rvec xj, rvec xk, rvec xl, rvec x,
 204                             real a, real b, real c, t_pbc *pbc)
 205 {
 206     rvec xij, xjk, xjl, temp;
 207     real d;
 208
 209     pbc_rvec_sub(pbc, xj, xi, xij);
 210     pbc_rvec_sub(pbc, xk, xj, xjk);
 211     pbc_rvec_sub(pbc, xl, xj, xjl);
 212     /* 9 flops */
 213
 214     /* temp goes from i to a point on the plane jkl */
 215     temp[XX] = xij[XX] + a*xjk[XX] + b*xjl[XX];
 216     temp[YY] = xij[YY] + a*xjk[YY] + b*xjl[YY];
 217     temp[ZZ] = xij[ZZ] + a*xjk[ZZ] + b*xjl[ZZ];
 218     /* 12 flops */
 219
 220     d = c*gmx_invsqrt(iprod(temp, temp));
 221     /* 6 + 10 flops */
 222
 223     x[XX] = xi[XX] + d*temp[XX];
 224     x[YY] = xi[YY] + d*temp[YY];
 225     x[ZZ] = xi[ZZ] + d*temp[ZZ];
 226     /* 6 Flops */
 227
 228     /* TOTAL: 43 flops */
 229 }
 230
 231
 232 static void constr_vsite4FDN(rvec xi, rvec xj, rvec xk, rvec xl, rvec x,
 233                              real a, real b, real c, t_pbc *pbc)
 234 {
 235     rvec xij, xik, xil, ra, rb, rja, rjb, rm;
 236     real d;
 237
 238     pbc_rvec_sub(pbc, xj, xi, xij);
 239     pbc_rvec_sub(pbc, xk, xi, xik);
 240     pbc_rvec_sub(pbc, xl, xi, xil);
 241     /* 9 flops */
 242
 243     ra[XX] = a*xik[XX];
 244     ra[YY] = a*xik[YY];
 245     ra[ZZ] = a*xik[ZZ];
 246
 247     rb[XX] = b*xil[XX];
 248     rb[YY] = b*xil[YY];
 249     rb[ZZ] = b*xil[ZZ];
 250
 251     /* 6 flops */
 252
 253     rvec_sub(ra, xij, rja);
 254     rvec_sub(rb, xij, rjb);
 255     /* 6 flops */
 256
 257     cprod(rja, rjb, rm);
 258     /* 9 flops */
 259
 260     d = c*gmx_invsqrt(norm2(rm));
 261     /* 5+5+1 flops */
 262
 263     x[XX] = xi[XX] + d*rm[XX];
 264     x[YY] = xi[YY] + d*rm[YY];
 265     x[ZZ] = xi[ZZ] + d*rm[ZZ];
 266     /* 6 Flops */
 267
 268     /* TOTAL: 47 flops */
 269 }
 270
 271
 272 static int constr_vsiten(t_iatom *ia, t_iparams ip[],
 273                          rvec *x, t_pbc *pbc)
 274 {
 275     rvec xs, x1, dx;
 276     dvec dsum;
 277     int  n3, av, ai, i;
 278     real a;
 279
 280     n3 = 3*ip[ia[0]].vsiten.n;
 281     av = ia[1];
 282     ai = ia[2];
 283     copy_rvec(x[ai], x1);
 284     clear_dvec(dsum);
 285     for (i = 3; i < n3; i += 3)
 286     {
 287         ai = ia[i+2];
 288         a  = ip[ia[i]].vsiten.a;
 289         if (pbc)
 290         {
 291             pbc_dx_aiuc(pbc, x[ai], x1, dx);
 292         }
 293         else
 294         {
 295             rvec_sub(x[ai], x1, dx);
 296         }
 297         dsum[XX] += a*dx[XX];
 298         dsum[YY] += a*dx[YY];
 299         dsum[ZZ] += a*dx[ZZ];
 300         /* 9 Flops */
 301     }
 302
 303     x[av][XX] = x1[XX] + dsum[XX];
 304     x[av][YY] = x1[YY] + dsum[YY];
 305     x[av][ZZ] = x1[ZZ] + dsum[ZZ];
 306
 307     return n3;
 308 }
 309
 310
 311 void construct_vsites_thread(gmx_vsite_t *vsite,
 312                              rvec x[],
 313                              real dt, rvec *v,
 314                              t_iparams ip[], t_ilist ilist[],
 315                              t_pbc *pbc_null)
 316 {
 317     gmx_bool   bPBCAll;
 318     rvec       xpbc, xv, vv, dx;
 319     real       a1, b1, c1, inv_dt;
 320     int        i, inc, ii, nra, nr, tp, ftype;
 321     t_iatom    avsite, ai, aj, ak, al, pbc_atom;
 322     t_iatom   *ia;
 323     t_pbc     *pbc_null2;
 324     int       *vsite_pbc, ishift;
 325     rvec       reftmp, vtmp, rtmp;
 326
 327     if (v != NULL)
 328     {
 329         inv_dt = 1.0/dt;
 330     }
 331     else
 332     {
 333         inv_dt = 1.0;
 334     }
 335
 336     bPBCAll = (pbc_null != NULL && !vsite->bHaveChargeGroups);
 337
 338     pbc_null2 = NULL;
 339     vsite_pbc = NULL;
 340     for (ftype = 0; (ftype < F_NRE); ftype++)
 341     {
 342         if ((interaction_function[ftype].flags & IF_VSITE) &&
 343             ilist[ftype].nr > 0)
 344         {
 345             nra    = interaction_function[ftype].nratoms;
 346             inc    = 1 + nra;
 347             nr     = ilist[ftype].nr;
 348             ia     = ilist[ftype].iatoms;
 349
 350             if (bPBCAll)
 351             {
 352                 pbc_null2 = pbc_null;
 353             }
 354             else if (pbc_null != NULL)
 355             {
 356                 vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2];
 357             }
 358
 359             for (i = 0; i < nr; )
 360             {
 361                 tp   = ia[0];
 362
 363                 /* The vsite and constructing atoms */
 364                 avsite = ia[1];
 365                 ai     = ia[2];
 366
 367                 /* Constants for constructing vsites */
 368                 a1   = ip[tp].vsite.a;
 369                 /* Check what kind of pbc we need to use */
 370                 if (bPBCAll)
 371                 {
 372                     /* No charge groups, vsite follows its own pbc */
 373                     pbc_atom = avsite;
 374                     copy_rvec(x[avsite], xpbc);
 375                 }
 376                 else if (vsite_pbc != NULL)
 377                 {
 378                     pbc_atom = vsite_pbc[i/(1+nra)];
 379                     if (pbc_atom > -2)
 380                     {
 381                         if (pbc_atom >= 0)
 382                         {
 383                             /* We need to copy the coordinates here,
 384                              * single for single atom cg's pbc_atom
 385                              * is the vsite itself.
 386                              */
 387                             copy_rvec(x[pbc_atom], xpbc);
 388                         }
 389                         pbc_null2 = pbc_null;
 390                     }
 391                     else
 392                     {
 393                         pbc_null2 = NULL;
 394                     }
 395                 }
 396                 else
 397                 {
 398                     pbc_atom = -2;
 399                 }
 400                 /* Copy the old position */
 401                 copy_rvec(x[avsite], xv);
 402
 403                 /* Construct the vsite depending on type */
 404                 switch (ftype)
 405                 {
 406                     case F_VSITE2:
 407                         aj = ia[3];
 408                         constr_vsite2(x[ai], x[aj], x[avsite], a1, pbc_null2);
 409                         break;
 410                     case F_VSITE3:
 411                         aj = ia[3];
 412                         ak = ia[4];
 413                         b1 = ip[tp].vsite.b;
 414                         constr_vsite3(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2);
 415                         break;
 416                     case F_VSITE3FD:
 417                         aj = ia[3];
 418                         ak = ia[4];
 419                         b1 = ip[tp].vsite.b;
 420                         constr_vsite3FD(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2);
 421                         break;
 422                     case F_VSITE3FAD:
 423                         aj = ia[3];
 424                         ak = ia[4];
 425                         b1 = ip[tp].vsite.b;
 426                         constr_vsite3FAD(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2);
 427                         break;
 428                     case F_VSITE3OUT:
 429                         aj = ia[3];
 430                         ak = ia[4];
 431                         b1 = ip[tp].vsite.b;
 432                         c1 = ip[tp].vsite.c;
 433                         constr_vsite3OUT(x[ai], x[aj], x[ak], x[avsite], a1, b1, c1, pbc_null2);
 434                         break;
 435                     case F_VSITE4FD:
 436                         aj = ia[3];
 437                         ak = ia[4];
 438                         al = ia[5];
 439                         b1 = ip[tp].vsite.b;
 440                         c1 = ip[tp].vsite.c;
 441                         constr_vsite4FD(x[ai], x[aj], x[ak], x[al], x[avsite], a1, b1, c1,
 442                                         pbc_null2);
 443                         break;
 444                     case F_VSITE4FDN:
 445                         aj = ia[3];
 446                         ak = ia[4];
 447                         al = ia[5];
 448                         b1 = ip[tp].vsite.b;
 449                         c1 = ip[tp].vsite.c;
 450                         constr_vsite4FDN(x[ai], x[aj], x[ak], x[al], x[avsite], a1, b1, c1,
 451                                          pbc_null2);
 452                         break;
 453                     case F_VSITEN:
 454                         inc = constr_vsiten(ia, ip, x, pbc_null2);
 455                         break;
 456                     default:
 457                         gmx_fatal(FARGS, "No such vsite type %d in %s, line %d",
 458                                   ftype, __FILE__, __LINE__);
 459                 }
 460
 461                 if (pbc_atom >= 0)
 462                 {
 463                     /* Match the pbc of this vsite to the rest of its charge group */
 464                     ishift = pbc_dx_aiuc(pbc_null, x[avsite], xpbc, dx);
 465                     if (ishift != CENTRAL)
 466                     {
 467                         rvec_add(xpbc, dx, x[avsite]);
 468                     }
 469                 }
 470                 if (v != NULL)
 471                 {
 472                     /* Calculate velocity of vsite... */
 473                     rvec_sub(x[avsite], xv, vv);
 474                     svmul(inv_dt, vv, v[avsite]);
 475                 }
 476
 477                 /* Increment loop variables */
 478                 i  += inc;
 479                 ia += inc;
 480             }
 481         }
 482     }
 483 }
 484
 485 void construct_vsites(gmx_vsite_t *vsite,
 486                       rvec x[],
 487                       real dt, rvec *v,
 488                       t_iparams ip[], t_ilist ilist[],
 489                       int ePBC, gmx_bool bMolPBC,
 490                       t_commrec *cr, matrix box)
 491 {
 492     t_pbc     pbc, *pbc_null;
 493     gmx_bool  bDomDec;
 494     int       nthreads;
 495
 496     bDomDec = cr && DOMAINDECOMP(cr);
 497
 498     /* We only need to do pbc when we have inter-cg vsites */
 499     if (ePBC != epbcNONE && (bDomDec || bMolPBC) && vsite->n_intercg_vsite)
 500     {
 501         /* This is wasting some CPU time as we now do this multiple times
 502          * per MD step. But how often do we have vsites with full pbc?
 503          */
 504         pbc_null = set_pbc_dd(&pbc, ePBC, cr != NULL ? cr->dd : NULL, FALSE, box);
 505     }
 506     else
 507     {
 508         pbc_null = NULL;
 509     }
 510
 511     if (bDomDec)
 512     {
 513         dd_move_x_vsites(cr->dd, box, x);
 514     }
 515
 516     if (vsite->nthreads == 1)
 517     {
 518         construct_vsites_thread(vsite,
 519                                 x, dt, v,
 520                                 ip, ilist,
 521                                 pbc_null);
 522     }
 523     else
 524     {
 525 #pragma omp parallel num_threads(vsite->nthreads)
 526         {
 527             construct_vsites_thread(vsite,
 528                                     x, dt, v,
 529                                     ip, vsite->tdata[gmx_omp_get_thread_num()].ilist,
 530                                     pbc_null);
 531         }
 532         /* Now we can construct the vsites that might depend on other vsites */
 533         construct_vsites_thread(vsite,
 534                                 x, dt, v,
 535                                 ip, vsite->tdata[vsite->nthreads].ilist,
 536                                 pbc_null);
 537     }
 538 }
 539
 540 static void spread_vsite2(t_iatom ia[], real a,
 541                           rvec x[], rvec f[], rvec fshift[],
 542                           t_pbc *pbc, t_graph *g)
 543 {
 544     rvec    fi, fj, dx;
 545     t_iatom av, ai, aj;
 546     ivec    di;
 547     real    b;
 548     int     siv, sij;
 549
 550     av = ia[1];
 551     ai = ia[2];
 552     aj = ia[3];
 553
 554     svmul(1-a, f[av], fi);
 555     svmul(  a, f[av], fj);
 556     /* 7 flop */
 557
 558     rvec_inc(f[ai], fi);
 559     rvec_inc(f[aj], fj);
 560     /* 6 Flops */
 561
 562     if (g)
 563     {
 564         ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, av), di);
 565         siv = IVEC2IS(di);
 566         ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), di);
 567         sij = IVEC2IS(di);
 568     }
 569     else if (pbc)
 570     {
 571         siv = pbc_dx_aiuc(pbc, x[ai], x[av], dx);
 572         sij = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
 573     }
 574     else
 575     {
 576         siv = CENTRAL;
 577         sij = CENTRAL;
 578     }
 579
 580     if (fshift && (siv != CENTRAL || sij != CENTRAL))
 581     {
 582         rvec_inc(fshift[siv], f[av]);
 583         rvec_dec(fshift[CENTRAL], fi);
 584         rvec_dec(fshift[sij], fj);
 585     }
 586
 587     /* TOTAL: 13 flops */
 588 }
 589
 590 void construct_vsites_mtop(gmx_vsite_t *vsite,
 591                            gmx_mtop_t *mtop, rvec x[])
 592 {
 593     int             as, mb, mol;
 594     gmx_molblock_t *molb;
 595     gmx_moltype_t  *molt;
 596
 597     as = 0;
 598     for (mb = 0; mb < mtop->nmolblock; mb++)
 599     {
 600         molb = &mtop->molblock[mb];
 601         molt = &mtop->moltype[molb->type];
 602         for (mol = 0; mol < molb->nmol; mol++)
 603         {
 604             construct_vsites(vsite, x+as, 0.0, NULL,
 605                              mtop->ffparams.iparams, molt->ilist,
 606                              epbcNONE, TRUE, NULL, NULL);
 607             as += molt->atoms.nr;
 608         }
 609     }
 610 }
 611
 612 static void spread_vsite3(t_iatom ia[], real a, real b,
 613                           rvec x[], rvec f[], rvec fshift[],
 614                           t_pbc *pbc, t_graph *g)
 615 {
 616     rvec    fi, fj, fk, dx;
 617     atom_id av, ai, aj, ak;
 618     ivec    di;
 619     int     siv, sij, sik;
 620
 621     av = ia[1];
 622     ai = ia[2];
 623     aj = ia[3];
 624     ak = ia[4];
 625
 626     svmul(1-a-b, f[av], fi);
 627     svmul(    a, f[av], fj);
 628     svmul(    b, f[av], fk);
 629     /* 11 flops */
 630
 631     rvec_inc(f[ai], fi);
 632     rvec_inc(f[aj], fj);
 633     rvec_inc(f[ak], fk);
 634     /* 9 Flops */
 635
 636     if (g)
 637     {
 638         ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, ia[1]), di);
 639         siv = IVEC2IS(di);
 640         ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), di);
 641         sij = IVEC2IS(di);
 642         ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, ak), di);
 643         sik = IVEC2IS(di);
 644     }
 645     else if (pbc)
 646     {
 647         siv = pbc_dx_aiuc(pbc, x[ai], x[av], dx);
 648         sij = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
 649         sik = pbc_dx_aiuc(pbc, x[ai], x[ak], dx);
 650     }
 651     else
 652     {
 653         siv = CENTRAL;
 654         sij = CENTRAL;
 655         sik = CENTRAL;
 656     }
 657
 658     if (fshift && (siv != CENTRAL || sij != CENTRAL || sik != CENTRAL))
 659     {
 660         rvec_inc(fshift[siv], f[av]);
 661         rvec_dec(fshift[CENTRAL], fi);
 662         rvec_dec(fshift[sij], fj);
 663         rvec_dec(fshift[sik], fk);
 664     }
 665
 666     /* TOTAL: 20 flops */
 667 }
 668
 669 static void spread_vsite3FD(t_iatom ia[], real a, real b,
 670                             rvec x[], rvec f[], rvec fshift[],
 671                             gmx_bool VirCorr, matrix dxdf,
 672                             t_pbc *pbc, t_graph *g)
 673 {
 674     real    fx, fy, fz, c, invl, fproj, a1;
 675     rvec    xvi, xij, xjk, xix, fv, temp;
 676     t_iatom av, ai, aj, ak;
 677     int     svi, sji, skj, d;
 678     ivec    di;
 679
 680     av = ia[1];
 681     ai = ia[2];
 682     aj = ia[3];
 683     ak = ia[4];
 684     copy_rvec(f[av], fv);
 685
 686     sji = pbc_rvec_sub(pbc, x[aj], x[ai], xij);
 687     skj = pbc_rvec_sub(pbc, x[ak], x[aj], xjk);
 688     /* 6 flops */
 689
 690     /* xix goes from i to point x on the line jk */
 691     xix[XX] = xij[XX]+a*xjk[XX];
 692     xix[YY] = xij[YY]+a*xjk[YY];
 693     xix[ZZ] = xij[ZZ]+a*xjk[ZZ];
 694     /* 6 flops */
 695
 696     invl = gmx_invsqrt(iprod(xix, xix));
 697     c    = b*invl;
 698     /* 4 + ?10? flops */
 699
 700     fproj = iprod(xix, fv)*invl*invl; /* = (xix . f)/(xix . xix) */
 701
 702     temp[XX] = c*(fv[XX]-fproj*xix[XX]);
 703     temp[YY] = c*(fv[YY]-fproj*xix[YY]);
 704     temp[ZZ] = c*(fv[ZZ]-fproj*xix[ZZ]);
 705     /* 16 */
 706
 707     /* c is already calculated in constr_vsite3FD
 708        storing c somewhere will save 26 flops!     */
 709
 710     a1         = 1-a;
 711     f[ai][XX] += fv[XX] - temp[XX];
 712     f[ai][YY] += fv[YY] - temp[YY];
 713     f[ai][ZZ] += fv[ZZ] - temp[ZZ];
 714     f[aj][XX] += a1*temp[XX];
 715     f[aj][YY] += a1*temp[YY];
 716     f[aj][ZZ] += a1*temp[ZZ];
 717     f[ak][XX] += a*temp[XX];
 718     f[ak][YY] += a*temp[YY];
 719     f[ak][ZZ] += a*temp[ZZ];
 720     /* 19 Flops */
 721
 722     if (g)
 723     {
 724         ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di);
 725         svi = IVEC2IS(di);
 726         ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di);
 727         sji = IVEC2IS(di);
 728         ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, aj), di);
 729         skj = IVEC2IS(di);
 730     }
 731     else if (pbc)
 732     {
 733         svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi);
 734     }
 735     else
 736     {
 737         svi = CENTRAL;
 738     }
 739
 740     if (fshift && (svi != CENTRAL || sji != CENTRAL || skj != CENTRAL))
 741     {
 742         rvec_dec(fshift[svi], fv);
 743         fshift[CENTRAL][XX] += fv[XX] - (1 + a)*temp[XX];
 744         fshift[CENTRAL][YY] += fv[YY] - (1 + a)*temp[YY];
 745         fshift[CENTRAL][ZZ] += fv[ZZ] - (1 + a)*temp[ZZ];
 746         fshift[    sji][XX] += temp[XX];
 747         fshift[    sji][YY] += temp[YY];
 748         fshift[    sji][ZZ] += temp[ZZ];
 749         fshift[    skj][XX] += a*temp[XX];
 750         fshift[    skj][YY] += a*temp[YY];
 751         fshift[    skj][ZZ] += a*temp[ZZ];
 752     }
 753
 754     if (VirCorr)
 755     {
 756         /* When VirCorr=TRUE, the virial for the current forces is not
 757          * calculated from the redistributed forces. This means that
 758          * the effect of non-linear virtual site constructions on the virial
 759          * needs to be added separately. This contribution can be calculated
 760          * in many ways, but the simplest and cheapest way is to use
 761          * the first constructing atom ai as a reference position in space:
 762          * subtract (xv-xi)*fv and add (xj-xi)*fj + (xk-xi)*fk.
 763          */
 764         rvec xiv;
 765         int  i, j;
 766
 767         pbc_rvec_sub(pbc, x[av], x[ai], xiv);
 768
 769         for (i = 0; i < DIM; i++)
 770         {
 771             for (j = 0; j < DIM; j++)
 772             {
 773                 /* As xix is a linear combination of j and k, use that here */
 774                 dxdf[i][j] += -xiv[i]*fv[j] + xix[i]*temp[j];
 775             }
 776         }
 777     }
 778
 779     /* TOTAL: 61 flops */
 780 }
 781
 782 static void spread_vsite3FAD(t_iatom ia[], real a, real b,
 783                              rvec x[], rvec f[], rvec fshift[],
 784                              gmx_bool VirCorr, matrix dxdf,
 785                              t_pbc *pbc, t_graph *g)
 786 {
 787     rvec    xvi, xij, xjk, xperp, Fpij, Fppp, fv, f1, f2, f3;
 788     real    a1, b1, c1, c2, invdij, invdij2, invdp, fproj;
 789     t_iatom av, ai, aj, ak;
 790     int     svi, sji, skj, d;
 791     ivec    di;
 792
 793     av = ia[1];
 794     ai = ia[2];
 795     aj = ia[3];
 796     ak = ia[4];
 797     copy_rvec(f[ia[1]], fv);
 798
 799     sji = pbc_rvec_sub(pbc, x[aj], x[ai], xij);
 800     skj = pbc_rvec_sub(pbc, x[ak], x[aj], xjk);
 801     /* 6 flops */
 802
 803     invdij    = gmx_invsqrt(iprod(xij, xij));
 804     invdij2   = invdij * invdij;
 805     c1        = iprod(xij, xjk) * invdij2;
 806     xperp[XX] = xjk[XX] - c1*xij[XX];
 807     xperp[YY] = xjk[YY] - c1*xij[YY];
 808     xperp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
 809     /* xperp in plane ijk, perp. to ij */
 810     invdp = gmx_invsqrt(iprod(xperp, xperp));
 811     a1    = a*invdij;
 812     b1    = b*invdp;
 813     /* 45 flops */
 814
 815     /* a1, b1 and c1 are already calculated in constr_vsite3FAD
 816        storing them somewhere will save 45 flops!     */
 817
 818     fproj = iprod(xij, fv)*invdij2;
 819     svmul(fproj,                      xij,  Fpij);    /* proj. f on xij */
 820     svmul(iprod(xperp, fv)*invdp*invdp, xperp, Fppp); /* proj. f on xperp */
 821     svmul(b1*fproj,                   xperp, f3);
 822     /* 23 flops */
 823
 824     rvec_sub(fv, Fpij, f1); /* f1 = f - Fpij */
 825     rvec_sub(f1, Fppp, f2); /* f2 = f - Fpij - Fppp */
 826     for (d = 0; (d < DIM); d++)
 827     {
 828         f1[d] *= a1;
 829         f2[d] *= b1;
 830     }
 831     /* 12 flops */
 832
 833     c2         = 1+c1;
 834     f[ai][XX] += fv[XX] - f1[XX] + c1*f2[XX] + f3[XX];
 835     f[ai][YY] += fv[YY] - f1[YY] + c1*f2[YY] + f3[YY];
 836     f[ai][ZZ] += fv[ZZ] - f1[ZZ] + c1*f2[ZZ] + f3[ZZ];
 837     f[aj][XX] +=          f1[XX] - c2*f2[XX] - f3[XX];
 838     f[aj][YY] +=          f1[YY] - c2*f2[YY] - f3[YY];
 839     f[aj][ZZ] +=          f1[ZZ] - c2*f2[ZZ] - f3[ZZ];
 840     f[ak][XX] +=                      f2[XX];
 841     f[ak][YY] +=                      f2[YY];
 842     f[ak][ZZ] +=                      f2[ZZ];
 843     /* 30 Flops */
 844
 845     if (g)
 846     {
 847         ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di);
 848         svi = IVEC2IS(di);
 849         ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di);
 850         sji = IVEC2IS(di);
 851         ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, aj), di);
 852         skj = IVEC2IS(di);
 853     }
 854     else if (pbc)
 855     {
 856         svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi);
 857     }
 858     else
 859     {
 860         svi = CENTRAL;
 861     }
 862
 863     if (fshift && (svi != CENTRAL || sji != CENTRAL || skj != CENTRAL))
 864     {
 865         rvec_dec(fshift[svi], fv);
 866         fshift[CENTRAL][XX] += fv[XX] - f1[XX] - (1-c1)*f2[XX] + f3[XX];
 867         fshift[CENTRAL][YY] += fv[YY] - f1[YY] - (1-c1)*f2[YY] + f3[YY];
 868         fshift[CENTRAL][ZZ] += fv[ZZ] - f1[ZZ] - (1-c1)*f2[ZZ] + f3[ZZ];
 869         fshift[    sji][XX] +=          f1[XX] -    c1 *f2[XX] - f3[XX];
 870         fshift[    sji][YY] +=          f1[YY] -    c1 *f2[YY] - f3[YY];
 871         fshift[    sji][ZZ] +=          f1[ZZ] -    c1 *f2[ZZ] - f3[ZZ];
 872         fshift[    skj][XX] +=                          f2[XX];
 873         fshift[    skj][YY] +=                          f2[YY];
 874         fshift[    skj][ZZ] +=                          f2[ZZ];
 875     }
 876
 877     if (VirCorr)
 878     {
 879         rvec xiv;
 880         int  i, j;
 881
 882         pbc_rvec_sub(pbc, x[av], x[ai], xiv);
 883
 884         for (i = 0; i < DIM; i++)
 885         {
 886             for (j = 0; j < DIM; j++)
 887             {
 888                 /* Note that xik=xij+xjk, so we have to add xij*f2 */
 889                 dxdf[i][j] +=
 890                     -xiv[i]*fv[j]
 891                     + xij[i]*(f1[j] + (1 - c2)*f2[j] - f3[j])
 892                     + xjk[i]*f2[j];
 893             }
 894         }
 895     }
 896
 897     /* TOTAL: 113 flops */
 898 }
 899
 900 static void spread_vsite3OUT(t_iatom ia[], real a, real b, real c,
 901                              rvec x[], rvec f[], rvec fshift[],
 902                              gmx_bool VirCorr, matrix dxdf,
 903                              t_pbc *pbc, t_graph *g)
 904 {
 905     rvec    xvi, xij, xik, fv, fj, fk;
 906     real    cfx, cfy, cfz;
 907     atom_id av, ai, aj, ak;
 908     ivec    di;
 909     int     svi, sji, ski;
 910
 911     av = ia[1];
 912     ai = ia[2];
 913     aj = ia[3];
 914     ak = ia[4];
 915
 916     sji = pbc_rvec_sub(pbc, x[aj], x[ai], xij);
 917     ski = pbc_rvec_sub(pbc, x[ak], x[ai], xik);
 918     /* 6 Flops */
 919
 920     copy_rvec(f[av], fv);
 921
 922     cfx = c*fv[XX];
 923     cfy = c*fv[YY];
 924     cfz = c*fv[ZZ];
 925     /* 3 Flops */
 926
 927     fj[XX] = a*fv[XX]     -  xik[ZZ]*cfy +  xik[YY]*cfz;
 928     fj[YY] =  xik[ZZ]*cfx + a*fv[YY]     -  xik[XX]*cfz;
 929     fj[ZZ] = -xik[YY]*cfx +  xik[XX]*cfy + a*fv[ZZ];
 930
 931     fk[XX] = b*fv[XX]     +  xij[ZZ]*cfy -  xij[YY]*cfz;
 932     fk[YY] = -xij[ZZ]*cfx + b*fv[YY]     +  xij[XX]*cfz;
 933     fk[ZZ] =  xij[YY]*cfx -  xij[XX]*cfy + b*fv[ZZ];
 934     /* 30 Flops */
 935
 936     f[ai][XX] += fv[XX] - fj[XX] - fk[XX];
 937     f[ai][YY] += fv[YY] - fj[YY] - fk[YY];
 938     f[ai][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
 939     rvec_inc(f[aj], fj);
 940     rvec_inc(f[ak], fk);
 941     /* 15 Flops */
 942
 943     if (g)
 944     {
 945         ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di);
 946         svi = IVEC2IS(di);
 947         ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di);
 948         sji = IVEC2IS(di);
 949         ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, ai), di);
 950         ski = IVEC2IS(di);
 951     }
 952     else if (pbc)
 953     {
 954         svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi);
 955     }
 956     else
 957     {
 958         svi = CENTRAL;
 959     }
 960
 961     if (fshift && (svi != CENTRAL || sji != CENTRAL || ski != CENTRAL))
 962     {
 963         rvec_dec(fshift[svi], fv);
 964         fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX];
 965         fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY];
 966         fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
 967         rvec_inc(fshift[sji], fj);
 968         rvec_inc(fshift[ski], fk);
 969     }
 970
 971     if (VirCorr)
 972     {
 973         rvec xiv;
 974         int  i, j;
 975
 976         pbc_rvec_sub(pbc, x[av], x[ai], xiv);
 977
 978         for (i = 0; i < DIM; i++)
 979         {
 980             for (j = 0; j < DIM; j++)
 981             {
 982                 dxdf[i][j] += -xiv[i]*fv[j] + xij[i]*fj[j] + xik[i]*fk[j];
 983             }
 984         }
 985     }
 986
 987     /* TOTAL: 54 flops */
 988 }
 989
 990 static void spread_vsite4FD(t_iatom ia[], real a, real b, real c,
 991                             rvec x[], rvec f[], rvec fshift[],
 992                             gmx_bool VirCorr, matrix dxdf,
 993                             t_pbc *pbc, t_graph *g)
 994 {
 995     real    d, invl, fproj, a1;
 996     rvec    xvi, xij, xjk, xjl, xix, fv, temp;
 997     atom_id av, ai, aj, ak, al;
 998     ivec    di;
 999     int     svi, sji, skj, slj, m;
1000
1001     av = ia[1];
1002     ai = ia[2];
1003     aj = ia[3];
1004     ak = ia[4];
1005     al = ia[5];
1006
1007     sji = pbc_rvec_sub(pbc, x[aj], x[ai], xij);
1008     skj = pbc_rvec_sub(pbc, x[ak], x[aj], xjk);
1009     slj = pbc_rvec_sub(pbc, x[al], x[aj], xjl);
1010     /* 9 flops */
1011
1012     /* xix goes from i to point x on the plane jkl */
1013     for (m = 0; m < DIM; m++)
1014     {
1015         xix[m] = xij[m] + a*xjk[m] + b*xjl[m];
1016     }
1017     /* 12 flops */
1018
1019     invl = gmx_invsqrt(iprod(xix, xix));
1020     d    = c*invl;
1021     /* 4 + ?10? flops */
1022
1023     copy_rvec(f[av], fv);
1024
1025     fproj = iprod(xix, fv)*invl*invl; /* = (xix . f)/(xix . xix) */
1026
1027     for (m = 0; m < DIM; m++)
1028     {
1029         temp[m] = d*(fv[m] - fproj*xix[m]);
1030     }
1031     /* 16 */
1032
1033     /* c is already calculated in constr_vsite3FD
1034        storing c somewhere will save 35 flops!     */
1035
1036     a1 = 1 - a - b;
1037     for (m = 0; m < DIM; m++)
1038     {
1039         f[ai][m] += fv[m] - temp[m];
1040         f[aj][m] += a1*temp[m];
1041         f[ak][m] += a*temp[m];
1042         f[al][m] += b*temp[m];
1043     }
1044     /* 26 Flops */
1045
1046     if (g)
1047     {
1048         ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di);
1049         svi = IVEC2IS(di);
1050         ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di);
1051         sji = IVEC2IS(di);
1052         ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, aj), di);
1053         skj = IVEC2IS(di);
1054         ivec_sub(SHIFT_IVEC(g, al), SHIFT_IVEC(g, aj), di);
1055         slj = IVEC2IS(di);
1056     }
1057     else if (pbc)
1058     {
1059         svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi);
1060     }
1061     else
1062     {
1063         svi = CENTRAL;
1064     }
1065
1066     if (fshift &&
1067         (svi != CENTRAL || sji != CENTRAL || skj != CENTRAL || slj != CENTRAL))
1068     {
1069         rvec_dec(fshift[svi], fv);
1070         for (m = 0; m < DIM; m++)
1071         {
1072             fshift[CENTRAL][m] += fv[m] - (1 + a + b)*temp[m];
1073             fshift[    sji][m] += temp[m];
1074             fshift[    skj][m] += a*temp[m];
1075             fshift[    slj][m] += b*temp[m];
1076         }
1077     }
1078
1079     if (VirCorr)
1080     {
1081         rvec xiv;
1082         int  i, j;
1083
1084         pbc_rvec_sub(pbc, x[av], x[ai], xiv);
1085
1086         for (i = 0; i < DIM; i++)
1087         {
1088             for (j = 0; j < DIM; j++)
1089             {
1090                 dxdf[i][j] += -xiv[i]*fv[j] + xix[i]*temp[j];
1091             }
1092         }
1093     }
1094
1095     /* TOTAL: 77 flops */
1096 }
1097
1098
1099 static void spread_vsite4FDN(t_iatom ia[], real a, real b, real c,
1100                              rvec x[], rvec f[], rvec fshift[],
1101                              gmx_bool VirCorr, matrix dxdf,
1102                              t_pbc *pbc, t_graph *g)
1103 {
1104     rvec xvi, xij, xik, xil, ra, rb, rja, rjb, rab, rm, rt;
1105     rvec fv, fj, fk, fl;
1106     real invrm, denom;
1107     real cfx, cfy, cfz;
1108     ivec di;
1109     int  av, ai, aj, ak, al;
1110     int  svi, sij, sik, sil;
1111
1112     /* DEBUG: check atom indices */
1113     av = ia[1];
1114     ai = ia[2];
1115     aj = ia[3];
1116     ak = ia[4];
1117     al = ia[5];
1118
1119     copy_rvec(f[av], fv);
1120
1121     sij = pbc_rvec_sub(pbc, x[aj], x[ai], xij);
1122     sik = pbc_rvec_sub(pbc, x[ak], x[ai], xik);
1123     sil = pbc_rvec_sub(pbc, x[al], x[ai], xil);
1124     /* 9 flops */
1125
1126     ra[XX] = a*xik[XX];
1127     ra[YY] = a*xik[YY];
1128     ra[ZZ] = a*xik[ZZ];
1129
1130     rb[XX] = b*xil[XX];
1131     rb[YY] = b*xil[YY];
1132     rb[ZZ] = b*xil[ZZ];
1133
1134     /* 6 flops */
1135
1136     rvec_sub(ra, xij, rja);
1137     rvec_sub(rb, xij, rjb);
1138     rvec_sub(rb, ra, rab);
1139     /* 9 flops */
1140
1141     cprod(rja, rjb, rm);
1142     /* 9 flops */
1143
1144     invrm = gmx_invsqrt(norm2(rm));
1145     denom = invrm*invrm;
1146     /* 5+5+2 flops */
1147
1148     cfx = c*invrm*fv[XX];
1149     cfy = c*invrm*fv[YY];
1150     cfz = c*invrm*fv[ZZ];
1151     /* 6 Flops */
1152
1153     cprod(rm, rab, rt);
1154     /* 9 flops */
1155
1156     rt[XX] *= denom;
1157     rt[YY] *= denom;
1158     rt[ZZ] *= denom;
1159     /* 3flops */
1160
1161     fj[XX] = (        -rm[XX]*rt[XX]) * cfx + ( rab[ZZ]-rm[YY]*rt[XX]) * cfy + (-rab[YY]-rm[ZZ]*rt[XX]) * cfz;
1162     fj[YY] = (-rab[ZZ]-rm[XX]*rt[YY]) * cfx + (        -rm[YY]*rt[YY]) * cfy + ( rab[XX]-rm[ZZ]*rt[YY]) * cfz;
1163     fj[ZZ] = ( rab[YY]-rm[XX]*rt[ZZ]) * cfx + (-rab[XX]-rm[YY]*rt[ZZ]) * cfy + (        -rm[ZZ]*rt[ZZ]) * cfz;
1164     /* 30 flops */
1165
1166     cprod(rjb, rm, rt);
1167     /* 9 flops */
1168
1169     rt[XX] *= denom*a;
1170     rt[YY] *= denom*a;
1171     rt[ZZ] *= denom*a;
1172     /* 3flops */
1173
1174     fk[XX] = (          -rm[XX]*rt[XX]) * cfx + (-a*rjb[ZZ]-rm[YY]*rt[XX]) * cfy + ( a*rjb[YY]-rm[ZZ]*rt[XX]) * cfz;
1175     fk[YY] = ( a*rjb[ZZ]-rm[XX]*rt[YY]) * cfx + (          -rm[YY]*rt[YY]) * cfy + (-a*rjb[XX]-rm[ZZ]*rt[YY]) * cfz;
1176     fk[ZZ] = (-a*rjb[YY]-rm[XX]*rt[ZZ]) * cfx + ( a*rjb[XX]-rm[YY]*rt[ZZ]) * cfy + (          -rm[ZZ]*rt[ZZ]) * cfz;
1177     /* 36 flops */
1178
1179     cprod(rm, rja, rt);
1180     /* 9 flops */
1181
1182     rt[XX] *= denom*b;
1183     rt[YY] *= denom*b;
1184     rt[ZZ] *= denom*b;
1185     /* 3flops */
1186
1187     fl[XX] = (          -rm[XX]*rt[XX]) * cfx + ( b*rja[ZZ]-rm[YY]*rt[XX]) * cfy + (-b*rja[YY]-rm[ZZ]*rt[XX]) * cfz;
1188     fl[YY] = (-b*rja[ZZ]-rm[XX]*rt[YY]) * cfx + (          -rm[YY]*rt[YY]) * cfy + ( b*rja[XX]-rm[ZZ]*rt[YY]) * cfz;
1189     fl[ZZ] = ( b*rja[YY]-rm[XX]*rt[ZZ]) * cfx + (-b*rja[XX]-rm[YY]*rt[ZZ]) * cfy + (          -rm[ZZ]*rt[ZZ]) * cfz;
1190     /* 36 flops */
1191
1192     f[ai][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
1193     f[ai][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
1194     f[ai][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
1195     rvec_inc(f[aj], fj);
1196     rvec_inc(f[ak], fk);
1197     rvec_inc(f[al], fl);
1198     /* 21 flops */
1199
1200     if (g)
1201     {
1202         ivec_sub(SHIFT_IVEC(g, av), SHIFT_IVEC(g, ai), di);
1203         svi = IVEC2IS(di);
1204         ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di);
1205         sij = IVEC2IS(di);
1206         ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, ai), di);
1207         sik = IVEC2IS(di);
1208         ivec_sub(SHIFT_IVEC(g, al), SHIFT_IVEC(g, ai), di);
1209         sil = IVEC2IS(di);
1210     }
1211     else if (pbc)
1212     {
1213         svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi);
1214     }
1215     else
1216     {
1217         svi = CENTRAL;
1218     }
1219
1220     if (fshift && (svi != CENTRAL || sij != CENTRAL || sik != CENTRAL || sil != CENTRAL))
1221     {
1222         rvec_dec(fshift[svi], fv);
1223         fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
1224         fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
1225         fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
1226         rvec_inc(fshift[sij], fj);
1227         rvec_inc(fshift[sik], fk);
1228         rvec_inc(fshift[sil], fl);
1229     }
1230
1231     if (VirCorr)
1232     {
1233         rvec xiv;
1234         int  i, j;
1235
1236         pbc_rvec_sub(pbc, x[av], x[ai], xiv);
1237
1238         for (i = 0; i < DIM; i++)
1239         {
1240             for (j = 0; j < DIM; j++)
1241             {
1242                 dxdf[i][j] += -xiv[i]*fv[j] + xij[i]*fj[j] + xik[i]*fk[j] + xil[i]*fl[j];
1243             }
1244         }
1245     }
1246
1247     /* Total: 207 flops (Yuck!) */
1248 }
1249
1250
1251 static int spread_vsiten(t_iatom ia[], t_iparams ip[],
1252                          rvec x[], rvec f[], rvec fshift[],
1253                          t_pbc *pbc, t_graph *g)
1254 {
1255     rvec xv, dx, fi;
1256     int  n3, av, i, ai;
1257     real a;
1258     ivec di;
1259     int  siv;
1260
1261     n3 = 3*ip[ia[0]].vsiten.n;
1262     av = ia[1];
1263     copy_rvec(x[av], xv);
1264
1265     for (i = 0; i < n3; i += 3)
1266     {
1267         ai = ia[i+2];
1268         if (g)
1269         {
1270             ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, av), di);
1271             siv = IVEC2IS(di);
1272         }
1273         else if (pbc)
1274         {
1275             siv = pbc_dx_aiuc(pbc, x[ai], xv, dx);
1276         }
1277         else
1278         {
1279             siv = CENTRAL;
1280         }
1281         a = ip[ia[i]].vsiten.a;
1282         svmul(a, f[av], fi);
1283         rvec_inc(f[ai], fi);
1284         if (fshift && siv != CENTRAL)
1285         {
1286             rvec_inc(fshift[siv], fi);
1287             rvec_dec(fshift[CENTRAL], fi);
1288         }
1289         /* 6 Flops */
1290     }
1291
1292     return n3;
1293 }
1294
1295
1296 static int vsite_count(const t_ilist *ilist, int ftype)
1297 {
1298     if (ftype == F_VSITEN)
1299     {
1300         return ilist[ftype].nr/3;
1301     }
1302     else
1303     {
1304         return ilist[ftype].nr/(1 + interaction_function[ftype].nratoms);
1305     }
1306 }
1307
1308 static void spread_vsite_f_thread(gmx_vsite_t *vsite,
1309                                   rvec x[], rvec f[], rvec *fshift,
1310                                   gmx_bool VirCorr, matrix dxdf,
1311                                   t_iparams ip[], t_ilist ilist[],
1312                                   t_graph *g, t_pbc *pbc_null)
1313 {
1314     gmx_bool   bPBCAll;
1315     real       a1, b1, c1;
1316     int        i, inc, m, nra, nr, tp, ftype;
1317     t_iatom   *ia;
1318     t_pbc     *pbc_null2;
1319     int       *vsite_pbc;
1320
1321     if (VirCorr)
1322     {
1323         clear_mat(dxdf);
1324     }
1325
1326     bPBCAll = (pbc_null != NULL && !vsite->bHaveChargeGroups);
1327
1328     /* this loop goes backwards to be able to build *
1329      * higher type vsites from lower types         */
1330     pbc_null2 = NULL;
1331     vsite_pbc = NULL;
1332     for (ftype = F_NRE-1; (ftype >= 0); ftype--)
1333     {
1334         if ((interaction_function[ftype].flags & IF_VSITE) &&
1335             ilist[ftype].nr > 0)
1336         {
1337             nra    = interaction_function[ftype].nratoms;
1338             inc    = 1 + nra;
1339             nr     = ilist[ftype].nr;
1340             ia     = ilist[ftype].iatoms;
1341
1342             if (bPBCAll)
1343             {
1344                 pbc_null2 = pbc_null;
1345             }
1346             else if (pbc_null != NULL)
1347             {
1348                 vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2];
1349             }
1350
1351             for (i = 0; i < nr; )
1352             {
1353                 if (vsite_pbc != NULL)
1354                 {
1355                     if (vsite_pbc[i/(1+nra)] > -2)
1356                     {
1357                         pbc_null2 = pbc_null;
1358                     }
1359                     else
1360                     {
1361                         pbc_null2 = NULL;
1362                     }
1363                 }
1364
1365                 tp   = ia[0];
1366
1367                 /* Constants for constructing */
1368                 a1   = ip[tp].vsite.a;
1369                 /* Construct the vsite depending on type */
1370                 switch (ftype)
1371                 {
1372                     case F_VSITE2:
1373                         spread_vsite2(ia, a1, x, f, fshift, pbc_null2, g);
1374                         break;
1375                     case F_VSITE3:
1376                         b1 = ip[tp].vsite.b;
1377                         spread_vsite3(ia, a1, b1, x, f, fshift, pbc_null2, g);
1378                         break;
1379                     case F_VSITE3FD:
1380                         b1 = ip[tp].vsite.b;
1381                         spread_vsite3FD(ia, a1, b1, x, f, fshift, VirCorr, dxdf, pbc_null2, g);
1382                         break;
1383                     case F_VSITE3FAD:
1384                         b1 = ip[tp].vsite.b;
1385                         spread_vsite3FAD(ia, a1, b1, x, f, fshift, VirCorr, dxdf, pbc_null2, g);
1386                         break;
1387                     case F_VSITE3OUT:
1388                         b1 = ip[tp].vsite.b;
1389                         c1 = ip[tp].vsite.c;
1390                         spread_vsite3OUT(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2, g);
1391                         break;
1392                     case F_VSITE4FD:
1393                         b1 = ip[tp].vsite.b;
1394                         c1 = ip[tp].vsite.c;
1395                         spread_vsite4FD(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2, g);
1396                         break;
1397                     case F_VSITE4FDN:
1398                         b1 = ip[tp].vsite.b;
1399                         c1 = ip[tp].vsite.c;
1400                         spread_vsite4FDN(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2, g);
1401                         break;
1402                     case F_VSITEN:
1403                         inc = spread_vsiten(ia, ip, x, f, fshift, pbc_null2, g);
1404                         break;
1405                     default:
1406                         gmx_fatal(FARGS, "No such vsite type %d in %s, line %d",
1407                                   ftype, __FILE__, __LINE__);
1408                 }
1409                 clear_rvec(f[ia[1]]);
1410
1411                 /* Increment loop variables */
1412                 i  += inc;
1413                 ia += inc;
1414             }
1415         }
1416     }
1417 }
1418
1419 void spread_vsite_f(gmx_vsite_t *vsite,
1420                     rvec x[], rvec f[], rvec *fshift,
1421                     gmx_bool VirCorr, matrix vir,
1422                     t_nrnb *nrnb, t_idef *idef,
1423                     int ePBC, gmx_bool bMolPBC, t_graph *g, matrix box,
1424                     t_commrec *cr)
1425 {
1426     t_pbc pbc, *pbc_null;
1427     int   th;
1428
1429     /* We only need to do pbc when we have inter-cg vsites */
1430     if ((DOMAINDECOMP(cr) || bMolPBC) && vsite->n_intercg_vsite)
1431     {
1432         /* This is wasting some CPU time as we now do this multiple times
1433          * per MD step. But how often do we have vsites with full pbc?
1434          */
1435         pbc_null = set_pbc_dd(&pbc, ePBC, cr->dd, FALSE, box);
1436     }
1437     else
1438     {
1439         pbc_null = NULL;
1440     }
1441
1442     if (DOMAINDECOMP(cr))
1443     {
1444         dd_clear_f_vsites(cr->dd, f);
1445     }
1446
1447     if (vsite->nthreads == 1)
1448     {
1449         spread_vsite_f_thread(vsite,
1450                               x, f, fshift,
1451                               VirCorr, vsite->tdata[0].dxdf,
1452                               idef->iparams, idef->il,
1453                               g, pbc_null);
1454     }
1455     else
1456     {
1457         /* First spread the vsites that might depend on other vsites */
1458         spread_vsite_f_thread(vsite,
1459                               x, f, fshift,
1460                               VirCorr, vsite->tdata[vsite->nthreads].dxdf,
1461                               idef->iparams,
1462                               vsite->tdata[vsite->nthreads].ilist,
1463                               g, pbc_null);
1464
1465 #pragma omp parallel num_threads(vsite->nthreads)
1466         {
1467             int   thread;
1468             rvec *fshift_t;
1469
1470             thread = gmx_omp_get_thread_num();
1471
1472             if (thread == 0 || fshift == NULL)
1473             {
1474                 fshift_t = fshift;
1475             }
1476             else
1477             {
1478                 int i;
1479
1480                 fshift_t = vsite->tdata[thread].fshift;
1481
1482                 for (i = 0; i < SHIFTS; i++)
1483                 {
1484                     clear_rvec(fshift_t[i]);
1485                 }
1486             }
1487
1488             spread_vsite_f_thread(vsite,
1489                                   x, f, fshift_t,
1490                                   VirCorr, vsite->tdata[thread].dxdf,
1491                                   idef->iparams,
1492                                   vsite->tdata[thread].ilist,
1493                                   g, pbc_null);
1494         }
1495
1496         if (fshift != NULL)
1497         {
1498             int i;
1499
1500             for (th = 1; th < vsite->nthreads; th++)
1501             {
1502                 for (i = 0; i < SHIFTS; i++)
1503                 {
1504                     rvec_inc(fshift[i], vsite->tdata[th].fshift[i]);
1505                 }
1506             }
1507         }
1508     }
1509
1510     if (VirCorr)
1511     {
1512         int i, j;
1513
1514         for (th = 0; th < (vsite->nthreads == 1 ? 1 : vsite->nthreads+1); th++)
1515         {
1516             for (i = 0; i < DIM; i++)
1517             {
1518                 for (j = 0; j < DIM; j++)
1519                 {
1520                     vir[i][j] += -0.5*vsite->tdata[th].dxdf[i][j];
1521                 }
1522             }
1523         }
1524     }
1525
1526     if (DOMAINDECOMP(cr))
1527     {
1528         dd_move_f_vsites(cr->dd, f, fshift);
1529     }
1530
1531     inc_nrnb(nrnb, eNR_VSITE2,   vsite_count(idef->il, F_VSITE2));
1532     inc_nrnb(nrnb, eNR_VSITE3,   vsite_count(idef->il, F_VSITE3));
1533     inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef->il, F_VSITE3FD));
1534     inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef->il, F_VSITE3FAD));
1535     inc_nrnb(nrnb, eNR_VSITE3OUT, vsite_count(idef->il, F_VSITE3OUT));
1536     inc_nrnb(nrnb, eNR_VSITE4FD, vsite_count(idef->il, F_VSITE4FD));
1537     inc_nrnb(nrnb, eNR_VSITE4FDN, vsite_count(idef->il, F_VSITE4FDN));
1538     inc_nrnb(nrnb, eNR_VSITEN,   vsite_count(idef->il, F_VSITEN));
1539 }
1540
1541 static int *atom2cg(t_block *cgs)
1542 {
1543     int *a2cg, cg, i;
1544
1545     snew(a2cg, cgs->index[cgs->nr]);
1546     for (cg = 0; cg < cgs->nr; cg++)
1547     {
1548         for (i = cgs->index[cg]; i < cgs->index[cg+1]; i++)
1549         {
1550             a2cg[i] = cg;
1551         }
1552     }
1553
1554     return a2cg;
1555 }
1556
1557 static int count_intercg_vsite(gmx_mtop_t *mtop,
1558                                gmx_bool   *bHaveChargeGroups)
1559 {
1560     int             mb, mt, ftype, nral, i, cg, a;
1561     gmx_molblock_t *molb;
1562     gmx_moltype_t  *molt;
1563     int            *a2cg;
1564     t_ilist        *il;
1565     t_iatom        *ia;
1566     int             n_intercg_vsite;
1567
1568     *bHaveChargeGroups = FALSE;
1569
1570     n_intercg_vsite = 0;
1571     for (mb = 0; mb < mtop->nmolblock; mb++)
1572     {
1573         molb = &mtop->molblock[mb];
1574         molt = &mtop->moltype[molb->type];
1575
1576         if (molt->cgs.nr < molt->atoms.nr)
1577         {
1578             *bHaveChargeGroups = TRUE;
1579         }
1580
1581         a2cg = atom2cg(&molt->cgs);
1582         for (ftype = 0; ftype < F_NRE; ftype++)
1583         {
1584             if (interaction_function[ftype].flags & IF_VSITE)
1585             {
1586                 nral = NRAL(ftype);
1587                 il   = &molt->ilist[ftype];
1588                 ia   = il->iatoms;
1589                 for (i = 0; i < il->nr; i += 1+nral)
1590                 {
1591                     cg = a2cg[ia[1+i]];
1592                     for (a = 1; a < nral; a++)
1593                     {
1594                         if (a2cg[ia[1+a]] != cg)
1595                         {
1596                             n_intercg_vsite += molb->nmol;
1597                             break;
1598                         }
1599                     }
1600                 }
1601             }
1602         }
1603         sfree(a2cg);
1604     }
1605
1606     return n_intercg_vsite;
1607 }
1608
1609 static int **get_vsite_pbc(t_iparams *iparams, t_ilist *ilist,
1610                            t_atom *atom, t_mdatoms *md,
1611                            t_block *cgs, int *a2cg)
1612 {
1613     int      ftype, nral, i, j, vsi, vsite, cg_v, cg_c, a, nc3 = 0;
1614     t_ilist *il;
1615     t_iatom *ia;
1616     int    **vsite_pbc, *vsite_pbc_f;
1617     char    *pbc_set;
1618     gmx_bool bViteOnlyCG_and_FirstAtom;
1619
1620     /* Make an array that tells if the pbc of an atom is set */
1621     snew(pbc_set, cgs->index[cgs->nr]);
1622     /* PBC is set for all non vsites */
1623     for (a = 0; a < cgs->index[cgs->nr]; a++)
1624     {
1625         if ((atom && atom[a].ptype != eptVSite) ||
1626             (md   && md->ptype[a]  != eptVSite))
1627         {
1628             pbc_set[a] = 1;
1629         }
1630     }
1631
1632     snew(vsite_pbc, F_VSITEN-F_VSITE2+1);
1633
1634     for (ftype = 0; ftype < F_NRE; ftype++)
1635     {
1636         if (interaction_function[ftype].flags & IF_VSITE)
1637         {
1638             nral = NRAL(ftype);
1639             il   = &ilist[ftype];
1640             ia   = il->iatoms;
1641
1642             snew(vsite_pbc[ftype-F_VSITE2], il->nr/(1+nral));
1643             vsite_pbc_f = vsite_pbc[ftype-F_VSITE2];
1644
1645             i = 0;
1646             while (i < il->nr)
1647             {
1648                 vsi   = i/(1+nral);
1649                 vsite = ia[i+1];
1650                 cg_v  = a2cg[vsite];
1651                 /* A value of -2 signals that this vsite and its contructing
1652                  * atoms are all within the same cg, so no pbc is required.
1653                  */
1654                 vsite_pbc_f[vsi] = -2;
1655                 /* Check if constructing atoms are outside the vsite's cg */
1656                 nc3 = 0;
1657                 if (ftype == F_VSITEN)
1658                 {
1659                     nc3 = 3*iparams[ia[i]].vsiten.n;
1660                     for (j = 0; j < nc3; j += 3)
1661                     {
1662                         if (a2cg[ia[i+j+2]] != cg_v)
1663                         {
1664                             vsite_pbc_f[vsi] = -1;
1665                         }
1666                     }
1667                 }
1668                 else
1669                 {
1670                     for (a = 1; a < nral; a++)
1671                     {
1672                         if (a2cg[ia[i+1+a]] != cg_v)
1673                         {
1674                             vsite_pbc_f[vsi] = -1;
1675                         }
1676                     }
1677                 }
1678                 if (vsite_pbc_f[vsi] == -1)
1679                 {
1680                     /* Check if this is the first processed atom of a vsite only cg */
1681                     bViteOnlyCG_and_FirstAtom = TRUE;
1682                     for (a = cgs->index[cg_v]; a < cgs->index[cg_v+1]; a++)
1683                     {
1684                         /* Non-vsites already have pbc set, so simply check for pbc_set */
1685                         if (pbc_set[a])
1686                         {
1687                             bViteOnlyCG_and_FirstAtom = FALSE;
1688                             break;
1689                         }
1690                     }
1691                     if (bViteOnlyCG_and_FirstAtom)
1692                     {
1693                         /* First processed atom of a vsite only charge group.
1694                          * The pbc of the input coordinates to construct_vsites
1695                          * should be preserved.
1696                          */
1697                         vsite_pbc_f[vsi] = vsite;
1698                     }
1699                     else if (cg_v != a2cg[ia[1+i+1]])
1700                     {
1701                         /* This vsite has a different charge group index
1702                          * than it's first constructing atom
1703                          * and the charge group has more than one atom,
1704                          * search for the first normal particle
1705                          * or vsite that already had its pbc defined.
1706                          * If nothing is found, use full pbc for this vsite.
1707                          */
1708                         for (a = cgs->index[cg_v]; a < cgs->index[cg_v+1]; a++)
1709                         {
1710                             if (a != vsite && pbc_set[a])
1711                             {
1712                                 vsite_pbc_f[vsi] = a;
1713                                 if (gmx_debug_at)
1714                                 {
1715                                     fprintf(debug, "vsite %d match pbc with atom %d\n",
1716                                             vsite+1, a+1);
1717                                 }
1718                                 break;
1719                             }
1720                         }
1721                         if (gmx_debug_at)
1722                         {
1723                             fprintf(debug, "vsite atom %d  cg %d - %d pbc atom %d\n",
1724                                     vsite+1, cgs->index[cg_v]+1, cgs->index[cg_v+1],
1725                                     vsite_pbc_f[vsi]+1);
1726                         }
1727                     }
1728                 }
1729                 if (ftype == F_VSITEN)
1730                 {
1731                     /* The other entries in vsite_pbc_f are not used for center vsites */
1732                     i += nc3;
1733                 }
1734                 else
1735                 {
1736                     i += 1+nral;
1737                 }
1738
1739                 /* This vsite now has its pbc defined */
1740                 pbc_set[vsite] = 1;
1741             }
1742         }
1743     }
1744
1745     sfree(pbc_set);
1746
1747     return vsite_pbc;
1748 }
1749
1750
1751 gmx_vsite_t *init_vsite(gmx_mtop_t *mtop, t_commrec *cr,
1752                         gmx_bool bSerial_NoPBC)
1753 {
1754     int            nvsite, i;
1755     int           *a2cg, cg;
1756     gmx_vsite_t   *vsite;
1757     int            mt;
1758     gmx_moltype_t *molt;
1759     int            nthreads;
1760
1761     /* check if there are vsites */
1762     nvsite = 0;
1763     for (i = 0; i < F_NRE; i++)
1764     {
1765         if (interaction_function[i].flags & IF_VSITE)
1766         {
1767             nvsite += gmx_mtop_ftype_count(mtop, i);
1768         }
1769     }
1770
1771     if (nvsite == 0)
1772     {
1773         return NULL;
1774     }
1775
1776     snew(vsite, 1);
1777
1778     vsite->n_intercg_vsite = count_intercg_vsite(mtop,
1779                                                  &vsite->bHaveChargeGroups);
1780
1781     /* If we don't have charge groups, the vsite follows its own pbc */
1782     if (!bSerial_NoPBC &&
1783         vsite->bHaveChargeGroups &&
1784         vsite->n_intercg_vsite > 0 && DOMAINDECOMP(cr))
1785     {
1786         vsite->nvsite_pbc_molt = mtop->nmoltype;
1787         snew(vsite->vsite_pbc_molt, vsite->nvsite_pbc_molt);
1788         for (mt = 0; mt < mtop->nmoltype; mt++)
1789         {
1790             molt = &mtop->moltype[mt];
1791             /* Make an atom to charge group index */
1792             a2cg = atom2cg(&molt->cgs);
1793             vsite->vsite_pbc_molt[mt] = get_vsite_pbc(mtop->ffparams.iparams,
1794                                                       molt->ilist,
1795                                                       molt->atoms.atom, NULL,
1796                                                       &molt->cgs, a2cg);
1797             sfree(a2cg);
1798         }
1799
1800         snew(vsite->vsite_pbc_loc_nalloc, F_VSITEN-F_VSITE2+1);
1801         snew(vsite->vsite_pbc_loc, F_VSITEN-F_VSITE2+1);
1802     }
1803
1804     if (bSerial_NoPBC)
1805     {
1806         vsite->nthreads = 1;
1807     }
1808     else
1809     {
1810         vsite->nthreads = gmx_omp_nthreads_get(emntVSITE);
1811     }
1812     if (!bSerial_NoPBC)
1813     {
1814         /* We need one extra thread data structure for the overlap vsites */
1815         snew(vsite->tdata, vsite->nthreads+1);
1816     }
1817
1818     vsite->th_ind        = NULL;
1819     vsite->th_ind_nalloc = 0;
1820
1821     return vsite;
1822 }
1823
1824 static void prepare_vsite_thread(const t_ilist      *ilist,
1825                                  gmx_vsite_thread_t *vsite_th)
1826 {
1827     int ftype;
1828
1829     for (ftype = 0; ftype < F_NRE; ftype++)
1830     {
1831         if (interaction_function[ftype].flags & IF_VSITE)
1832         {
1833             if (ilist[ftype].nr > vsite_th->ilist[ftype].nalloc)
1834             {
1835                 vsite_th->ilist[ftype].nalloc = over_alloc_large(ilist[ftype].nr);
1836                 srenew(vsite_th->ilist[ftype].iatoms, vsite_th->ilist[ftype].nalloc);
1837             }
1838
1839             vsite_th->ilist[ftype].nr = 0;
1840         }
1841     }
1842 }
1843
1844 void split_vsites_over_threads(const t_ilist   *ilist,
1845                                const t_mdatoms *mdatoms,
1846                                gmx_bool         bLimitRange,
1847                                gmx_vsite_t     *vsite)
1848 {
1849     int      th;
1850     int      vsite_atom_range, natperthread;
1851     int     *th_ind;
1852     int      ftype;
1853     t_iatom *iat;
1854     t_ilist *il_th;
1855     int      nral1, inc, i, j;
1856
1857     if (vsite->nthreads == 1)
1858     {
1859         /* Nothing to do */
1860         return;
1861     }
1862
1863 #pragma omp parallel for num_threads(vsite->nthreads) schedule(static)
1864     for (th = 0; th < vsite->nthreads; th++)
1865     {
1866         prepare_vsite_thread(ilist, &vsite->tdata[th]);
1867     }
1868     /* Master threads does the (potential) overlap vsites */
1869     prepare_vsite_thread(ilist, &vsite->tdata[vsite->nthreads]);
1870
1871     /* The current way of distributing the vsites over threads in primitive.
1872      * We divide the atom range 0 - natoms_in_vsite uniformly over threads,
1873      * without taking into account how the vsites are distributed.
1874      * Without domain decomposition we bLimitRange=TRUE and we at least
1875      * tighten the upper bound of the range (useful for common systems
1876      * such as a vsite-protein in 3-site water).
1877      */
1878     if (bLimitRange)
1879     {
1880         vsite_atom_range = -1;
1881         for (ftype = 0; ftype < F_NRE; ftype++)
1882         {
1883             if ((interaction_function[ftype].flags & IF_VSITE) &&
1884                 ftype != F_VSITEN)
1885             {
1886                 nral1 = 1 + NRAL(ftype);
1887                 iat   = ilist[ftype].iatoms;
1888                 for (i = 0; i < ilist[ftype].nr; i += nral1)
1889                 {
1890                     for (j = i+1; j < i+nral1; j++)
1891                     {
1892                         vsite_atom_range = max(vsite_atom_range, iat[j]);
1893                     }
1894                 }
1895             }
1896         }
1897         vsite_atom_range++;
1898     }
1899     else
1900     {
1901         vsite_atom_range = mdatoms->homenr;
1902     }
1903     natperthread = (vsite_atom_range + vsite->nthreads - 1)/vsite->nthreads;
1904
1905     if (debug)
1906     {
1907         fprintf(debug, "virtual site thread dist: natoms %d, range %d, natperthread %d\n", mdatoms->nr, vsite_atom_range, natperthread);
1908     }
1909
1910     /* To simplify the vsite assignment, we make an index which tells us
1911      * to which thread particles, both non-vsites and vsites, are assigned.
1912      */
1913     if (mdatoms->nr > vsite->th_ind_nalloc)
1914     {
1915         vsite->th_ind_nalloc = over_alloc_large(mdatoms->nr);
1916         srenew(vsite->th_ind, vsite->th_ind_nalloc);
1917     }
1918     th_ind = vsite->th_ind;
1919     th     = 0;
1920     for (i = 0; i < mdatoms->nr; i++)
1921     {
1922         if (mdatoms->ptype[i] == eptVSite)
1923         {
1924             /* vsites are not assigned to a thread yet */
1925             th_ind[i] = -1;
1926         }
1927         else
1928         {
1929             /* assign non-vsite particles to thread th */
1930             th_ind[i] = th;
1931         }
1932         if (i == (th + 1)*natperthread && th < vsite->nthreads)
1933         {
1934             th++;
1935         }
1936     }
1937
1938     for (ftype = 0; ftype < F_NRE; ftype++)
1939     {
1940         if ((interaction_function[ftype].flags & IF_VSITE) &&
1941             ftype != F_VSITEN)
1942         {
1943             nral1 = 1 + NRAL(ftype);
1944             inc   = nral1;
1945             iat   = ilist[ftype].iatoms;
1946             for (i = 0; i < ilist[ftype].nr; )
1947             {
1948                 th = iat[1+i]/natperthread;
1949                 /* We would like to assign this vsite the thread th,
1950                  * but it might depend on atoms outside the atom range of th
1951                  * or on another vsite not assigned to thread th.
1952                  */
1953                 if (ftype != F_VSITEN)
1954                 {
1955                     for (j = i+2; j < i+nral1; j++)
1956                     {
1957                         if (th_ind[iat[j]] != th)
1958                         {
1959                             /* Some constructing atoms are not assigned to
1960                              * thread th, move this vsite to a separate batch.
1961                              */
1962                             th = vsite->nthreads;
1963                         }
1964                     }
1965                 }
1966                 else
1967                 {
1968                     inc = iat[i];
1969                     for (j = i+2; j < i+inc; j += 3)
1970                     {
1971                         if (th_ind[iat[j]] != th)
1972                         {
1973                             th = vsite->nthreads;
1974                         }
1975                     }
1976                 }
1977                 /* Copy this vsite to the thread data struct of thread th */
1978                 il_th = &vsite->tdata[th].ilist[ftype];
1979                 for (j = i; j < i+inc; j++)
1980                 {
1981                     il_th->iatoms[il_th->nr++] = iat[j];
1982                 }
1983                 /* Update this vsite's thread index entry */
1984                 th_ind[iat[1+i]] = th;
1985
1986                 i += inc;
1987             }
1988         }
1989     }
1990
1991     if (debug)
1992     {
1993         for (ftype = 0; ftype < F_NRE; ftype++)
1994         {
1995             if ((interaction_function[ftype].flags & IF_VSITE) &&
1996                 ilist[ftype].nr > 0)
1997             {
1998                 fprintf(debug, "%-20s thread dist:",
1999                         interaction_function[ftype].longname);
2000                 for (th = 0; th < vsite->nthreads+1; th++)
2001                 {
2002                     fprintf(debug, " %4d", vsite->tdata[th].ilist[ftype].nr);
2003                 }
2004                 fprintf(debug, "\n");
2005             }
2006         }
2007     }
2008 }
2009
2010 void set_vsite_top(gmx_vsite_t *vsite, gmx_localtop_t *top, t_mdatoms *md,
2011                    t_commrec *cr)
2012 {
2013     int *a2cg;
2014
2015     if (vsite->n_intercg_vsite > 0)
2016     {
2017         if (vsite->bHaveChargeGroups)
2018         {
2019             /* Make an atom to charge group index */
2020             a2cg                 = atom2cg(&top->cgs);
2021             vsite->vsite_pbc_loc = get_vsite_pbc(top->idef.iparams,
2022                                                  top->idef.il, NULL, md,
2023                                                  &top->cgs, a2cg);
2024             sfree(a2cg);
2025         }
2026
2027     }
2028
2029     if (vsite->nthreads > 1)
2030     {
2031         if (vsite->bHaveChargeGroups)
2032         {
2033             gmx_fatal(FARGS, "The combination of threading, virtual sites and charge groups is not implemented");
2034         }
2035
2036         split_vsites_over_threads(top->idef.il, md, !DOMAINDECOMP(cr), vsite);
2037     }
2038 }