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