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