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