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