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