src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel223_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_kernel223_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel223_adress_cg
  33  * Coulomb interaction:     Reaction field
  34  * VdW interaction:         Buckingham
  35  * water optimization:      TIP4P - other atoms
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel223_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         jq;
  80     real         qq,vcoul,vctot;
  81     int           nti;
  82     int           tj;
  83     real         rinvsix;
  84     real         Vvdw6,Vvdwtot;
  85     real         krsq;
  86     real         Vvdwexp,br;
  87     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
  88     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
  89     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
  90     real         ix4,iy4,iz4,fix4,fiy4,fiz4;
  91     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
  92     real         dx11,dy11,dz11,rsq11,rinv11;
  93     real         dx21,dy21,dz21,rsq21,rinv21;
  94     real         dx31,dy31,dz31,rsq31,rinv31;
  95     real         dx41,dy41,dz41,rsq41,rinv41;
  96     real         qH,qM;
  97     real         c6,cexp1,cexp2;
  98     real         weight_cg1, weight_cg2, weight_product;
  99     real         hybscal;
 100
 101     nri              = *p_nri;
 102     ntype            = *p_ntype;
 103     nthreads         = *p_nthreads;
 104     facel            = *p_facel;
 105     krf              = *p_krf;
 106     crf              = *p_crf;
 107     tabscale         = *p_tabscale;
 108     ii               = iinr[0];
 109     qH               = facel*charge[ii+1];
 110     qM               = facel*charge[ii+3];
 111     nti              = 3*ntype*type[ii];
 112
 113     nouter           = 0;
 114     ninner           = 0;
 115
 116     do
 117     {
 118         #ifdef GMX_THREAD_SHM_FDECOMP
 119         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 120         nn0              = *count;
 121         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 122         *count           = nn1;
 123         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 124         if(nn1>nri) nn1=nri;
 125         #else
 126         nn0 = 0;
 127         nn1 = nri;
 128         #endif
 129
 130         for(n=nn0; (n<nn1); n++)
 131         {
 132             is3              = 3*shift[n];
 133             shX              = shiftvec[is3];
 134             shY              = shiftvec[is3+1];
 135             shZ              = shiftvec[is3+2];
 136             nj0              = jindex[n];
 137             nj1              = jindex[n+1];
 138             ii               = iinr[n];
 139             ii3              = 3*ii;
 140             ix1              = shX + pos[ii3+0];
 141             iy1              = shY + pos[ii3+1];
 142             iz1              = shZ + pos[ii3+2];
 143             ix2              = shX + pos[ii3+3];
 144             iy2              = shY + pos[ii3+4];
 145             iz2              = shZ + pos[ii3+5];
 146             ix3              = shX + pos[ii3+6];
 147             iy3              = shY + pos[ii3+7];
 148             iz3              = shZ + pos[ii3+8];
 149             ix4              = shX + pos[ii3+9];
 150             iy4              = shY + pos[ii3+10];
 151             iz4              = shZ + pos[ii3+11];
 152             weight_cg1       = wf[ii];
 153             vctot            = 0;
 154             Vvdwtot          = 0;
 155             fix1             = 0;
 156             fiy1             = 0;
 157             fiz1             = 0;
 158             fix2             = 0;
 159             fiy2             = 0;
 160             fiz2             = 0;
 161             fix3             = 0;
 162             fiy3             = 0;
 163             fiz3             = 0;
 164             fix4             = 0;
 165             fiy4             = 0;
 166             fiz4             = 0;
 167
 168             for(k=nj0; (k<nj1); k++)
 169             {
 170                 jnr              = jjnr[k];
 171                 weight_cg2       = wf[jnr];
 172                 weight_product   = weight_cg1*weight_cg2;
 173                 if (weight_product < ALMOST_ZERO) {
 174                        hybscal = 1.0;
 175                 }
 176                 else if (weight_product >= ALMOST_ONE)
 177                 {
 178                   /* force is zero, skip this molecule */
 179                        continue;
 180                 }
 181                 else
 182                 {
 183                    hybscal = 1.0 - weight_product;
 184                 }
 185                 j3               = 3*jnr;
 186                 jx1              = pos[j3+0];
 187                 jy1              = pos[j3+1];
 188                 jz1              = pos[j3+2];
 189                 dx11             = ix1 - jx1;
 190                 dy11             = iy1 - jy1;
 191                 dz11             = iz1 - jz1;
 192                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 193                 dx21             = ix2 - jx1;
 194                 dy21             = iy2 - jy1;
 195                 dz21             = iz2 - jz1;
 196                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 197                 dx31             = ix3 - jx1;
 198                 dy31             = iy3 - jy1;
 199                 dz31             = iz3 - jz1;
 200                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 201                 dx41             = ix4 - jx1;
 202                 dy41             = iy4 - jy1;
 203                 dz41             = iz4 - jz1;
 204                 rsq41            = dx41*dx41+dy41*dy41+dz41*dz41;
 205                 rinv11           = 1.0/sqrt(rsq11);
 206                 rinv21           = 1.0/sqrt(rsq21);
 207                 rinv31           = 1.0/sqrt(rsq31);
 208                 rinv41           = 1.0/sqrt(rsq41);
 209                 tj               = nti+3*type[jnr];
 210                 c6               = vdwparam[tj];
 211                 cexp1            = vdwparam[tj+1];
 212                 cexp2            = vdwparam[tj+2];
 213                 rinvsq           = rinv11*rinv11;
 214                 rinvsix          = rinvsq*rinvsq*rinvsq;
 215                 Vvdw6            = c6*rinvsix;
 216                 br               = cexp2*rsq11*rinv11;
 217                 Vvdwexp          = cexp1*exp(-br);
 218                 Vvdwtot          = Vvdwtot+Vvdwexp-Vvdw6;
 219                 fscal            = (br*Vvdwexp-6.0*Vvdw6)*rinvsq;
 220                 fscal *= hybscal;
 221                 tx               = fscal*dx11;
 222                 ty               = fscal*dy11;
 223                 tz               = fscal*dz11;
 224                 fix1             = fix1 + tx;
 225                 fiy1             = fiy1 + ty;
 226                 fiz1             = fiz1 + tz;
 227                 fjx1             = faction[j3+0] - tx;
 228                 fjy1             = faction[j3+1] - ty;
 229                 fjz1             = faction[j3+2] - tz;
 230                 jq               = charge[jnr+0];
 231                 qq               = qH*jq;
 232                 rinvsq           = rinv21*rinv21;
 233                 krsq             = krf*rsq21;
 234                 vcoul            = qq*(rinv21+krsq-crf);
 235                 vctot            = vctot+vcoul;
 236                 fscal            = (qq*(rinv21-2.0*krsq))*rinvsq;
 237                 fscal *= hybscal;
 238                 tx               = fscal*dx21;
 239                 ty               = fscal*dy21;
 240                 tz               = fscal*dz21;
 241                 fix2             = fix2 + tx;
 242                 fiy2             = fiy2 + ty;
 243                 fiz2             = fiz2 + tz;
 244                 fjx1             = fjx1 - tx;
 245                 fjy1             = fjy1 - ty;
 246                 fjz1             = fjz1 - tz;
 247                 rinvsq           = rinv31*rinv31;
 248                 krsq             = krf*rsq31;
 249                 vcoul            = qq*(rinv31+krsq-crf);
 250                 vctot            = vctot+vcoul;
 251                 fscal            = (qq*(rinv31-2.0*krsq))*rinvsq;
 252                 fscal *= hybscal;
 253                 tx               = fscal*dx31;
 254                 ty               = fscal*dy31;
 255                 tz               = fscal*dz31;
 256                 fix3             = fix3 + tx;
 257                 fiy3             = fiy3 + ty;
 258                 fiz3             = fiz3 + tz;
 259                 fjx1             = fjx1 - tx;
 260                 fjy1             = fjy1 - ty;
 261                 fjz1             = fjz1 - tz;
 262                 qq               = qM*jq;
 263                 rinvsq           = rinv41*rinv41;
 264                 krsq             = krf*rsq41;
 265                 vcoul            = qq*(rinv41+krsq-crf);
 266                 vctot            = vctot+vcoul;
 267                 fscal            = (qq*(rinv41-2.0*krsq))*rinvsq;
 268                 fscal *= hybscal;
 269                 tx               = fscal*dx41;
 270                 ty               = fscal*dy41;
 271                 tz               = fscal*dz41;
 272                 fix4             = fix4 + tx;
 273                 fiy4             = fiy4 + ty;
 274                 fiz4             = fiz4 + tz;
 275                 faction[j3+0]    = fjx1 - tx;
 276                 faction[j3+1]    = fjy1 - ty;
 277                 faction[j3+2]    = fjz1 - tz;
 278             }
 279
 280             faction[ii3+0]   = faction[ii3+0] + fix1;
 281             faction[ii3+1]   = faction[ii3+1] + fiy1;
 282             faction[ii3+2]   = faction[ii3+2] + fiz1;
 283             faction[ii3+3]   = faction[ii3+3] + fix2;
 284             faction[ii3+4]   = faction[ii3+4] + fiy2;
 285             faction[ii3+5]   = faction[ii3+5] + fiz2;
 286             faction[ii3+6]   = faction[ii3+6] + fix3;
 287             faction[ii3+7]   = faction[ii3+7] + fiy3;
 288             faction[ii3+8]   = faction[ii3+8] + fiz3;
 289             faction[ii3+9]   = faction[ii3+9] + fix4;
 290             faction[ii3+10]  = faction[ii3+10] + fiy4;
 291             faction[ii3+11]  = faction[ii3+11] + fiz4;
 292             fshift[is3]      = fshift[is3]+fix1+fix2+fix3+fix4;
 293             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3+fiy4;
 294             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3+fiz4;
 295             ggid             = gid[n];
 296             Vc[ggid]         = Vc[ggid] + vctot;
 297             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 298             ninner           = ninner + nj1 - nj0;
 299         }
 300
 301         nouter           = nouter + nn1 - nn0;
 302     }
 303     while (nn1<nri);
 304
 305     *outeriter       = nouter;
 306     *inneriter       = ninner;
 307 }
 308
 309
 310
 311
 312
 313 /*
 314  * Gromacs nonbonded kernel nb_kernel223_adress_ex
 315  * Coulomb interaction:     Reaction field
 316  * VdW interaction:         Buckingham
 317  * water optimization:      TIP4P - other atoms
 318  * Calculate forces:        yes
 319  */
 320 void nb_kernel223_adress_ex(
 321                     int *           p_nri,
 322                     int *           iinr,
 323                     int *           jindex,
 324                     int *           jjnr,
 325                     int *           shift,
 326                     real *         shiftvec,
 327                     real *         fshift,
 328                     int *           gid,
 329                     real *         pos,
 330                     real *         faction,
 331                     real *         charge,
 332                     real *         p_facel,
 333                     real *         p_krf,
 334                     real *         p_crf,
 335                     real *         Vc,
 336                     int *           type,
 337                     int *           p_ntype,
 338                     real *         vdwparam,
 339                     real *         Vvdw,
 340                     real *         p_tabscale,
 341                     real *         VFtab,
 342                     real *         invsqrta,
 343                     real *         dvda,
 344                     real *         p_gbtabscale,
 345                     real *         GBtab,
 346                     int *           p_nthreads,
 347                     int *           count,
 348                     void *          mtx,
 349                     int *           outeriter,
 350                     int *           inneriter,
 351                     real           force_cap,
 352                     real *         wf)
 353 {
 354     int           nri,ntype,nthreads;
 355     real         facel,krf,crf,tabscale,gbtabscale;
 356     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 357     int           nn0,nn1,nouter,ninner;
 358     real         shX,shY,shZ;
 359     real         fscal,tx,ty,tz;
 360     real         rinvsq;
 361     real         jq;
 362     real         qq,vcoul,vctot;
 363     int           nti;
 364     int           tj;
 365     real         rinvsix;
 366     real         Vvdw6,Vvdwtot;
 367     real         krsq;
 368     real         Vvdwexp,br;
 369     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 370     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
 371     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
 372     real         ix4,iy4,iz4,fix4,fiy4,fiz4;
 373     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
 374     real         dx11,dy11,dz11,rsq11,rinv11;
 375     real         dx21,dy21,dz21,rsq21,rinv21;
 376     real         dx31,dy31,dz31,rsq31,rinv31;
 377     real         dx41,dy41,dz41,rsq41,rinv41;
 378     real         qH,qM;
 379     real         c6,cexp1,cexp2;
 380     real         weight_cg1, weight_cg2, weight_product;
 381     real         hybscal;
 382
 383     nri              = *p_nri;
 384     ntype            = *p_ntype;
 385     nthreads         = *p_nthreads;
 386     facel            = *p_facel;
 387     krf              = *p_krf;
 388     crf              = *p_crf;
 389     tabscale         = *p_tabscale;
 390     ii               = iinr[0];
 391     qH               = facel*charge[ii+1];
 392     qM               = facel*charge[ii+3];
 393     nti              = 3*ntype*type[ii];
 394
 395     nouter           = 0;
 396     ninner           = 0;
 397
 398     do
 399     {
 400         #ifdef GMX_THREAD_SHM_FDECOMP
 401         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 402         nn0              = *count;
 403         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 404         *count           = nn1;
 405         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 406         if(nn1>nri) nn1=nri;
 407         #else
 408         nn0 = 0;
 409         nn1 = nri;
 410         #endif
 411
 412         for(n=nn0; (n<nn1); n++)
 413         {
 414             is3              = 3*shift[n];
 415             shX              = shiftvec[is3];
 416             shY              = shiftvec[is3+1];
 417             shZ              = shiftvec[is3+2];
 418             nj0              = jindex[n];
 419             nj1              = jindex[n+1];
 420             ii               = iinr[n];
 421             ii3              = 3*ii;
 422             ix1              = shX + pos[ii3+0];
 423             iy1              = shY + pos[ii3+1];
 424             iz1              = shZ + pos[ii3+2];
 425             ix2              = shX + pos[ii3+3];
 426             iy2              = shY + pos[ii3+4];
 427             iz2              = shZ + pos[ii3+5];
 428             ix3              = shX + pos[ii3+6];
 429             iy3              = shY + pos[ii3+7];
 430             iz3              = shZ + pos[ii3+8];
 431             ix4              = shX + pos[ii3+9];
 432             iy4              = shY + pos[ii3+10];
 433             iz4              = shZ + pos[ii3+11];
 434             weight_cg1       = wf[ii];
 435             vctot            = 0;
 436             Vvdwtot          = 0;
 437             fix1             = 0;
 438             fiy1             = 0;
 439             fiz1             = 0;
 440             fix2             = 0;
 441             fiy2             = 0;
 442             fiz2             = 0;
 443             fix3             = 0;
 444             fiy3             = 0;
 445             fiz3             = 0;
 446             fix4             = 0;
 447             fiy4             = 0;
 448             fiz4             = 0;
 449
 450             for(k=nj0; (k<nj1); k++)
 451             {
 452                 jnr              = jjnr[k];
 453                 weight_cg2       = wf[jnr];
 454                 weight_product   = weight_cg1*weight_cg2;
 455                 if (weight_product < ALMOST_ZERO) {
 456                 /* force is zero, skip this molecule */
 457                  continue;
 458                 }
 459                 else if (weight_product >= ALMOST_ONE)
 460                 {
 461                        hybscal = 1.0;
 462                 }
 463                 else
 464                 {
 465                    hybscal = weight_product;
 466                 }
 467                 j3               = 3*jnr;
 468                 jx1              = pos[j3+0];
 469                 jy1              = pos[j3+1];
 470                 jz1              = pos[j3+2];
 471                 dx11             = ix1 - jx1;
 472                 dy11             = iy1 - jy1;
 473                 dz11             = iz1 - jz1;
 474                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 475                 dx21             = ix2 - jx1;
 476                 dy21             = iy2 - jy1;
 477                 dz21             = iz2 - jz1;
 478                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 479                 dx31             = ix3 - jx1;
 480                 dy31             = iy3 - jy1;
 481                 dz31             = iz3 - jz1;
 482                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 483                 dx41             = ix4 - jx1;
 484                 dy41             = iy4 - jy1;
 485                 dz41             = iz4 - jz1;
 486                 rsq41            = dx41*dx41+dy41*dy41+dz41*dz41;
 487                 rinv11           = 1.0/sqrt(rsq11);
 488                 rinv21           = 1.0/sqrt(rsq21);
 489                 rinv31           = 1.0/sqrt(rsq31);
 490                 rinv41           = 1.0/sqrt(rsq41);
 491                 tj               = nti+3*type[jnr];
 492                 c6               = vdwparam[tj];
 493                 cexp1            = vdwparam[tj+1];
 494                 cexp2            = vdwparam[tj+2];
 495                 rinvsq           = rinv11*rinv11;
 496                 rinvsix          = rinvsq*rinvsq*rinvsq;
 497                 Vvdw6            = c6*rinvsix;
 498                 br               = cexp2*rsq11*rinv11;
 499                 Vvdwexp          = cexp1*exp(-br);
 500                 Vvdwtot          = Vvdwtot+Vvdwexp-Vvdw6;
 501                 fscal            = (br*Vvdwexp-6.0*Vvdw6)*rinvsq;
 502                 fscal *= hybscal;
 503                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 504                 fscal=force_cap*fscal/fabs(fscal);
 505                 }
 506                 tx               = fscal*dx11;
 507                 ty               = fscal*dy11;
 508                 tz               = fscal*dz11;
 509                 fix1             = fix1 + tx;
 510                 fiy1             = fiy1 + ty;
 511                 fiz1             = fiz1 + tz;
 512                 fjx1             = faction[j3+0] - tx;
 513                 fjy1             = faction[j3+1] - ty;
 514                 fjz1             = faction[j3+2] - tz;
 515                 jq               = charge[jnr+0];
 516                 qq               = qH*jq;
 517                 rinvsq           = rinv21*rinv21;
 518                 krsq             = krf*rsq21;
 519                 vcoul            = qq*(rinv21+krsq-crf);
 520                 vctot            = vctot+vcoul;
 521                 fscal            = (qq*(rinv21-2.0*krsq))*rinvsq;
 522                 fscal *= hybscal;
 523                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 524                 fscal=force_cap*fscal/fabs(fscal);
 525                 }
 526                 tx               = fscal*dx21;
 527                 ty               = fscal*dy21;
 528                 tz               = fscal*dz21;
 529                 fix2             = fix2 + tx;
 530                 fiy2             = fiy2 + ty;
 531                 fiz2             = fiz2 + tz;
 532                 fjx1             = fjx1 - tx;
 533                 fjy1             = fjy1 - ty;
 534                 fjz1             = fjz1 - tz;
 535                 rinvsq           = rinv31*rinv31;
 536                 krsq             = krf*rsq31;
 537                 vcoul            = qq*(rinv31+krsq-crf);
 538                 vctot            = vctot+vcoul;
 539                 fscal            = (qq*(rinv31-2.0*krsq))*rinvsq;
 540                 fscal *= hybscal;
 541                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 542                 fscal=force_cap*fscal/fabs(fscal);
 543                 }
 544                 tx               = fscal*dx31;
 545                 ty               = fscal*dy31;
 546                 tz               = fscal*dz31;
 547                 fix3             = fix3 + tx;
 548                 fiy3             = fiy3 + ty;
 549                 fiz3             = fiz3 + tz;
 550                 fjx1             = fjx1 - tx;
 551                 fjy1             = fjy1 - ty;
 552                 fjz1             = fjz1 - tz;
 553                 qq               = qM*jq;
 554                 rinvsq           = rinv41*rinv41;
 555                 krsq             = krf*rsq41;
 556                 vcoul            = qq*(rinv41+krsq-crf);
 557                 vctot            = vctot+vcoul;
 558                 fscal            = (qq*(rinv41-2.0*krsq))*rinvsq;
 559                 fscal *= hybscal;
 560                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 561                 fscal=force_cap*fscal/fabs(fscal);
 562                 }
 563                 tx               = fscal*dx41;
 564                 ty               = fscal*dy41;
 565                 tz               = fscal*dz41;
 566                 fix4             = fix4 + tx;
 567                 fiy4             = fiy4 + ty;
 568                 fiz4             = fiz4 + tz;
 569                 faction[j3+0]    = fjx1 - tx;
 570                 faction[j3+1]    = fjy1 - ty;
 571                 faction[j3+2]    = fjz1 - tz;
 572             }
 573
 574             faction[ii3+0]   = faction[ii3+0] + fix1;
 575             faction[ii3+1]   = faction[ii3+1] + fiy1;
 576             faction[ii3+2]   = faction[ii3+2] + fiz1;
 577             faction[ii3+3]   = faction[ii3+3] + fix2;
 578             faction[ii3+4]   = faction[ii3+4] + fiy2;
 579             faction[ii3+5]   = faction[ii3+5] + fiz2;
 580             faction[ii3+6]   = faction[ii3+6] + fix3;
 581             faction[ii3+7]   = faction[ii3+7] + fiy3;
 582             faction[ii3+8]   = faction[ii3+8] + fiz3;
 583             faction[ii3+9]   = faction[ii3+9] + fix4;
 584             faction[ii3+10]  = faction[ii3+10] + fiy4;
 585             faction[ii3+11]  = faction[ii3+11] + fiz4;
 586             fshift[is3]      = fshift[is3]+fix1+fix2+fix3+fix4;
 587             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3+fiy4;
 588             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3+fiz4;
 589             ggid             = gid[n];
 590             Vc[ggid]         = Vc[ggid] + vctot;
 591             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 592             ninner           = ninner + nj1 - nj0;
 593         }
 594
 595         nouter           = nouter + nn1 - nn0;
 596     }
 597     while (nn1<nri);
 598
 599     *outeriter       = nouter;
 600     *inneriter       = ninner;
 601 }
 602
 603