src/gromacs/mdlib/vsite.c

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