src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel314_c_adress.c

   1 /*
   2  * Copyright (c) Erik Lindahl, David van der Spoel 2003
   3  *
   4  * This file is generated automatically at compile time
   5  * by the program mknb in the Gromacs distribution.
   6  *
   7  * Options used when generation this file:
   8  * Language:         c
   9  * Precision:        single
  10  * Threads:          yes
  11  * Software invsqrt: no
  12  * PowerPC invsqrt:  no
  13  * Prefetch forces:  no
  14  * Adress kernel:  yes
  15  * Comments:         no
  16  */
  17 #ifdef HAVE_CONFIG_H
  18 #include<config.h>
  19 #endif
  20 #ifdef GMX_THREAD_SHM_FDECOMP
  21 #include<thread_mpi.h>
  22 #endif
  23 #define ALMOST_ZERO 1e-30
  24 #define ALMOST_ONE 1-(1e-30)
  25 #include<math.h>
  26
  27 #include "nb_kernel314_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel314_adress_cg
  33  * Coulomb interaction:     Tabulated
  34  * VdW interaction:         Lennard-Jones
  35  * water optimization:      pairs of TIP4P interactions
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel314_adress_cg(
  39                     int *           p_nri,
  40                     int *           iinr,
  41                     int *           jindex,
  42                     int *           jjnr,
  43                     int *           shift,
  44                     real *         shiftvec,
  45                     real *         fshift,
  46                     int *           gid,
  47                     real *         pos,
  48                     real *         faction,
  49                     real *         charge,
  50                     real *         p_facel,
  51                     real *         p_krf,
  52                     real *         p_crf,
  53                     real *         Vc,
  54                     int *           type,
  55                     int *           p_ntype,
  56                     real *         vdwparam,
  57                     real *         Vvdw,
  58                     real *         p_tabscale,
  59                     real *         VFtab,
  60                     real *         invsqrta,
  61                     real *         dvda,
  62                     real *         p_gbtabscale,
  63                     real *         GBtab,
  64                     int *           p_nthreads,
  65                     int *           count,
  66                     void *          mtx,
  67                     int *           outeriter,
  68                     int *           inneriter,
  69                     real           force_cap,
  70                     real *         wf)
  71 {
  72     int           nri,ntype,nthreads;
  73     real         facel,krf,crf,tabscale,gbtabscale;
  74     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
  75     int           nn0,nn1,nouter,ninner;
  76     real         shX,shY,shZ;
  77     real         fscal,tx,ty,tz;
  78     real         rinvsq;
  79     real         qq,vcoul,vctot;
  80     int           tj;
  81     real         rinvsix;
  82     real         Vvdw6,Vvdwtot;
  83     real         Vvdw12;
  84     real         r,rt,eps,eps2;
  85     int           n0,nnn;
  86     real         Y,F,Geps,Heps2,Fp,VV;
  87     real         FF;
  88     real         fijC;
  89     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
  90     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
  91     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
  92     real         ix4,iy4,iz4,fix4,fiy4,fiz4;
  93     real         jx1,jy1,jz1;
  94     real         jx2,jy2,jz2,fjx2,fjy2,fjz2;
  95     real         jx3,jy3,jz3,fjx3,fjy3,fjz3;
  96     real         jx4,jy4,jz4,fjx4,fjy4,fjz4;
  97     real         dx11,dy11,dz11,rsq11;
  98     real         dx22,dy22,dz22,rsq22,rinv22;
  99     real         dx23,dy23,dz23,rsq23,rinv23;
 100     real         dx24,dy24,dz24,rsq24,rinv24;
 101     real         dx32,dy32,dz32,rsq32,rinv32;
 102     real         dx33,dy33,dz33,rsq33,rinv33;
 103     real         dx34,dy34,dz34,rsq34,rinv34;
 104     real         dx42,dy42,dz42,rsq42,rinv42;
 105     real         dx43,dy43,dz43,rsq43,rinv43;
 106     real         dx44,dy44,dz44,rsq44,rinv44;
 107     real         qH,qM,qqMM,qqMH,qqHH;
 108     real         c6,c12;
 109     real         weight_cg1, weight_cg2, weight_product;
 110     real         hybscal;
 111
 112     nri              = *p_nri;
 113     ntype            = *p_ntype;
 114     nthreads         = *p_nthreads;
 115     facel            = *p_facel;
 116     krf              = *p_krf;
 117     crf              = *p_crf;
 118     tabscale         = *p_tabscale;
 119     ii               = iinr[0];
 120     qH               = charge[ii+1];
 121     qM               = charge[ii+3];
 122     qqMM             = facel*qM*qM;
 123     qqMH             = facel*qM*qH;
 124     qqHH             = facel*qH*qH;
 125     tj               = 2*(ntype+1)*type[ii];
 126     c6               = vdwparam[tj];
 127     c12              = vdwparam[tj+1];
 128
 129     nouter           = 0;
 130     ninner           = 0;
 131
 132     do
 133     {
 134         #ifdef GMX_THREAD_SHM_FDECOMP
 135         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 136         nn0              = *count;
 137         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 138         *count           = nn1;
 139         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 140         if(nn1>nri) nn1=nri;
 141         #else
 142         nn0 = 0;
 143         nn1 = nri;
 144         #endif
 145
 146         for(n=nn0; (n<nn1); n++)
 147         {
 148             is3              = 3*shift[n];
 149             shX              = shiftvec[is3];
 150             shY              = shiftvec[is3+1];
 151             shZ              = shiftvec[is3+2];
 152             nj0              = jindex[n];
 153             nj1              = jindex[n+1];
 154             ii               = iinr[n];
 155             ii3              = 3*ii;
 156             ix1              = shX + pos[ii3+0];
 157             iy1              = shY + pos[ii3+1];
 158             iz1              = shZ + pos[ii3+2];
 159             ix2              = shX + pos[ii3+3];
 160             iy2              = shY + pos[ii3+4];
 161             iz2              = shZ + pos[ii3+5];
 162             ix3              = shX + pos[ii3+6];
 163             iy3              = shY + pos[ii3+7];
 164             iz3              = shZ + pos[ii3+8];
 165             ix4              = shX + pos[ii3+9];
 166             iy4              = shY + pos[ii3+10];
 167             iz4              = shZ + pos[ii3+11];
 168             weight_cg1       = wf[ii];
 169             vctot            = 0;
 170             Vvdwtot          = 0;
 171             fix1             = 0;
 172             fiy1             = 0;
 173             fiz1             = 0;
 174             fix2             = 0;
 175             fiy2             = 0;
 176             fiz2             = 0;
 177             fix3             = 0;
 178             fiy3             = 0;
 179             fiz3             = 0;
 180             fix4             = 0;
 181             fiy4             = 0;
 182             fiz4             = 0;
 183
 184             for(k=nj0; (k<nj1); k++)
 185             {
 186                 jnr              = jjnr[k];
 187                 weight_cg2       = wf[jnr];
 188                 weight_product   = weight_cg1*weight_cg2;
 189                 if (weight_product < ALMOST_ZERO) {
 190                        hybscal = 1.0;
 191                 }
 192                 else if (weight_product >= ALMOST_ONE)
 193                 {
 194                   /* force is zero, skip this molecule */
 195                        continue;
 196                 }
 197                 else
 198                 {
 199                    hybscal = 1.0 - weight_product;
 200                 }
 201                 j3               = 3*jnr;
 202                 jx1              = pos[j3+0];
 203                 jy1              = pos[j3+1];
 204                 jz1              = pos[j3+2];
 205                 jx2              = pos[j3+3];
 206                 jy2              = pos[j3+4];
 207                 jz2              = pos[j3+5];
 208                 jx3              = pos[j3+6];
 209                 jy3              = pos[j3+7];
 210                 jz3              = pos[j3+8];
 211                 jx4              = pos[j3+9];
 212                 jy4              = pos[j3+10];
 213                 jz4              = pos[j3+11];
 214                 dx11             = ix1 - jx1;
 215                 dy11             = iy1 - jy1;
 216                 dz11             = iz1 - jz1;
 217                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 218                 dx22             = ix2 - jx2;
 219                 dy22             = iy2 - jy2;
 220                 dz22             = iz2 - jz2;
 221                 rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 222                 dx23             = ix2 - jx3;
 223                 dy23             = iy2 - jy3;
 224                 dz23             = iz2 - jz3;
 225                 rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 226                 dx24             = ix2 - jx4;
 227                 dy24             = iy2 - jy4;
 228                 dz24             = iz2 - jz4;
 229                 rsq24            = dx24*dx24+dy24*dy24+dz24*dz24;
 230                 dx32             = ix3 - jx2;
 231                 dy32             = iy3 - jy2;
 232                 dz32             = iz3 - jz2;
 233                 rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 234                 dx33             = ix3 - jx3;
 235                 dy33             = iy3 - jy3;
 236                 dz33             = iz3 - jz3;
 237                 rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 238                 dx34             = ix3 - jx4;
 239                 dy34             = iy3 - jy4;
 240                 dz34             = iz3 - jz4;
 241                 rsq34            = dx34*dx34+dy34*dy34+dz34*dz34;
 242                 dx42             = ix4 - jx2;
 243                 dy42             = iy4 - jy2;
 244                 dz42             = iz4 - jz2;
 245                 rsq42            = dx42*dx42+dy42*dy42+dz42*dz42;
 246                 dx43             = ix4 - jx3;
 247                 dy43             = iy4 - jy3;
 248                 dz43             = iz4 - jz3;
 249                 rsq43            = dx43*dx43+dy43*dy43+dz43*dz43;
 250                 dx44             = ix4 - jx4;
 251                 dy44             = iy4 - jy4;
 252                 dz44             = iz4 - jz4;
 253                 rsq44            = dx44*dx44+dy44*dy44+dz44*dz44;
 254                 rinvsq           = 1.0/rsq11;
 255                 rinv22           = 1.0/sqrt(rsq22);
 256                 rinv23           = 1.0/sqrt(rsq23);
 257                 rinv24           = 1.0/sqrt(rsq24);
 258                 rinv32           = 1.0/sqrt(rsq32);
 259                 rinv33           = 1.0/sqrt(rsq33);
 260                 rinv34           = 1.0/sqrt(rsq34);
 261                 rinv42           = 1.0/sqrt(rsq42);
 262                 rinv43           = 1.0/sqrt(rsq43);
 263                 rinv44           = 1.0/sqrt(rsq44);
 264                 rinvsix          = rinvsq*rinvsq*rinvsq;
 265                 Vvdw6            = c6*rinvsix;
 266                 Vvdw12           = c12*rinvsix*rinvsix;
 267                 Vvdwtot          = Vvdwtot+Vvdw12-Vvdw6;
 268                 fscal            = (12.0*Vvdw12-6.0*Vvdw6)*rinvsq;
 269                 fscal *= hybscal;
 270                 tx               = fscal*dx11;
 271                 ty               = fscal*dy11;
 272                 tz               = fscal*dz11;
 273                 fix1             = fix1 + tx;
 274                 fiy1             = fiy1 + ty;
 275                 fiz1             = fiz1 + tz;
 276                 faction[j3+0]    = faction[j3+0] - tx;
 277                 faction[j3+1]    = faction[j3+1] - ty;
 278                 faction[j3+2]    = faction[j3+2] - tz;
 279                 qq               = qqHH;
 280                 r                = rsq22*rinv22;
 281                 rt               = r*tabscale;
 282                 n0               = rt;
 283                 eps              = rt-n0;
 284                 eps2             = eps*eps;
 285                 nnn              = 4*n0;
 286                 Y                = VFtab[nnn];
 287                 F                = VFtab[nnn+1];
 288                 Geps             = eps*VFtab[nnn+2];
 289                 Heps2            = eps2*VFtab[nnn+3];
 290                 Fp               = F+Geps+Heps2;
 291                 VV               = Y+eps*Fp;
 292                 FF               = Fp+Geps+2.0*Heps2;
 293                 vcoul            = qq*VV;
 294                 fijC             = qq*FF;
 295                 vctot            = vctot + vcoul;
 296                 fscal            = -((fijC)*tabscale)*rinv22;
 297                 fscal *= hybscal;
 298                 tx               = fscal*dx22;
 299                 ty               = fscal*dy22;
 300                 tz               = fscal*dz22;
 301                 fix2             = fix2 + tx;
 302                 fiy2             = fiy2 + ty;
 303                 fiz2             = fiz2 + tz;
 304                 fjx2             = faction[j3+3] - tx;
 305                 fjy2             = faction[j3+4] - ty;
 306                 fjz2             = faction[j3+5] - tz;
 307                 qq               = qqHH;
 308                 r                = rsq23*rinv23;
 309                 rt               = r*tabscale;
 310                 n0               = rt;
 311                 eps              = rt-n0;
 312                 eps2             = eps*eps;
 313                 nnn              = 4*n0;
 314                 Y                = VFtab[nnn];
 315                 F                = VFtab[nnn+1];
 316                 Geps             = eps*VFtab[nnn+2];
 317                 Heps2            = eps2*VFtab[nnn+3];
 318                 Fp               = F+Geps+Heps2;
 319                 VV               = Y+eps*Fp;
 320                 FF               = Fp+Geps+2.0*Heps2;
 321                 vcoul            = qq*VV;
 322                 fijC             = qq*FF;
 323                 vctot            = vctot + vcoul;
 324                 fscal            = -((fijC)*tabscale)*rinv23;
 325                 fscal *= hybscal;
 326                 tx               = fscal*dx23;
 327                 ty               = fscal*dy23;
 328                 tz               = fscal*dz23;
 329                 fix2             = fix2 + tx;
 330                 fiy2             = fiy2 + ty;
 331                 fiz2             = fiz2 + tz;
 332                 fjx3             = faction[j3+6] - tx;
 333                 fjy3             = faction[j3+7] - ty;
 334                 fjz3             = faction[j3+8] - tz;
 335                 qq               = qqMH;
 336                 r                = rsq24*rinv24;
 337                 rt               = r*tabscale;
 338                 n0               = rt;
 339                 eps              = rt-n0;
 340                 eps2             = eps*eps;
 341                 nnn              = 4*n0;
 342                 Y                = VFtab[nnn];
 343                 F                = VFtab[nnn+1];
 344                 Geps             = eps*VFtab[nnn+2];
 345                 Heps2            = eps2*VFtab[nnn+3];
 346                 Fp               = F+Geps+Heps2;
 347                 VV               = Y+eps*Fp;
 348                 FF               = Fp+Geps+2.0*Heps2;
 349                 vcoul            = qq*VV;
 350                 fijC             = qq*FF;
 351                 vctot            = vctot + vcoul;
 352                 fscal            = -((fijC)*tabscale)*rinv24;
 353                 fscal *= hybscal;
 354                 tx               = fscal*dx24;
 355                 ty               = fscal*dy24;
 356                 tz               = fscal*dz24;
 357                 fix2             = fix2 + tx;
 358                 fiy2             = fiy2 + ty;
 359                 fiz2             = fiz2 + tz;
 360                 fjx4             = faction[j3+9] - tx;
 361                 fjy4             = faction[j3+10] - ty;
 362                 fjz4             = faction[j3+11] - tz;
 363                 qq               = qqHH;
 364                 r                = rsq32*rinv32;
 365                 rt               = r*tabscale;
 366                 n0               = rt;
 367                 eps              = rt-n0;
 368                 eps2             = eps*eps;
 369                 nnn              = 4*n0;
 370                 Y                = VFtab[nnn];
 371                 F                = VFtab[nnn+1];
 372                 Geps             = eps*VFtab[nnn+2];
 373                 Heps2            = eps2*VFtab[nnn+3];
 374                 Fp               = F+Geps+Heps2;
 375                 VV               = Y+eps*Fp;
 376                 FF               = Fp+Geps+2.0*Heps2;
 377                 vcoul            = qq*VV;
 378                 fijC             = qq*FF;
 379                 vctot            = vctot + vcoul;
 380                 fscal            = -((fijC)*tabscale)*rinv32;
 381                 fscal *= hybscal;
 382                 tx               = fscal*dx32;
 383                 ty               = fscal*dy32;
 384                 tz               = fscal*dz32;
 385                 fix3             = fix3 + tx;
 386                 fiy3             = fiy3 + ty;
 387                 fiz3             = fiz3 + tz;
 388                 fjx2             = fjx2 - tx;
 389                 fjy2             = fjy2 - ty;
 390                 fjz2             = fjz2 - tz;
 391                 qq               = qqHH;
 392                 r                = rsq33*rinv33;
 393                 rt               = r*tabscale;
 394                 n0               = rt;
 395                 eps              = rt-n0;
 396                 eps2             = eps*eps;
 397                 nnn              = 4*n0;
 398                 Y                = VFtab[nnn];
 399                 F                = VFtab[nnn+1];
 400                 Geps             = eps*VFtab[nnn+2];
 401                 Heps2            = eps2*VFtab[nnn+3];
 402                 Fp               = F+Geps+Heps2;
 403                 VV               = Y+eps*Fp;
 404                 FF               = Fp+Geps+2.0*Heps2;
 405                 vcoul            = qq*VV;
 406                 fijC             = qq*FF;
 407                 vctot            = vctot + vcoul;
 408                 fscal            = -((fijC)*tabscale)*rinv33;
 409                 fscal *= hybscal;
 410                 tx               = fscal*dx33;
 411                 ty               = fscal*dy33;
 412                 tz               = fscal*dz33;
 413                 fix3             = fix3 + tx;
 414                 fiy3             = fiy3 + ty;
 415                 fiz3             = fiz3 + tz;
 416                 fjx3             = fjx3 - tx;
 417                 fjy3             = fjy3 - ty;
 418                 fjz3             = fjz3 - tz;
 419                 qq               = qqMH;
 420                 r                = rsq34*rinv34;
 421                 rt               = r*tabscale;
 422                 n0               = rt;
 423                 eps              = rt-n0;
 424                 eps2             = eps*eps;
 425                 nnn              = 4*n0;
 426                 Y                = VFtab[nnn];
 427                 F                = VFtab[nnn+1];
 428                 Geps             = eps*VFtab[nnn+2];
 429                 Heps2            = eps2*VFtab[nnn+3];
 430                 Fp               = F+Geps+Heps2;
 431                 VV               = Y+eps*Fp;
 432                 FF               = Fp+Geps+2.0*Heps2;
 433                 vcoul            = qq*VV;
 434                 fijC             = qq*FF;
 435                 vctot            = vctot + vcoul;
 436                 fscal            = -((fijC)*tabscale)*rinv34;
 437                 fscal *= hybscal;
 438                 tx               = fscal*dx34;
 439                 ty               = fscal*dy34;
 440                 tz               = fscal*dz34;
 441                 fix3             = fix3 + tx;
 442                 fiy3             = fiy3 + ty;
 443                 fiz3             = fiz3 + tz;
 444                 fjx4             = fjx4 - tx;
 445                 fjy4             = fjy4 - ty;
 446                 fjz4             = fjz4 - tz;
 447                 qq               = qqMH;
 448                 r                = rsq42*rinv42;
 449                 rt               = r*tabscale;
 450                 n0               = rt;
 451                 eps              = rt-n0;
 452                 eps2             = eps*eps;
 453                 nnn              = 4*n0;
 454                 Y                = VFtab[nnn];
 455                 F                = VFtab[nnn+1];
 456                 Geps             = eps*VFtab[nnn+2];
 457                 Heps2            = eps2*VFtab[nnn+3];
 458                 Fp               = F+Geps+Heps2;
 459                 VV               = Y+eps*Fp;
 460                 FF               = Fp+Geps+2.0*Heps2;
 461                 vcoul            = qq*VV;
 462                 fijC             = qq*FF;
 463                 vctot            = vctot + vcoul;
 464                 fscal            = -((fijC)*tabscale)*rinv42;
 465                 fscal *= hybscal;
 466                 tx               = fscal*dx42;
 467                 ty               = fscal*dy42;
 468                 tz               = fscal*dz42;
 469                 fix4             = fix4 + tx;
 470                 fiy4             = fiy4 + ty;
 471                 fiz4             = fiz4 + tz;
 472                 faction[j3+3]    = fjx2 - tx;
 473                 faction[j3+4]    = fjy2 - ty;
 474                 faction[j3+5]    = fjz2 - tz;
 475                 qq               = qqMH;
 476                 r                = rsq43*rinv43;
 477                 rt               = r*tabscale;
 478                 n0               = rt;
 479                 eps              = rt-n0;
 480                 eps2             = eps*eps;
 481                 nnn              = 4*n0;
 482                 Y                = VFtab[nnn];
 483                 F                = VFtab[nnn+1];
 484                 Geps             = eps*VFtab[nnn+2];
 485                 Heps2            = eps2*VFtab[nnn+3];
 486                 Fp               = F+Geps+Heps2;
 487                 VV               = Y+eps*Fp;
 488                 FF               = Fp+Geps+2.0*Heps2;
 489                 vcoul            = qq*VV;
 490                 fijC             = qq*FF;
 491                 vctot            = vctot + vcoul;
 492                 fscal            = -((fijC)*tabscale)*rinv43;
 493                 fscal *= hybscal;
 494                 tx               = fscal*dx43;
 495                 ty               = fscal*dy43;
 496                 tz               = fscal*dz43;
 497                 fix4             = fix4 + tx;
 498                 fiy4             = fiy4 + ty;
 499                 fiz4             = fiz4 + tz;
 500                 faction[j3+6]    = fjx3 - tx;
 501                 faction[j3+7]    = fjy3 - ty;
 502                 faction[j3+8]    = fjz3 - tz;
 503                 qq               = qqMM;
 504                 r                = rsq44*rinv44;
 505                 rt               = r*tabscale;
 506                 n0               = rt;
 507                 eps              = rt-n0;
 508                 eps2             = eps*eps;
 509                 nnn              = 4*n0;
 510                 Y                = VFtab[nnn];
 511                 F                = VFtab[nnn+1];
 512                 Geps             = eps*VFtab[nnn+2];
 513                 Heps2            = eps2*VFtab[nnn+3];
 514                 Fp               = F+Geps+Heps2;
 515                 VV               = Y+eps*Fp;
 516                 FF               = Fp+Geps+2.0*Heps2;
 517                 vcoul            = qq*VV;
 518                 fijC             = qq*FF;
 519                 vctot            = vctot + vcoul;
 520                 fscal            = -((fijC)*tabscale)*rinv44;
 521                 fscal *= hybscal;
 522                 tx               = fscal*dx44;
 523                 ty               = fscal*dy44;
 524                 tz               = fscal*dz44;
 525                 fix4             = fix4 + tx;
 526                 fiy4             = fiy4 + ty;
 527                 fiz4             = fiz4 + tz;
 528                 faction[j3+9]    = fjx4 - tx;
 529                 faction[j3+10]   = fjy4 - ty;
 530                 faction[j3+11]   = fjz4 - tz;
 531             }
 532
 533             faction[ii3+0]   = faction[ii3+0] + fix1;
 534             faction[ii3+1]   = faction[ii3+1] + fiy1;
 535             faction[ii3+2]   = faction[ii3+2] + fiz1;
 536             faction[ii3+3]   = faction[ii3+3] + fix2;
 537             faction[ii3+4]   = faction[ii3+4] + fiy2;
 538             faction[ii3+5]   = faction[ii3+5] + fiz2;
 539             faction[ii3+6]   = faction[ii3+6] + fix3;
 540             faction[ii3+7]   = faction[ii3+7] + fiy3;
 541             faction[ii3+8]   = faction[ii3+8] + fiz3;
 542             faction[ii3+9]   = faction[ii3+9] + fix4;
 543             faction[ii3+10]  = faction[ii3+10] + fiy4;
 544             faction[ii3+11]  = faction[ii3+11] + fiz4;
 545             fshift[is3]      = fshift[is3]+fix1+fix2+fix3+fix4;
 546             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3+fiy4;
 547             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3+fiz4;
 548             ggid             = gid[n];
 549             Vc[ggid]         = Vc[ggid] + vctot;
 550             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 551             ninner           = ninner + nj1 - nj0;
 552         }
 553
 554         nouter           = nouter + nn1 - nn0;
 555     }
 556     while (nn1<nri);
 557
 558     *outeriter       = nouter;
 559     *inneriter       = ninner;
 560 }
 561
 562
 563
 564
 565
 566 /*
 567  * Gromacs nonbonded kernel nb_kernel314_adress_ex
 568  * Coulomb interaction:     Tabulated
 569  * VdW interaction:         Lennard-Jones
 570  * water optimization:      pairs of TIP4P interactions
 571  * Calculate forces:        yes
 572  */
 573 void nb_kernel314_adress_ex(
 574                     int *           p_nri,
 575                     int *           iinr,
 576                     int *           jindex,
 577                     int *           jjnr,
 578                     int *           shift,
 579                     real *         shiftvec,
 580                     real *         fshift,
 581                     int *           gid,
 582                     real *         pos,
 583                     real *         faction,
 584                     real *         charge,
 585                     real *         p_facel,
 586                     real *         p_krf,
 587                     real *         p_crf,
 588                     real *         Vc,
 589                     int *           type,
 590                     int *           p_ntype,
 591                     real *         vdwparam,
 592                     real *         Vvdw,
 593                     real *         p_tabscale,
 594                     real *         VFtab,
 595                     real *         invsqrta,
 596                     real *         dvda,
 597                     real *         p_gbtabscale,
 598                     real *         GBtab,
 599                     int *           p_nthreads,
 600                     int *           count,
 601                     void *          mtx,
 602                     int *           outeriter,
 603                     int *           inneriter,
 604                     real           force_cap,
 605                     real *         wf)
 606 {
 607     int           nri,ntype,nthreads;
 608     real         facel,krf,crf,tabscale,gbtabscale;
 609     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 610     int           nn0,nn1,nouter,ninner;
 611     real         shX,shY,shZ;
 612     real         fscal,tx,ty,tz;
 613     real         rinvsq;
 614     real         qq,vcoul,vctot;
 615     int           tj;
 616     real         rinvsix;
 617     real         Vvdw6,Vvdwtot;
 618     real         Vvdw12;
 619     real         r,rt,eps,eps2;
 620     int           n0,nnn;
 621     real         Y,F,Geps,Heps2,Fp,VV;
 622     real         FF;
 623     real         fijC;
 624     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 625     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
 626     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
 627     real         ix4,iy4,iz4,fix4,fiy4,fiz4;
 628     real         jx1,jy1,jz1;
 629     real         jx2,jy2,jz2,fjx2,fjy2,fjz2;
 630     real         jx3,jy3,jz3,fjx3,fjy3,fjz3;
 631     real         jx4,jy4,jz4,fjx4,fjy4,fjz4;
 632     real         dx11,dy11,dz11,rsq11;
 633     real         dx22,dy22,dz22,rsq22,rinv22;
 634     real         dx23,dy23,dz23,rsq23,rinv23;
 635     real         dx24,dy24,dz24,rsq24,rinv24;
 636     real         dx32,dy32,dz32,rsq32,rinv32;
 637     real         dx33,dy33,dz33,rsq33,rinv33;
 638     real         dx34,dy34,dz34,rsq34,rinv34;
 639     real         dx42,dy42,dz42,rsq42,rinv42;
 640     real         dx43,dy43,dz43,rsq43,rinv43;
 641     real         dx44,dy44,dz44,rsq44,rinv44;
 642     real         qH,qM,qqMM,qqMH,qqHH;
 643     real         c6,c12;
 644     real         weight_cg1, weight_cg2, weight_product;
 645     real         hybscal;
 646
 647     nri              = *p_nri;
 648     ntype            = *p_ntype;
 649     nthreads         = *p_nthreads;
 650     facel            = *p_facel;
 651     krf              = *p_krf;
 652     crf              = *p_crf;
 653     tabscale         = *p_tabscale;
 654     ii               = iinr[0];
 655     qH               = charge[ii+1];
 656     qM               = charge[ii+3];
 657     qqMM             = facel*qM*qM;
 658     qqMH             = facel*qM*qH;
 659     qqHH             = facel*qH*qH;
 660     tj               = 2*(ntype+1)*type[ii];
 661     c6               = vdwparam[tj];
 662     c12              = vdwparam[tj+1];
 663
 664     nouter           = 0;
 665     ninner           = 0;
 666
 667     do
 668     {
 669         #ifdef GMX_THREAD_SHM_FDECOMP
 670         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 671         nn0              = *count;
 672         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 673         *count           = nn1;
 674         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 675         if(nn1>nri) nn1=nri;
 676         #else
 677         nn0 = 0;
 678         nn1 = nri;
 679         #endif
 680
 681         for(n=nn0; (n<nn1); n++)
 682         {
 683             is3              = 3*shift[n];
 684             shX              = shiftvec[is3];
 685             shY              = shiftvec[is3+1];
 686             shZ              = shiftvec[is3+2];
 687             nj0              = jindex[n];
 688             nj1              = jindex[n+1];
 689             ii               = iinr[n];
 690             ii3              = 3*ii;
 691             ix1              = shX + pos[ii3+0];
 692             iy1              = shY + pos[ii3+1];
 693             iz1              = shZ + pos[ii3+2];
 694             ix2              = shX + pos[ii3+3];
 695             iy2              = shY + pos[ii3+4];
 696             iz2              = shZ + pos[ii3+5];
 697             ix3              = shX + pos[ii3+6];
 698             iy3              = shY + pos[ii3+7];
 699             iz3              = shZ + pos[ii3+8];
 700             ix4              = shX + pos[ii3+9];
 701             iy4              = shY + pos[ii3+10];
 702             iz4              = shZ + pos[ii3+11];
 703             weight_cg1       = wf[ii];
 704             vctot            = 0;
 705             Vvdwtot          = 0;
 706             fix1             = 0;
 707             fiy1             = 0;
 708             fiz1             = 0;
 709             fix2             = 0;
 710             fiy2             = 0;
 711             fiz2             = 0;
 712             fix3             = 0;
 713             fiy3             = 0;
 714             fiz3             = 0;
 715             fix4             = 0;
 716             fiy4             = 0;
 717             fiz4             = 0;
 718
 719             for(k=nj0; (k<nj1); k++)
 720             {
 721                 jnr              = jjnr[k];
 722                 weight_cg2       = wf[jnr];
 723                 weight_product   = weight_cg1*weight_cg2;
 724                 if (weight_product < ALMOST_ZERO) {
 725                 /* force is zero, skip this molecule */
 726                  continue;
 727                 }
 728                 else if (weight_product >= ALMOST_ONE)
 729                 {
 730                        hybscal = 1.0;
 731                 }
 732                 else
 733                 {
 734                    hybscal = weight_product;
 735                 }
 736                 j3               = 3*jnr;
 737                 jx1              = pos[j3+0];
 738                 jy1              = pos[j3+1];
 739                 jz1              = pos[j3+2];
 740                 jx2              = pos[j3+3];
 741                 jy2              = pos[j3+4];
 742                 jz2              = pos[j3+5];
 743                 jx3              = pos[j3+6];
 744                 jy3              = pos[j3+7];
 745                 jz3              = pos[j3+8];
 746                 jx4              = pos[j3+9];
 747                 jy4              = pos[j3+10];
 748                 jz4              = pos[j3+11];
 749                 dx11             = ix1 - jx1;
 750                 dy11             = iy1 - jy1;
 751                 dz11             = iz1 - jz1;
 752                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 753                 dx22             = ix2 - jx2;
 754                 dy22             = iy2 - jy2;
 755                 dz22             = iz2 - jz2;
 756                 rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 757                 dx23             = ix2 - jx3;
 758                 dy23             = iy2 - jy3;
 759                 dz23             = iz2 - jz3;
 760                 rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 761                 dx24             = ix2 - jx4;
 762                 dy24             = iy2 - jy4;
 763                 dz24             = iz2 - jz4;
 764                 rsq24            = dx24*dx24+dy24*dy24+dz24*dz24;
 765                 dx32             = ix3 - jx2;
 766                 dy32             = iy3 - jy2;
 767                 dz32             = iz3 - jz2;
 768                 rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 769                 dx33             = ix3 - jx3;
 770                 dy33             = iy3 - jy3;
 771                 dz33             = iz3 - jz3;
 772                 rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 773                 dx34             = ix3 - jx4;
 774                 dy34             = iy3 - jy4;
 775                 dz34             = iz3 - jz4;
 776                 rsq34            = dx34*dx34+dy34*dy34+dz34*dz34;
 777                 dx42             = ix4 - jx2;
 778                 dy42             = iy4 - jy2;
 779                 dz42             = iz4 - jz2;
 780                 rsq42            = dx42*dx42+dy42*dy42+dz42*dz42;
 781                 dx43             = ix4 - jx3;
 782                 dy43             = iy4 - jy3;
 783                 dz43             = iz4 - jz3;
 784                 rsq43            = dx43*dx43+dy43*dy43+dz43*dz43;
 785                 dx44             = ix4 - jx4;
 786                 dy44             = iy4 - jy4;
 787                 dz44             = iz4 - jz4;
 788                 rsq44            = dx44*dx44+dy44*dy44+dz44*dz44;
 789                 rinvsq           = 1.0/rsq11;
 790                 rinv22           = 1.0/sqrt(rsq22);
 791                 rinv23           = 1.0/sqrt(rsq23);
 792                 rinv24           = 1.0/sqrt(rsq24);
 793                 rinv32           = 1.0/sqrt(rsq32);
 794                 rinv33           = 1.0/sqrt(rsq33);
 795                 rinv34           = 1.0/sqrt(rsq34);
 796                 rinv42           = 1.0/sqrt(rsq42);
 797                 rinv43           = 1.0/sqrt(rsq43);
 798                 rinv44           = 1.0/sqrt(rsq44);
 799                 rinvsix          = rinvsq*rinvsq*rinvsq;
 800                 Vvdw6            = c6*rinvsix;
 801                 Vvdw12           = c12*rinvsix*rinvsix;
 802                 Vvdwtot          = Vvdwtot+Vvdw12-Vvdw6;
 803                 fscal            = (12.0*Vvdw12-6.0*Vvdw6)*rinvsq;
 804                 fscal *= hybscal;
 805                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 806                 fscal=force_cap*fscal/fabs(fscal);
 807                 }
 808                 tx               = fscal*dx11;
 809                 ty               = fscal*dy11;
 810                 tz               = fscal*dz11;
 811                 fix1             = fix1 + tx;
 812                 fiy1             = fiy1 + ty;
 813                 fiz1             = fiz1 + tz;
 814                 faction[j3+0]    = faction[j3+0] - tx;
 815                 faction[j3+1]    = faction[j3+1] - ty;
 816                 faction[j3+2]    = faction[j3+2] - tz;
 817                 qq               = qqHH;
 818                 r                = rsq22*rinv22;
 819                 rt               = r*tabscale;
 820                 n0               = rt;
 821                 eps              = rt-n0;
 822                 eps2             = eps*eps;
 823                 nnn              = 4*n0;
 824                 Y                = VFtab[nnn];
 825                 F                = VFtab[nnn+1];
 826                 Geps             = eps*VFtab[nnn+2];
 827                 Heps2            = eps2*VFtab[nnn+3];
 828                 Fp               = F+Geps+Heps2;
 829                 VV               = Y+eps*Fp;
 830                 FF               = Fp+Geps+2.0*Heps2;
 831                 vcoul            = qq*VV;
 832                 fijC             = qq*FF;
 833                 vctot            = vctot + vcoul;
 834                 fscal            = -((fijC)*tabscale)*rinv22;
 835                 fscal *= hybscal;
 836                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 837                 fscal=force_cap*fscal/fabs(fscal);
 838                 }
 839                 tx               = fscal*dx22;
 840                 ty               = fscal*dy22;
 841                 tz               = fscal*dz22;
 842                 fix2             = fix2 + tx;
 843                 fiy2             = fiy2 + ty;
 844                 fiz2             = fiz2 + tz;
 845                 fjx2             = faction[j3+3] - tx;
 846                 fjy2             = faction[j3+4] - ty;
 847                 fjz2             = faction[j3+5] - tz;
 848                 qq               = qqHH;
 849                 r                = rsq23*rinv23;
 850                 rt               = r*tabscale;
 851                 n0               = rt;
 852                 eps              = rt-n0;
 853                 eps2             = eps*eps;
 854                 nnn              = 4*n0;
 855                 Y                = VFtab[nnn];
 856                 F                = VFtab[nnn+1];
 857                 Geps             = eps*VFtab[nnn+2];
 858                 Heps2            = eps2*VFtab[nnn+3];
 859                 Fp               = F+Geps+Heps2;
 860                 VV               = Y+eps*Fp;
 861                 FF               = Fp+Geps+2.0*Heps2;
 862                 vcoul            = qq*VV;
 863                 fijC             = qq*FF;
 864                 vctot            = vctot + vcoul;
 865                 fscal            = -((fijC)*tabscale)*rinv23;
 866                 fscal *= hybscal;
 867                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 868                 fscal=force_cap*fscal/fabs(fscal);
 869                 }
 870                 tx               = fscal*dx23;
 871                 ty               = fscal*dy23;
 872                 tz               = fscal*dz23;
 873                 fix2             = fix2 + tx;
 874                 fiy2             = fiy2 + ty;
 875                 fiz2             = fiz2 + tz;
 876                 fjx3             = faction[j3+6] - tx;
 877                 fjy3             = faction[j3+7] - ty;
 878                 fjz3             = faction[j3+8] - tz;
 879                 qq               = qqMH;
 880                 r                = rsq24*rinv24;
 881                 rt               = r*tabscale;
 882                 n0               = rt;
 883                 eps              = rt-n0;
 884                 eps2             = eps*eps;
 885                 nnn              = 4*n0;
 886                 Y                = VFtab[nnn];
 887                 F                = VFtab[nnn+1];
 888                 Geps             = eps*VFtab[nnn+2];
 889                 Heps2            = eps2*VFtab[nnn+3];
 890                 Fp               = F+Geps+Heps2;
 891                 VV               = Y+eps*Fp;
 892                 FF               = Fp+Geps+2.0*Heps2;
 893                 vcoul            = qq*VV;
 894                 fijC             = qq*FF;
 895                 vctot            = vctot + vcoul;
 896                 fscal            = -((fijC)*tabscale)*rinv24;
 897                 fscal *= hybscal;
 898                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 899                 fscal=force_cap*fscal/fabs(fscal);
 900                 }
 901                 tx               = fscal*dx24;
 902                 ty               = fscal*dy24;
 903                 tz               = fscal*dz24;
 904                 fix2             = fix2 + tx;
 905                 fiy2             = fiy2 + ty;
 906                 fiz2             = fiz2 + tz;
 907                 fjx4             = faction[j3+9] - tx;
 908                 fjy4             = faction[j3+10] - ty;
 909                 fjz4             = faction[j3+11] - tz;
 910                 qq               = qqHH;
 911                 r                = rsq32*rinv32;
 912                 rt               = r*tabscale;
 913                 n0               = rt;
 914                 eps              = rt-n0;
 915                 eps2             = eps*eps;
 916                 nnn              = 4*n0;
 917                 Y                = VFtab[nnn];
 918                 F                = VFtab[nnn+1];
 919                 Geps             = eps*VFtab[nnn+2];
 920                 Heps2            = eps2*VFtab[nnn+3];
 921                 Fp               = F+Geps+Heps2;
 922                 VV               = Y+eps*Fp;
 923                 FF               = Fp+Geps+2.0*Heps2;
 924                 vcoul            = qq*VV;
 925                 fijC             = qq*FF;
 926                 vctot            = vctot + vcoul;
 927                 fscal            = -((fijC)*tabscale)*rinv32;
 928                 fscal *= hybscal;
 929                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 930                 fscal=force_cap*fscal/fabs(fscal);
 931                 }
 932                 tx               = fscal*dx32;
 933                 ty               = fscal*dy32;
 934                 tz               = fscal*dz32;
 935                 fix3             = fix3 + tx;
 936                 fiy3             = fiy3 + ty;
 937                 fiz3             = fiz3 + tz;
 938                 fjx2             = fjx2 - tx;
 939                 fjy2             = fjy2 - ty;
 940                 fjz2             = fjz2 - tz;
 941                 qq               = qqHH;
 942                 r                = rsq33*rinv33;
 943                 rt               = r*tabscale;
 944                 n0               = rt;
 945                 eps              = rt-n0;
 946                 eps2             = eps*eps;
 947                 nnn              = 4*n0;
 948                 Y                = VFtab[nnn];
 949                 F                = VFtab[nnn+1];
 950                 Geps             = eps*VFtab[nnn+2];
 951                 Heps2            = eps2*VFtab[nnn+3];
 952                 Fp               = F+Geps+Heps2;
 953                 VV               = Y+eps*Fp;
 954                 FF               = Fp+Geps+2.0*Heps2;
 955                 vcoul            = qq*VV;
 956                 fijC             = qq*FF;
 957                 vctot            = vctot + vcoul;
 958                 fscal            = -((fijC)*tabscale)*rinv33;
 959                 fscal *= hybscal;
 960                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 961                 fscal=force_cap*fscal/fabs(fscal);
 962                 }
 963                 tx               = fscal*dx33;
 964                 ty               = fscal*dy33;
 965                 tz               = fscal*dz33;
 966                 fix3             = fix3 + tx;
 967                 fiy3             = fiy3 + ty;
 968                 fiz3             = fiz3 + tz;
 969                 fjx3             = fjx3 - tx;
 970                 fjy3             = fjy3 - ty;
 971                 fjz3             = fjz3 - tz;
 972                 qq               = qqMH;
 973                 r                = rsq34*rinv34;
 974                 rt               = r*tabscale;
 975                 n0               = rt;
 976                 eps              = rt-n0;
 977                 eps2             = eps*eps;
 978                 nnn              = 4*n0;
 979                 Y                = VFtab[nnn];
 980                 F                = VFtab[nnn+1];
 981                 Geps             = eps*VFtab[nnn+2];
 982                 Heps2            = eps2*VFtab[nnn+3];
 983                 Fp               = F+Geps+Heps2;
 984                 VV               = Y+eps*Fp;
 985                 FF               = Fp+Geps+2.0*Heps2;
 986                 vcoul            = qq*VV;
 987                 fijC             = qq*FF;
 988                 vctot            = vctot + vcoul;
 989                 fscal            = -((fijC)*tabscale)*rinv34;
 990                 fscal *= hybscal;
 991                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 992                 fscal=force_cap*fscal/fabs(fscal);
 993                 }
 994                 tx               = fscal*dx34;
 995                 ty               = fscal*dy34;
 996                 tz               = fscal*dz34;
 997                 fix3             = fix3 + tx;
 998                 fiy3             = fiy3 + ty;
 999                 fiz3             = fiz3 + tz;
1000                 fjx4             = fjx4 - tx;
1001                 fjy4             = fjy4 - ty;
1002                 fjz4             = fjz4 - tz;
1003                 qq               = qqMH;
1004                 r                = rsq42*rinv42;
1005                 rt               = r*tabscale;
1006                 n0               = rt;
1007                 eps              = rt-n0;
1008                 eps2             = eps*eps;
1009                 nnn              = 4*n0;
1010                 Y                = VFtab[nnn];
1011                 F                = VFtab[nnn+1];
1012                 Geps             = eps*VFtab[nnn+2];
1013                 Heps2            = eps2*VFtab[nnn+3];
1014                 Fp               = F+Geps+Heps2;
1015                 VV               = Y+eps*Fp;
1016                 FF               = Fp+Geps+2.0*Heps2;
1017                 vcoul            = qq*VV;
1018                 fijC             = qq*FF;
1019                 vctot            = vctot + vcoul;
1020                 fscal            = -((fijC)*tabscale)*rinv42;
1021                 fscal *= hybscal;
1022                 if(force_cap>0 && (fabs(fscal)> force_cap)){
1023                 fscal=force_cap*fscal/fabs(fscal);
1024                 }
1025                 tx               = fscal*dx42;
1026                 ty               = fscal*dy42;
1027                 tz               = fscal*dz42;
1028                 fix4             = fix4 + tx;
1029                 fiy4             = fiy4 + ty;
1030                 fiz4             = fiz4 + tz;
1031                 faction[j3+3]    = fjx2 - tx;
1032                 faction[j3+4]    = fjy2 - ty;
1033                 faction[j3+5]    = fjz2 - tz;
1034                 qq               = qqMH;
1035                 r                = rsq43*rinv43;
1036                 rt               = r*tabscale;
1037                 n0               = rt;
1038                 eps              = rt-n0;
1039                 eps2             = eps*eps;
1040                 nnn              = 4*n0;
1041                 Y                = VFtab[nnn];
1042                 F                = VFtab[nnn+1];
1043                 Geps             = eps*VFtab[nnn+2];
1044                 Heps2            = eps2*VFtab[nnn+3];
1045                 Fp               = F+Geps+Heps2;
1046                 VV               = Y+eps*Fp;
1047                 FF               = Fp+Geps+2.0*Heps2;
1048                 vcoul            = qq*VV;
1049                 fijC             = qq*FF;
1050                 vctot            = vctot + vcoul;
1051                 fscal            = -((fijC)*tabscale)*rinv43;
1052                 fscal *= hybscal;
1053                 if(force_cap>0 && (fabs(fscal)> force_cap)){
1054                 fscal=force_cap*fscal/fabs(fscal);
1055                 }
1056                 tx               = fscal*dx43;
1057                 ty               = fscal*dy43;
1058                 tz               = fscal*dz43;
1059                 fix4             = fix4 + tx;
1060                 fiy4             = fiy4 + ty;
1061                 fiz4             = fiz4 + tz;
1062                 faction[j3+6]    = fjx3 - tx;
1063                 faction[j3+7]    = fjy3 - ty;
1064                 faction[j3+8]    = fjz3 - tz;
1065                 qq               = qqMM;
1066                 r                = rsq44*rinv44;
1067                 rt               = r*tabscale;
1068                 n0               = rt;
1069                 eps              = rt-n0;
1070                 eps2             = eps*eps;
1071                 nnn              = 4*n0;
1072                 Y                = VFtab[nnn];
1073                 F                = VFtab[nnn+1];
1074                 Geps             = eps*VFtab[nnn+2];
1075                 Heps2            = eps2*VFtab[nnn+3];
1076                 Fp               = F+Geps+Heps2;
1077                 VV               = Y+eps*Fp;
1078                 FF               = Fp+Geps+2.0*Heps2;
1079                 vcoul            = qq*VV;
1080                 fijC             = qq*FF;
1081                 vctot            = vctot + vcoul;
1082                 fscal            = -((fijC)*tabscale)*rinv44;
1083                 fscal *= hybscal;
1084                 if(force_cap>0 && (fabs(fscal)> force_cap)){
1085                 fscal=force_cap*fscal/fabs(fscal);
1086                 }
1087                 tx               = fscal*dx44;
1088                 ty               = fscal*dy44;
1089                 tz               = fscal*dz44;
1090                 fix4             = fix4 + tx;
1091                 fiy4             = fiy4 + ty;
1092                 fiz4             = fiz4 + tz;
1093                 faction[j3+9]    = fjx4 - tx;
1094                 faction[j3+10]   = fjy4 - ty;
1095                 faction[j3+11]   = fjz4 - tz;
1096             }
1097
1098             faction[ii3+0]   = faction[ii3+0] + fix1;
1099             faction[ii3+1]   = faction[ii3+1] + fiy1;
1100             faction[ii3+2]   = faction[ii3+2] + fiz1;
1101             faction[ii3+3]   = faction[ii3+3] + fix2;
1102             faction[ii3+4]   = faction[ii3+4] + fiy2;
1103             faction[ii3+5]   = faction[ii3+5] + fiz2;
1104             faction[ii3+6]   = faction[ii3+6] + fix3;
1105             faction[ii3+7]   = faction[ii3+7] + fiy3;
1106             faction[ii3+8]   = faction[ii3+8] + fiz3;
1107             faction[ii3+9]   = faction[ii3+9] + fix4;
1108             faction[ii3+10]  = faction[ii3+10] + fiy4;
1109             faction[ii3+11]  = faction[ii3+11] + fiz4;
1110             fshift[is3]      = fshift[is3]+fix1+fix2+fix3+fix4;
1111             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3+fiy4;
1112             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3+fiz4;
1113             ggid             = gid[n];
1114             Vc[ggid]         = Vc[ggid] + vctot;
1115             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
1116             ninner           = ninner + nj1 - nj0;
1117         }
1118
1119         nouter           = nouter + nn1 - nn0;
1120     }
1121     while (nn1<nri);
1122
1123     *outeriter       = nouter;
1124     *inneriter       = ninner;
1125 }
1126
1127