src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel134_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_kernel134_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel134_adress_cg
  33  * Coulomb interaction:     Normal Coulomb
  34  * VdW interaction:         Tabulated
  35  * water optimization:      pairs of TIP4P interactions
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel134_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         qq,vcoul,vctot;
  80     int           tj;
  81     real         Vvdw6,Vvdwtot;
  82     real         Vvdw12;
  83     real         r,rt,eps,eps2;
  84     int           n0,nnn;
  85     real         Y,F,Geps,Heps2,Fp,VV;
  86     real         FF;
  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              = 8*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                 Vvdw6            = c6*VV;
 277                 fijD             = c6*FF;
 278                 nnn              = nnn+4;
 279                 Y                = VFtab[nnn];
 280                 F                = VFtab[nnn+1];
 281                 Geps             = eps*VFtab[nnn+2];
 282                 Heps2            = eps2*VFtab[nnn+3];
 283                 Fp               = F+Geps+Heps2;
 284                 VV               = Y+eps*Fp;
 285                 FF               = Fp+Geps+2.0*Heps2;
 286                 Vvdw12           = c12*VV;
 287                 fijR             = c12*FF;
 288                 Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12;
 289                 fscal            = -((fijD+fijR)*tabscale)*rinv11;
 290                 fscal *= hybscal;
 291                 tx               = fscal*dx11;
 292                 ty               = fscal*dy11;
 293                 tz               = fscal*dz11;
 294                 fix1             = fix1 + tx;
 295                 fiy1             = fiy1 + ty;
 296                 fiz1             = fiz1 + tz;
 297                 faction[j3+0]    = faction[j3+0] - tx;
 298                 faction[j3+1]    = faction[j3+1] - ty;
 299                 faction[j3+2]    = faction[j3+2] - tz;
 300                 qq               = qqHH;
 301                 rinvsq           = rinv22*rinv22;
 302                 vcoul            = qq*rinv22;
 303                 vctot            = vctot+vcoul;
 304                 fscal            = (vcoul)*rinvsq;
 305                 fscal *= hybscal;
 306                 tx               = fscal*dx22;
 307                 ty               = fscal*dy22;
 308                 tz               = fscal*dz22;
 309                 fix2             = fix2 + tx;
 310                 fiy2             = fiy2 + ty;
 311                 fiz2             = fiz2 + tz;
 312                 fjx2             = faction[j3+3] - tx;
 313                 fjy2             = faction[j3+4] - ty;
 314                 fjz2             = faction[j3+5] - tz;
 315                 qq               = qqHH;
 316                 rinvsq           = rinv23*rinv23;
 317                 vcoul            = qq*rinv23;
 318                 vctot            = vctot+vcoul;
 319                 fscal            = (vcoul)*rinvsq;
 320                 fscal *= hybscal;
 321                 tx               = fscal*dx23;
 322                 ty               = fscal*dy23;
 323                 tz               = fscal*dz23;
 324                 fix2             = fix2 + tx;
 325                 fiy2             = fiy2 + ty;
 326                 fiz2             = fiz2 + tz;
 327                 fjx3             = faction[j3+6] - tx;
 328                 fjy3             = faction[j3+7] - ty;
 329                 fjz3             = faction[j3+8] - tz;
 330                 qq               = qqMH;
 331                 rinvsq           = rinv24*rinv24;
 332                 vcoul            = qq*rinv24;
 333                 vctot            = vctot+vcoul;
 334                 fscal            = (vcoul)*rinvsq;
 335                 fscal *= hybscal;
 336                 tx               = fscal*dx24;
 337                 ty               = fscal*dy24;
 338                 tz               = fscal*dz24;
 339                 fix2             = fix2 + tx;
 340                 fiy2             = fiy2 + ty;
 341                 fiz2             = fiz2 + tz;
 342                 fjx4             = faction[j3+9] - tx;
 343                 fjy4             = faction[j3+10] - ty;
 344                 fjz4             = faction[j3+11] - tz;
 345                 qq               = qqHH;
 346                 rinvsq           = rinv32*rinv32;
 347                 vcoul            = qq*rinv32;
 348                 vctot            = vctot+vcoul;
 349                 fscal            = (vcoul)*rinvsq;
 350                 fscal *= hybscal;
 351                 tx               = fscal*dx32;
 352                 ty               = fscal*dy32;
 353                 tz               = fscal*dz32;
 354                 fix3             = fix3 + tx;
 355                 fiy3             = fiy3 + ty;
 356                 fiz3             = fiz3 + tz;
 357                 fjx2             = fjx2 - tx;
 358                 fjy2             = fjy2 - ty;
 359                 fjz2             = fjz2 - tz;
 360                 qq               = qqHH;
 361                 rinvsq           = rinv33*rinv33;
 362                 vcoul            = qq*rinv33;
 363                 vctot            = vctot+vcoul;
 364                 fscal            = (vcoul)*rinvsq;
 365                 fscal *= hybscal;
 366                 tx               = fscal*dx33;
 367                 ty               = fscal*dy33;
 368                 tz               = fscal*dz33;
 369                 fix3             = fix3 + tx;
 370                 fiy3             = fiy3 + ty;
 371                 fiz3             = fiz3 + tz;
 372                 fjx3             = fjx3 - tx;
 373                 fjy3             = fjy3 - ty;
 374                 fjz3             = fjz3 - tz;
 375                 qq               = qqMH;
 376                 rinvsq           = rinv34*rinv34;
 377                 vcoul            = qq*rinv34;
 378                 vctot            = vctot+vcoul;
 379                 fscal            = (vcoul)*rinvsq;
 380                 fscal *= hybscal;
 381                 tx               = fscal*dx34;
 382                 ty               = fscal*dy34;
 383                 tz               = fscal*dz34;
 384                 fix3             = fix3 + tx;
 385                 fiy3             = fiy3 + ty;
 386                 fiz3             = fiz3 + tz;
 387                 fjx4             = fjx4 - tx;
 388                 fjy4             = fjy4 - ty;
 389                 fjz4             = fjz4 - tz;
 390                 qq               = qqMH;
 391                 rinvsq           = rinv42*rinv42;
 392                 vcoul            = qq*rinv42;
 393                 vctot            = vctot+vcoul;
 394                 fscal            = (vcoul)*rinvsq;
 395                 fscal *= hybscal;
 396                 tx               = fscal*dx42;
 397                 ty               = fscal*dy42;
 398                 tz               = fscal*dz42;
 399                 fix4             = fix4 + tx;
 400                 fiy4             = fiy4 + ty;
 401                 fiz4             = fiz4 + tz;
 402                 faction[j3+3]    = fjx2 - tx;
 403                 faction[j3+4]    = fjy2 - ty;
 404                 faction[j3+5]    = fjz2 - tz;
 405                 qq               = qqMH;
 406                 rinvsq           = rinv43*rinv43;
 407                 vcoul            = qq*rinv43;
 408                 vctot            = vctot+vcoul;
 409                 fscal            = (vcoul)*rinvsq;
 410                 fscal *= hybscal;
 411                 tx               = fscal*dx43;
 412                 ty               = fscal*dy43;
 413                 tz               = fscal*dz43;
 414                 fix4             = fix4 + tx;
 415                 fiy4             = fiy4 + ty;
 416                 fiz4             = fiz4 + tz;
 417                 faction[j3+6]    = fjx3 - tx;
 418                 faction[j3+7]    = fjy3 - ty;
 419                 faction[j3+8]    = fjz3 - tz;
 420                 qq               = qqMM;
 421                 rinvsq           = rinv44*rinv44;
 422                 vcoul            = qq*rinv44;
 423                 vctot            = vctot+vcoul;
 424                 fscal            = (vcoul)*rinvsq;
 425                 fscal *= hybscal;
 426                 tx               = fscal*dx44;
 427                 ty               = fscal*dy44;
 428                 tz               = fscal*dz44;
 429                 fix4             = fix4 + tx;
 430                 fiy4             = fiy4 + ty;
 431                 fiz4             = fiz4 + tz;
 432                 faction[j3+9]    = fjx4 - tx;
 433                 faction[j3+10]   = fjy4 - ty;
 434                 faction[j3+11]   = fjz4 - tz;
 435             }
 436
 437             faction[ii3+0]   = faction[ii3+0] + fix1;
 438             faction[ii3+1]   = faction[ii3+1] + fiy1;
 439             faction[ii3+2]   = faction[ii3+2] + fiz1;
 440             faction[ii3+3]   = faction[ii3+3] + fix2;
 441             faction[ii3+4]   = faction[ii3+4] + fiy2;
 442             faction[ii3+5]   = faction[ii3+5] + fiz2;
 443             faction[ii3+6]   = faction[ii3+6] + fix3;
 444             faction[ii3+7]   = faction[ii3+7] + fiy3;
 445             faction[ii3+8]   = faction[ii3+8] + fiz3;
 446             faction[ii3+9]   = faction[ii3+9] + fix4;
 447             faction[ii3+10]  = faction[ii3+10] + fiy4;
 448             faction[ii3+11]  = faction[ii3+11] + fiz4;
 449             fshift[is3]      = fshift[is3]+fix1+fix2+fix3+fix4;
 450             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3+fiy4;
 451             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3+fiz4;
 452             ggid             = gid[n];
 453             Vc[ggid]         = Vc[ggid] + vctot;
 454             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 455             ninner           = ninner + nj1 - nj0;
 456         }
 457
 458         nouter           = nouter + nn1 - nn0;
 459     }
 460     while (nn1<nri);
 461
 462     *outeriter       = nouter;
 463     *inneriter       = ninner;
 464 }
 465
 466
 467
 468
 469
 470 /*
 471  * Gromacs nonbonded kernel nb_kernel134_adress_ex
 472  * Coulomb interaction:     Normal Coulomb
 473  * VdW interaction:         Tabulated
 474  * water optimization:      pairs of TIP4P interactions
 475  * Calculate forces:        yes
 476  */
 477 void nb_kernel134_adress_ex(
 478                     int *           p_nri,
 479                     int *           iinr,
 480                     int *           jindex,
 481                     int *           jjnr,
 482                     int *           shift,
 483                     real *         shiftvec,
 484                     real *         fshift,
 485                     int *           gid,
 486                     real *         pos,
 487                     real *         faction,
 488                     real *         charge,
 489                     real *         p_facel,
 490                     real *         p_krf,
 491                     real *         p_crf,
 492                     real *         Vc,
 493                     int *           type,
 494                     int *           p_ntype,
 495                     real *         vdwparam,
 496                     real *         Vvdw,
 497                     real *         p_tabscale,
 498                     real *         VFtab,
 499                     real *         invsqrta,
 500                     real *         dvda,
 501                     real *         p_gbtabscale,
 502                     real *         GBtab,
 503                     int *           p_nthreads,
 504                     int *           count,
 505                     void *          mtx,
 506                     int *           outeriter,
 507                     int *           inneriter,
 508                     real           force_cap,
 509                     real *         wf)
 510 {
 511     int           nri,ntype,nthreads;
 512     real         facel,krf,crf,tabscale,gbtabscale;
 513     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 514     int           nn0,nn1,nouter,ninner;
 515     real         shX,shY,shZ;
 516     real         fscal,tx,ty,tz;
 517     real         rinvsq;
 518     real         qq,vcoul,vctot;
 519     int           tj;
 520     real         Vvdw6,Vvdwtot;
 521     real         Vvdw12;
 522     real         r,rt,eps,eps2;
 523     int           n0,nnn;
 524     real         Y,F,Geps,Heps2,Fp,VV;
 525     real         FF;
 526     real         fijD,fijR;
 527     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 528     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
 529     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
 530     real         ix4,iy4,iz4,fix4,fiy4,fiz4;
 531     real         jx1,jy1,jz1;
 532     real         jx2,jy2,jz2,fjx2,fjy2,fjz2;
 533     real         jx3,jy3,jz3,fjx3,fjy3,fjz3;
 534     real         jx4,jy4,jz4,fjx4,fjy4,fjz4;
 535     real         dx11,dy11,dz11,rsq11,rinv11;
 536     real         dx22,dy22,dz22,rsq22,rinv22;
 537     real         dx23,dy23,dz23,rsq23,rinv23;
 538     real         dx24,dy24,dz24,rsq24,rinv24;
 539     real         dx32,dy32,dz32,rsq32,rinv32;
 540     real         dx33,dy33,dz33,rsq33,rinv33;
 541     real         dx34,dy34,dz34,rsq34,rinv34;
 542     real         dx42,dy42,dz42,rsq42,rinv42;
 543     real         dx43,dy43,dz43,rsq43,rinv43;
 544     real         dx44,dy44,dz44,rsq44,rinv44;
 545     real         qH,qM,qqMM,qqMH,qqHH;
 546     real         c6,c12;
 547     real         weight_cg1, weight_cg2, weight_product;
 548     real         hybscal;
 549
 550     nri              = *p_nri;
 551     ntype            = *p_ntype;
 552     nthreads         = *p_nthreads;
 553     facel            = *p_facel;
 554     krf              = *p_krf;
 555     crf              = *p_crf;
 556     tabscale         = *p_tabscale;
 557     ii               = iinr[0];
 558     qH               = charge[ii+1];
 559     qM               = charge[ii+3];
 560     qqMM             = facel*qM*qM;
 561     qqMH             = facel*qM*qH;
 562     qqHH             = facel*qH*qH;
 563     tj               = 2*(ntype+1)*type[ii];
 564     c6               = vdwparam[tj];
 565     c12              = vdwparam[tj+1];
 566
 567     nouter           = 0;
 568     ninner           = 0;
 569
 570     do
 571     {
 572         #ifdef GMX_THREAD_SHM_FDECOMP
 573         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 574         nn0              = *count;
 575         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 576         *count           = nn1;
 577         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 578         if(nn1>nri) nn1=nri;
 579         #else
 580         nn0 = 0;
 581         nn1 = nri;
 582         #endif
 583
 584         for(n=nn0; (n<nn1); n++)
 585         {
 586             is3              = 3*shift[n];
 587             shX              = shiftvec[is3];
 588             shY              = shiftvec[is3+1];
 589             shZ              = shiftvec[is3+2];
 590             nj0              = jindex[n];
 591             nj1              = jindex[n+1];
 592             ii               = iinr[n];
 593             ii3              = 3*ii;
 594             ix1              = shX + pos[ii3+0];
 595             iy1              = shY + pos[ii3+1];
 596             iz1              = shZ + pos[ii3+2];
 597             ix2              = shX + pos[ii3+3];
 598             iy2              = shY + pos[ii3+4];
 599             iz2              = shZ + pos[ii3+5];
 600             ix3              = shX + pos[ii3+6];
 601             iy3              = shY + pos[ii3+7];
 602             iz3              = shZ + pos[ii3+8];
 603             ix4              = shX + pos[ii3+9];
 604             iy4              = shY + pos[ii3+10];
 605             iz4              = shZ + pos[ii3+11];
 606             weight_cg1       = wf[ii];
 607             vctot            = 0;
 608             Vvdwtot          = 0;
 609             fix1             = 0;
 610             fiy1             = 0;
 611             fiz1             = 0;
 612             fix2             = 0;
 613             fiy2             = 0;
 614             fiz2             = 0;
 615             fix3             = 0;
 616             fiy3             = 0;
 617             fiz3             = 0;
 618             fix4             = 0;
 619             fiy4             = 0;
 620             fiz4             = 0;
 621
 622             for(k=nj0; (k<nj1); k++)
 623             {
 624                 jnr              = jjnr[k];
 625                 weight_cg2       = wf[jnr];
 626                 weight_product   = weight_cg1*weight_cg2;
 627                 if (weight_product < ALMOST_ZERO) {
 628                 /* force is zero, skip this molecule */
 629                  continue;
 630                 }
 631                 else if (weight_product >= ALMOST_ONE)
 632                 {
 633                        hybscal = 1.0;
 634                 }
 635                 else
 636                 {
 637                    hybscal = weight_product;
 638                 }
 639                 j3               = 3*jnr;
 640                 jx1              = pos[j3+0];
 641                 jy1              = pos[j3+1];
 642                 jz1              = pos[j3+2];
 643                 jx2              = pos[j3+3];
 644                 jy2              = pos[j3+4];
 645                 jz2              = pos[j3+5];
 646                 jx3              = pos[j3+6];
 647                 jy3              = pos[j3+7];
 648                 jz3              = pos[j3+8];
 649                 jx4              = pos[j3+9];
 650                 jy4              = pos[j3+10];
 651                 jz4              = pos[j3+11];
 652                 dx11             = ix1 - jx1;
 653                 dy11             = iy1 - jy1;
 654                 dz11             = iz1 - jz1;
 655                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 656                 dx22             = ix2 - jx2;
 657                 dy22             = iy2 - jy2;
 658                 dz22             = iz2 - jz2;
 659                 rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 660                 dx23             = ix2 - jx3;
 661                 dy23             = iy2 - jy3;
 662                 dz23             = iz2 - jz3;
 663                 rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 664                 dx24             = ix2 - jx4;
 665                 dy24             = iy2 - jy4;
 666                 dz24             = iz2 - jz4;
 667                 rsq24            = dx24*dx24+dy24*dy24+dz24*dz24;
 668                 dx32             = ix3 - jx2;
 669                 dy32             = iy3 - jy2;
 670                 dz32             = iz3 - jz2;
 671                 rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 672                 dx33             = ix3 - jx3;
 673                 dy33             = iy3 - jy3;
 674                 dz33             = iz3 - jz3;
 675                 rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 676                 dx34             = ix3 - jx4;
 677                 dy34             = iy3 - jy4;
 678                 dz34             = iz3 - jz4;
 679                 rsq34            = dx34*dx34+dy34*dy34+dz34*dz34;
 680                 dx42             = ix4 - jx2;
 681                 dy42             = iy4 - jy2;
 682                 dz42             = iz4 - jz2;
 683                 rsq42            = dx42*dx42+dy42*dy42+dz42*dz42;
 684                 dx43             = ix4 - jx3;
 685                 dy43             = iy4 - jy3;
 686                 dz43             = iz4 - jz3;
 687                 rsq43            = dx43*dx43+dy43*dy43+dz43*dz43;
 688                 dx44             = ix4 - jx4;
 689                 dy44             = iy4 - jy4;
 690                 dz44             = iz4 - jz4;
 691                 rsq44            = dx44*dx44+dy44*dy44+dz44*dz44;
 692                 rinv11           = 1.0/sqrt(rsq11);
 693                 rinv22           = 1.0/sqrt(rsq22);
 694                 rinv23           = 1.0/sqrt(rsq23);
 695                 rinv24           = 1.0/sqrt(rsq24);
 696                 rinv32           = 1.0/sqrt(rsq32);
 697                 rinv33           = 1.0/sqrt(rsq33);
 698                 rinv34           = 1.0/sqrt(rsq34);
 699                 rinv42           = 1.0/sqrt(rsq42);
 700                 rinv43           = 1.0/sqrt(rsq43);
 701                 rinv44           = 1.0/sqrt(rsq44);
 702                 r                = rsq11*rinv11;
 703                 rt               = r*tabscale;
 704                 n0               = rt;
 705                 eps              = rt-n0;
 706                 eps2             = eps*eps;
 707                 nnn              = 8*n0;
 708                 Y                = VFtab[nnn];
 709                 F                = VFtab[nnn+1];
 710                 Geps             = eps*VFtab[nnn+2];
 711                 Heps2            = eps2*VFtab[nnn+3];
 712                 Fp               = F+Geps+Heps2;
 713                 VV               = Y+eps*Fp;
 714                 FF               = Fp+Geps+2.0*Heps2;
 715                 Vvdw6            = c6*VV;
 716                 fijD             = c6*FF;
 717                 nnn              = nnn+4;
 718                 Y                = VFtab[nnn];
 719                 F                = VFtab[nnn+1];
 720                 Geps             = eps*VFtab[nnn+2];
 721                 Heps2            = eps2*VFtab[nnn+3];
 722                 Fp               = F+Geps+Heps2;
 723                 VV               = Y+eps*Fp;
 724                 FF               = Fp+Geps+2.0*Heps2;
 725                 Vvdw12           = c12*VV;
 726                 fijR             = c12*FF;
 727                 Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12;
 728                 fscal            = -((fijD+fijR)*tabscale)*rinv11;
 729                 fscal *= hybscal;
 730                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 731                 fscal=force_cap*fscal/fabs(fscal);
 732                 }
 733                 tx               = fscal*dx11;
 734                 ty               = fscal*dy11;
 735                 tz               = fscal*dz11;
 736                 fix1             = fix1 + tx;
 737                 fiy1             = fiy1 + ty;
 738                 fiz1             = fiz1 + tz;
 739                 faction[j3+0]    = faction[j3+0] - tx;
 740                 faction[j3+1]    = faction[j3+1] - ty;
 741                 faction[j3+2]    = faction[j3+2] - tz;
 742                 qq               = qqHH;
 743                 rinvsq           = rinv22*rinv22;
 744                 vcoul            = qq*rinv22;
 745                 vctot            = vctot+vcoul;
 746                 fscal            = (vcoul)*rinvsq;
 747                 fscal *= hybscal;
 748                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 749                 fscal=force_cap*fscal/fabs(fscal);
 750                 }
 751                 tx               = fscal*dx22;
 752                 ty               = fscal*dy22;
 753                 tz               = fscal*dz22;
 754                 fix2             = fix2 + tx;
 755                 fiy2             = fiy2 + ty;
 756                 fiz2             = fiz2 + tz;
 757                 fjx2             = faction[j3+3] - tx;
 758                 fjy2             = faction[j3+4] - ty;
 759                 fjz2             = faction[j3+5] - tz;
 760                 qq               = qqHH;
 761                 rinvsq           = rinv23*rinv23;
 762                 vcoul            = qq*rinv23;
 763                 vctot            = vctot+vcoul;
 764                 fscal            = (vcoul)*rinvsq;
 765                 fscal *= hybscal;
 766                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 767                 fscal=force_cap*fscal/fabs(fscal);
 768                 }
 769                 tx               = fscal*dx23;
 770                 ty               = fscal*dy23;
 771                 tz               = fscal*dz23;
 772                 fix2             = fix2 + tx;
 773                 fiy2             = fiy2 + ty;
 774                 fiz2             = fiz2 + tz;
 775                 fjx3             = faction[j3+6] - tx;
 776                 fjy3             = faction[j3+7] - ty;
 777                 fjz3             = faction[j3+8] - tz;
 778                 qq               = qqMH;
 779                 rinvsq           = rinv24*rinv24;
 780                 vcoul            = qq*rinv24;
 781                 vctot            = vctot+vcoul;
 782                 fscal            = (vcoul)*rinvsq;
 783                 fscal *= hybscal;
 784                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 785                 fscal=force_cap*fscal/fabs(fscal);
 786                 }
 787                 tx               = fscal*dx24;
 788                 ty               = fscal*dy24;
 789                 tz               = fscal*dz24;
 790                 fix2             = fix2 + tx;
 791                 fiy2             = fiy2 + ty;
 792                 fiz2             = fiz2 + tz;
 793                 fjx4             = faction[j3+9] - tx;
 794                 fjy4             = faction[j3+10] - ty;
 795                 fjz4             = faction[j3+11] - tz;
 796                 qq               = qqHH;
 797                 rinvsq           = rinv32*rinv32;
 798                 vcoul            = qq*rinv32;
 799                 vctot            = vctot+vcoul;
 800                 fscal            = (vcoul)*rinvsq;
 801                 fscal *= hybscal;
 802                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 803                 fscal=force_cap*fscal/fabs(fscal);
 804                 }
 805                 tx               = fscal*dx32;
 806                 ty               = fscal*dy32;
 807                 tz               = fscal*dz32;
 808                 fix3             = fix3 + tx;
 809                 fiy3             = fiy3 + ty;
 810                 fiz3             = fiz3 + tz;
 811                 fjx2             = fjx2 - tx;
 812                 fjy2             = fjy2 - ty;
 813                 fjz2             = fjz2 - tz;
 814                 qq               = qqHH;
 815                 rinvsq           = rinv33*rinv33;
 816                 vcoul            = qq*rinv33;
 817                 vctot            = vctot+vcoul;
 818                 fscal            = (vcoul)*rinvsq;
 819                 fscal *= hybscal;
 820                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 821                 fscal=force_cap*fscal/fabs(fscal);
 822                 }
 823                 tx               = fscal*dx33;
 824                 ty               = fscal*dy33;
 825                 tz               = fscal*dz33;
 826                 fix3             = fix3 + tx;
 827                 fiy3             = fiy3 + ty;
 828                 fiz3             = fiz3 + tz;
 829                 fjx3             = fjx3 - tx;
 830                 fjy3             = fjy3 - ty;
 831                 fjz3             = fjz3 - tz;
 832                 qq               = qqMH;
 833                 rinvsq           = rinv34*rinv34;
 834                 vcoul            = qq*rinv34;
 835                 vctot            = vctot+vcoul;
 836                 fscal            = (vcoul)*rinvsq;
 837                 fscal *= hybscal;
 838                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 839                 fscal=force_cap*fscal/fabs(fscal);
 840                 }
 841                 tx               = fscal*dx34;
 842                 ty               = fscal*dy34;
 843                 tz               = fscal*dz34;
 844                 fix3             = fix3 + tx;
 845                 fiy3             = fiy3 + ty;
 846                 fiz3             = fiz3 + tz;
 847                 fjx4             = fjx4 - tx;
 848                 fjy4             = fjy4 - ty;
 849                 fjz4             = fjz4 - tz;
 850                 qq               = qqMH;
 851                 rinvsq           = rinv42*rinv42;
 852                 vcoul            = qq*rinv42;
 853                 vctot            = vctot+vcoul;
 854                 fscal            = (vcoul)*rinvsq;
 855                 fscal *= hybscal;
 856                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 857                 fscal=force_cap*fscal/fabs(fscal);
 858                 }
 859                 tx               = fscal*dx42;
 860                 ty               = fscal*dy42;
 861                 tz               = fscal*dz42;
 862                 fix4             = fix4 + tx;
 863                 fiy4             = fiy4 + ty;
 864                 fiz4             = fiz4 + tz;
 865                 faction[j3+3]    = fjx2 - tx;
 866                 faction[j3+4]    = fjy2 - ty;
 867                 faction[j3+5]    = fjz2 - tz;
 868                 qq               = qqMH;
 869                 rinvsq           = rinv43*rinv43;
 870                 vcoul            = qq*rinv43;
 871                 vctot            = vctot+vcoul;
 872                 fscal            = (vcoul)*rinvsq;
 873                 fscal *= hybscal;
 874                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 875                 fscal=force_cap*fscal/fabs(fscal);
 876                 }
 877                 tx               = fscal*dx43;
 878                 ty               = fscal*dy43;
 879                 tz               = fscal*dz43;
 880                 fix4             = fix4 + tx;
 881                 fiy4             = fiy4 + ty;
 882                 fiz4             = fiz4 + tz;
 883                 faction[j3+6]    = fjx3 - tx;
 884                 faction[j3+7]    = fjy3 - ty;
 885                 faction[j3+8]    = fjz3 - tz;
 886                 qq               = qqMM;
 887                 rinvsq           = rinv44*rinv44;
 888                 vcoul            = qq*rinv44;
 889                 vctot            = vctot+vcoul;
 890                 fscal            = (vcoul)*rinvsq;
 891                 fscal *= hybscal;
 892                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 893                 fscal=force_cap*fscal/fabs(fscal);
 894                 }
 895                 tx               = fscal*dx44;
 896                 ty               = fscal*dy44;
 897                 tz               = fscal*dz44;
 898                 fix4             = fix4 + tx;
 899                 fiy4             = fiy4 + ty;
 900                 fiz4             = fiz4 + tz;
 901                 faction[j3+9]    = fjx4 - tx;
 902                 faction[j3+10]   = fjy4 - ty;
 903                 faction[j3+11]   = fjz4 - tz;
 904             }
 905
 906             faction[ii3+0]   = faction[ii3+0] + fix1;
 907             faction[ii3+1]   = faction[ii3+1] + fiy1;
 908             faction[ii3+2]   = faction[ii3+2] + fiz1;
 909             faction[ii3+3]   = faction[ii3+3] + fix2;
 910             faction[ii3+4]   = faction[ii3+4] + fiy2;
 911             faction[ii3+5]   = faction[ii3+5] + fiz2;
 912             faction[ii3+6]   = faction[ii3+6] + fix3;
 913             faction[ii3+7]   = faction[ii3+7] + fiy3;
 914             faction[ii3+8]   = faction[ii3+8] + fiz3;
 915             faction[ii3+9]   = faction[ii3+9] + fix4;
 916             faction[ii3+10]  = faction[ii3+10] + fiy4;
 917             faction[ii3+11]  = faction[ii3+11] + fiz4;
 918             fshift[is3]      = fshift[is3]+fix1+fix2+fix3+fix4;
 919             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3+fiy4;
 920             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3+fiz4;
 921             ggid             = gid[n];
 922             Vc[ggid]         = Vc[ggid] + vctot;
 923             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 924             ninner           = ninner + nj1 - nj0;
 925         }
 926
 927         nouter           = nouter + nn1 - nn0;
 928     }
 929     while (nn1<nri);
 930
 931     *outeriter       = nouter;
 932     *inneriter       = ninner;
 933 }
 934
 935