src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel302_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_kernel302_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel302_adress_cg
  33  * Coulomb interaction:     Tabulated
  34  * VdW interaction:         Not calculated
  35  * water optimization:      pairs of SPC/TIP3P interactions
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel302_adress_cg(
  39                     int *           p_nri,
  40                     int *           iinr,
  41                     int *           jindex,
  42                     int *           jjnr,
  43                     int *           shift,
  44                     real *         shiftvec,
  45                     real *         fshift,
  46                     int *           gid,
  47                     real *         pos,
  48                     real *         faction,
  49                     real *         charge,
  50                     real *         p_facel,
  51                     real *         p_krf,
  52                     real *         p_crf,
  53                     real *         Vc,
  54                     int *           type,
  55                     int *           p_ntype,
  56                     real *         vdwparam,
  57                     real *         Vvdw,
  58                     real *         p_tabscale,
  59                     real *         VFtab,
  60                     real *         invsqrta,
  61                     real *         dvda,
  62                     real *         p_gbtabscale,
  63                     real *         GBtab,
  64                     int *           p_nthreads,
  65                     int *           count,
  66                     void *          mtx,
  67                     int *           outeriter,
  68                     int *           inneriter,
  69                     real           force_cap,
  70                     real *         wf)
  71 {
  72     int           nri,ntype,nthreads;
  73     real         facel,krf,crf,tabscale,gbtabscale;
  74     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
  75     int           nn0,nn1,nouter,ninner;
  76     real         shX,shY,shZ;
  77     real         fscal,tx,ty,tz;
  78     real         qq,vcoul,vctot;
  79     real         r,rt,eps,eps2;
  80     int           n0,nnn;
  81     real         Y,F,Geps,Heps2,Fp,VV;
  82     real         FF;
  83     real         fijC;
  84     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
  85     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
  86     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
  87     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
  88     real         jx2,jy2,jz2,fjx2,fjy2,fjz2;
  89     real         jx3,jy3,jz3,fjx3,fjy3,fjz3;
  90     real         dx11,dy11,dz11,rsq11,rinv11;
  91     real         dx12,dy12,dz12,rsq12,rinv12;
  92     real         dx13,dy13,dz13,rsq13,rinv13;
  93     real         dx21,dy21,dz21,rsq21,rinv21;
  94     real         dx22,dy22,dz22,rsq22,rinv22;
  95     real         dx23,dy23,dz23,rsq23,rinv23;
  96     real         dx31,dy31,dz31,rsq31,rinv31;
  97     real         dx32,dy32,dz32,rsq32,rinv32;
  98     real         dx33,dy33,dz33,rsq33,rinv33;
  99     real         qO,qH,qqOO,qqOH,qqHH;
 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               = charge[ii];
 112     qH               = charge[ii+1];
 113     qqOO             = facel*qO*qO;
 114     qqOH             = facel*qO*qH;
 115     qqHH             = facel*qH*qH;
 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             weight_cg1       = wf[ii];
 154             vctot            = 0;
 155             fix1             = 0;
 156             fiy1             = 0;
 157             fiz1             = 0;
 158             fix2             = 0;
 159             fiy2             = 0;
 160             fiz2             = 0;
 161             fix3             = 0;
 162             fiy3             = 0;
 163             fiz3             = 0;
 164
 165             for(k=nj0; (k<nj1); k++)
 166             {
 167                 jnr              = jjnr[k];
 168                 weight_cg2       = wf[jnr];
 169                 weight_product   = weight_cg1*weight_cg2;
 170                 if (weight_product < ALMOST_ZERO) {
 171                        hybscal = 1.0;
 172                 }
 173                 else if (weight_product >= ALMOST_ONE)
 174                 {
 175                   /* force is zero, skip this molecule */
 176                        continue;
 177                 }
 178                 else
 179                 {
 180                    hybscal = 1.0 - weight_product;
 181                 }
 182                 j3               = 3*jnr;
 183                 jx1              = pos[j3+0];
 184                 jy1              = pos[j3+1];
 185                 jz1              = pos[j3+2];
 186                 jx2              = pos[j3+3];
 187                 jy2              = pos[j3+4];
 188                 jz2              = pos[j3+5];
 189                 jx3              = pos[j3+6];
 190                 jy3              = pos[j3+7];
 191                 jz3              = pos[j3+8];
 192                 dx11             = ix1 - jx1;
 193                 dy11             = iy1 - jy1;
 194                 dz11             = iz1 - jz1;
 195                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 196                 dx12             = ix1 - jx2;
 197                 dy12             = iy1 - jy2;
 198                 dz12             = iz1 - jz2;
 199                 rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 200                 dx13             = ix1 - jx3;
 201                 dy13             = iy1 - jy3;
 202                 dz13             = iz1 - jz3;
 203                 rsq13            = dx13*dx13+dy13*dy13+dz13*dz13;
 204                 dx21             = ix2 - jx1;
 205                 dy21             = iy2 - jy1;
 206                 dz21             = iz2 - jz1;
 207                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 208                 dx22             = ix2 - jx2;
 209                 dy22             = iy2 - jy2;
 210                 dz22             = iz2 - jz2;
 211                 rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 212                 dx23             = ix2 - jx3;
 213                 dy23             = iy2 - jy3;
 214                 dz23             = iz2 - jz3;
 215                 rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 216                 dx31             = ix3 - jx1;
 217                 dy31             = iy3 - jy1;
 218                 dz31             = iz3 - jz1;
 219                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 220                 dx32             = ix3 - jx2;
 221                 dy32             = iy3 - jy2;
 222                 dz32             = iz3 - jz2;
 223                 rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 224                 dx33             = ix3 - jx3;
 225                 dy33             = iy3 - jy3;
 226                 dz33             = iz3 - jz3;
 227                 rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 228                 rinv11           = 1.0/sqrt(rsq11);
 229                 rinv12           = 1.0/sqrt(rsq12);
 230                 rinv13           = 1.0/sqrt(rsq13);
 231                 rinv21           = 1.0/sqrt(rsq21);
 232                 rinv22           = 1.0/sqrt(rsq22);
 233                 rinv23           = 1.0/sqrt(rsq23);
 234                 rinv31           = 1.0/sqrt(rsq31);
 235                 rinv32           = 1.0/sqrt(rsq32);
 236                 rinv33           = 1.0/sqrt(rsq33);
 237                 qq               = qqOO;
 238                 r                = rsq11*rinv11;
 239                 rt               = r*tabscale;
 240                 n0               = rt;
 241                 eps              = rt-n0;
 242                 eps2             = eps*eps;
 243                 nnn              = 4*n0;
 244                 Y                = VFtab[nnn];
 245                 F                = VFtab[nnn+1];
 246                 Geps             = eps*VFtab[nnn+2];
 247                 Heps2            = eps2*VFtab[nnn+3];
 248                 Fp               = F+Geps+Heps2;
 249                 VV               = Y+eps*Fp;
 250                 FF               = Fp+Geps+2.0*Heps2;
 251                 vcoul            = qq*VV;
 252                 fijC             = qq*FF;
 253                 vctot            = vctot + vcoul;
 254                 fscal            = -((fijC)*tabscale)*rinv11;
 255                 fscal *= hybscal;
 256                 tx               = fscal*dx11;
 257                 ty               = fscal*dy11;
 258                 tz               = fscal*dz11;
 259                 fix1             = fix1 + tx;
 260                 fiy1             = fiy1 + ty;
 261                 fiz1             = fiz1 + tz;
 262                 fjx1             = faction[j3+0] - tx;
 263                 fjy1             = faction[j3+1] - ty;
 264                 fjz1             = faction[j3+2] - tz;
 265                 qq               = qqOH;
 266                 r                = rsq12*rinv12;
 267                 rt               = r*tabscale;
 268                 n0               = rt;
 269                 eps              = rt-n0;
 270                 eps2             = eps*eps;
 271                 nnn              = 4*n0;
 272                 Y                = VFtab[nnn];
 273                 F                = VFtab[nnn+1];
 274                 Geps             = eps*VFtab[nnn+2];
 275                 Heps2            = eps2*VFtab[nnn+3];
 276                 Fp               = F+Geps+Heps2;
 277                 VV               = Y+eps*Fp;
 278                 FF               = Fp+Geps+2.0*Heps2;
 279                 vcoul            = qq*VV;
 280                 fijC             = qq*FF;
 281                 vctot            = vctot + vcoul;
 282                 fscal            = -((fijC)*tabscale)*rinv12;
 283                 fscal *= hybscal;
 284                 tx               = fscal*dx12;
 285                 ty               = fscal*dy12;
 286                 tz               = fscal*dz12;
 287                 fix1             = fix1 + tx;
 288                 fiy1             = fiy1 + ty;
 289                 fiz1             = fiz1 + tz;
 290                 fjx2             = faction[j3+3] - tx;
 291                 fjy2             = faction[j3+4] - ty;
 292                 fjz2             = faction[j3+5] - tz;
 293                 qq               = qqOH;
 294                 r                = rsq13*rinv13;
 295                 rt               = r*tabscale;
 296                 n0               = rt;
 297                 eps              = rt-n0;
 298                 eps2             = eps*eps;
 299                 nnn              = 4*n0;
 300                 Y                = VFtab[nnn];
 301                 F                = VFtab[nnn+1];
 302                 Geps             = eps*VFtab[nnn+2];
 303                 Heps2            = eps2*VFtab[nnn+3];
 304                 Fp               = F+Geps+Heps2;
 305                 VV               = Y+eps*Fp;
 306                 FF               = Fp+Geps+2.0*Heps2;
 307                 vcoul            = qq*VV;
 308                 fijC             = qq*FF;
 309                 vctot            = vctot + vcoul;
 310                 fscal            = -((fijC)*tabscale)*rinv13;
 311                 fscal *= hybscal;
 312                 tx               = fscal*dx13;
 313                 ty               = fscal*dy13;
 314                 tz               = fscal*dz13;
 315                 fix1             = fix1 + tx;
 316                 fiy1             = fiy1 + ty;
 317                 fiz1             = fiz1 + tz;
 318                 fjx3             = faction[j3+6] - tx;
 319                 fjy3             = faction[j3+7] - ty;
 320                 fjz3             = faction[j3+8] - tz;
 321                 qq               = qqOH;
 322                 r                = rsq21*rinv21;
 323                 rt               = r*tabscale;
 324                 n0               = rt;
 325                 eps              = rt-n0;
 326                 eps2             = eps*eps;
 327                 nnn              = 4*n0;
 328                 Y                = VFtab[nnn];
 329                 F                = VFtab[nnn+1];
 330                 Geps             = eps*VFtab[nnn+2];
 331                 Heps2            = eps2*VFtab[nnn+3];
 332                 Fp               = F+Geps+Heps2;
 333                 VV               = Y+eps*Fp;
 334                 FF               = Fp+Geps+2.0*Heps2;
 335                 vcoul            = qq*VV;
 336                 fijC             = qq*FF;
 337                 vctot            = vctot + vcoul;
 338                 fscal            = -((fijC)*tabscale)*rinv21;
 339                 fscal *= hybscal;
 340                 tx               = fscal*dx21;
 341                 ty               = fscal*dy21;
 342                 tz               = fscal*dz21;
 343                 fix2             = fix2 + tx;
 344                 fiy2             = fiy2 + ty;
 345                 fiz2             = fiz2 + tz;
 346                 fjx1             = fjx1 - tx;
 347                 fjy1             = fjy1 - ty;
 348                 fjz1             = fjz1 - tz;
 349                 qq               = qqHH;
 350                 r                = rsq22*rinv22;
 351                 rt               = r*tabscale;
 352                 n0               = rt;
 353                 eps              = rt-n0;
 354                 eps2             = eps*eps;
 355                 nnn              = 4*n0;
 356                 Y                = VFtab[nnn];
 357                 F                = VFtab[nnn+1];
 358                 Geps             = eps*VFtab[nnn+2];
 359                 Heps2            = eps2*VFtab[nnn+3];
 360                 Fp               = F+Geps+Heps2;
 361                 VV               = Y+eps*Fp;
 362                 FF               = Fp+Geps+2.0*Heps2;
 363                 vcoul            = qq*VV;
 364                 fijC             = qq*FF;
 365                 vctot            = vctot + vcoul;
 366                 fscal            = -((fijC)*tabscale)*rinv22;
 367                 fscal *= hybscal;
 368                 tx               = fscal*dx22;
 369                 ty               = fscal*dy22;
 370                 tz               = fscal*dz22;
 371                 fix2             = fix2 + tx;
 372                 fiy2             = fiy2 + ty;
 373                 fiz2             = fiz2 + tz;
 374                 fjx2             = fjx2 - tx;
 375                 fjy2             = fjy2 - ty;
 376                 fjz2             = fjz2 - tz;
 377                 qq               = qqHH;
 378                 r                = rsq23*rinv23;
 379                 rt               = r*tabscale;
 380                 n0               = rt;
 381                 eps              = rt-n0;
 382                 eps2             = eps*eps;
 383                 nnn              = 4*n0;
 384                 Y                = VFtab[nnn];
 385                 F                = VFtab[nnn+1];
 386                 Geps             = eps*VFtab[nnn+2];
 387                 Heps2            = eps2*VFtab[nnn+3];
 388                 Fp               = F+Geps+Heps2;
 389                 VV               = Y+eps*Fp;
 390                 FF               = Fp+Geps+2.0*Heps2;
 391                 vcoul            = qq*VV;
 392                 fijC             = qq*FF;
 393                 vctot            = vctot + vcoul;
 394                 fscal            = -((fijC)*tabscale)*rinv23;
 395                 fscal *= hybscal;
 396                 tx               = fscal*dx23;
 397                 ty               = fscal*dy23;
 398                 tz               = fscal*dz23;
 399                 fix2             = fix2 + tx;
 400                 fiy2             = fiy2 + ty;
 401                 fiz2             = fiz2 + tz;
 402                 fjx3             = fjx3 - tx;
 403                 fjy3             = fjy3 - ty;
 404                 fjz3             = fjz3 - tz;
 405                 qq               = qqOH;
 406                 r                = rsq31*rinv31;
 407                 rt               = r*tabscale;
 408                 n0               = rt;
 409                 eps              = rt-n0;
 410                 eps2             = eps*eps;
 411                 nnn              = 4*n0;
 412                 Y                = VFtab[nnn];
 413                 F                = VFtab[nnn+1];
 414                 Geps             = eps*VFtab[nnn+2];
 415                 Heps2            = eps2*VFtab[nnn+3];
 416                 Fp               = F+Geps+Heps2;
 417                 VV               = Y+eps*Fp;
 418                 FF               = Fp+Geps+2.0*Heps2;
 419                 vcoul            = qq*VV;
 420                 fijC             = qq*FF;
 421                 vctot            = vctot + vcoul;
 422                 fscal            = -((fijC)*tabscale)*rinv31;
 423                 fscal *= hybscal;
 424                 tx               = fscal*dx31;
 425                 ty               = fscal*dy31;
 426                 tz               = fscal*dz31;
 427                 fix3             = fix3 + tx;
 428                 fiy3             = fiy3 + ty;
 429                 fiz3             = fiz3 + tz;
 430                 faction[j3+0]    = fjx1 - tx;
 431                 faction[j3+1]    = fjy1 - ty;
 432                 faction[j3+2]    = fjz1 - tz;
 433                 qq               = qqHH;
 434                 r                = rsq32*rinv32;
 435                 rt               = r*tabscale;
 436                 n0               = rt;
 437                 eps              = rt-n0;
 438                 eps2             = eps*eps;
 439                 nnn              = 4*n0;
 440                 Y                = VFtab[nnn];
 441                 F                = VFtab[nnn+1];
 442                 Geps             = eps*VFtab[nnn+2];
 443                 Heps2            = eps2*VFtab[nnn+3];
 444                 Fp               = F+Geps+Heps2;
 445                 VV               = Y+eps*Fp;
 446                 FF               = Fp+Geps+2.0*Heps2;
 447                 vcoul            = qq*VV;
 448                 fijC             = qq*FF;
 449                 vctot            = vctot + vcoul;
 450                 fscal            = -((fijC)*tabscale)*rinv32;
 451                 fscal *= hybscal;
 452                 tx               = fscal*dx32;
 453                 ty               = fscal*dy32;
 454                 tz               = fscal*dz32;
 455                 fix3             = fix3 + tx;
 456                 fiy3             = fiy3 + ty;
 457                 fiz3             = fiz3 + tz;
 458                 faction[j3+3]    = fjx2 - tx;
 459                 faction[j3+4]    = fjy2 - ty;
 460                 faction[j3+5]    = fjz2 - tz;
 461                 qq               = qqHH;
 462                 r                = rsq33*rinv33;
 463                 rt               = r*tabscale;
 464                 n0               = rt;
 465                 eps              = rt-n0;
 466                 eps2             = eps*eps;
 467                 nnn              = 4*n0;
 468                 Y                = VFtab[nnn];
 469                 F                = VFtab[nnn+1];
 470                 Geps             = eps*VFtab[nnn+2];
 471                 Heps2            = eps2*VFtab[nnn+3];
 472                 Fp               = F+Geps+Heps2;
 473                 VV               = Y+eps*Fp;
 474                 FF               = Fp+Geps+2.0*Heps2;
 475                 vcoul            = qq*VV;
 476                 fijC             = qq*FF;
 477                 vctot            = vctot + vcoul;
 478                 fscal            = -((fijC)*tabscale)*rinv33;
 479                 fscal *= hybscal;
 480                 tx               = fscal*dx33;
 481                 ty               = fscal*dy33;
 482                 tz               = fscal*dz33;
 483                 fix3             = fix3 + tx;
 484                 fiy3             = fiy3 + ty;
 485                 fiz3             = fiz3 + tz;
 486                 faction[j3+6]    = fjx3 - tx;
 487                 faction[j3+7]    = fjy3 - ty;
 488                 faction[j3+8]    = fjz3 - tz;
 489             }
 490
 491             faction[ii3+0]   = faction[ii3+0] + fix1;
 492             faction[ii3+1]   = faction[ii3+1] + fiy1;
 493             faction[ii3+2]   = faction[ii3+2] + fiz1;
 494             faction[ii3+3]   = faction[ii3+3] + fix2;
 495             faction[ii3+4]   = faction[ii3+4] + fiy2;
 496             faction[ii3+5]   = faction[ii3+5] + fiz2;
 497             faction[ii3+6]   = faction[ii3+6] + fix3;
 498             faction[ii3+7]   = faction[ii3+7] + fiy3;
 499             faction[ii3+8]   = faction[ii3+8] + fiz3;
 500             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
 501             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
 502             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
 503             ggid             = gid[n];
 504             Vc[ggid]         = Vc[ggid] + vctot;
 505             ninner           = ninner + nj1 - nj0;
 506         }
 507
 508         nouter           = nouter + nn1 - nn0;
 509     }
 510     while (nn1<nri);
 511
 512     *outeriter       = nouter;
 513     *inneriter       = ninner;
 514 }
 515
 516
 517
 518
 519
 520 /*
 521  * Gromacs nonbonded kernel nb_kernel302_adress_ex
 522  * Coulomb interaction:     Tabulated
 523  * VdW interaction:         Not calculated
 524  * water optimization:      pairs of SPC/TIP3P interactions
 525  * Calculate forces:        yes
 526  */
 527 void nb_kernel302_adress_ex(
 528                     int *           p_nri,
 529                     int *           iinr,
 530                     int *           jindex,
 531                     int *           jjnr,
 532                     int *           shift,
 533                     real *         shiftvec,
 534                     real *         fshift,
 535                     int *           gid,
 536                     real *         pos,
 537                     real *         faction,
 538                     real *         charge,
 539                     real *         p_facel,
 540                     real *         p_krf,
 541                     real *         p_crf,
 542                     real *         Vc,
 543                     int *           type,
 544                     int *           p_ntype,
 545                     real *         vdwparam,
 546                     real *         Vvdw,
 547                     real *         p_tabscale,
 548                     real *         VFtab,
 549                     real *         invsqrta,
 550                     real *         dvda,
 551                     real *         p_gbtabscale,
 552                     real *         GBtab,
 553                     int *           p_nthreads,
 554                     int *           count,
 555                     void *          mtx,
 556                     int *           outeriter,
 557                     int *           inneriter,
 558                     real           force_cap,
 559                     real *         wf)
 560 {
 561     int           nri,ntype,nthreads;
 562     real         facel,krf,crf,tabscale,gbtabscale;
 563     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 564     int           nn0,nn1,nouter,ninner;
 565     real         shX,shY,shZ;
 566     real         fscal,tx,ty,tz;
 567     real         qq,vcoul,vctot;
 568     real         r,rt,eps,eps2;
 569     int           n0,nnn;
 570     real         Y,F,Geps,Heps2,Fp,VV;
 571     real         FF;
 572     real         fijC;
 573     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 574     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
 575     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
 576     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
 577     real         jx2,jy2,jz2,fjx2,fjy2,fjz2;
 578     real         jx3,jy3,jz3,fjx3,fjy3,fjz3;
 579     real         dx11,dy11,dz11,rsq11,rinv11;
 580     real         dx12,dy12,dz12,rsq12,rinv12;
 581     real         dx13,dy13,dz13,rsq13,rinv13;
 582     real         dx21,dy21,dz21,rsq21,rinv21;
 583     real         dx22,dy22,dz22,rsq22,rinv22;
 584     real         dx23,dy23,dz23,rsq23,rinv23;
 585     real         dx31,dy31,dz31,rsq31,rinv31;
 586     real         dx32,dy32,dz32,rsq32,rinv32;
 587     real         dx33,dy33,dz33,rsq33,rinv33;
 588     real         qO,qH,qqOO,qqOH,qqHH;
 589     real         weight_cg1, weight_cg2, weight_product;
 590     real         hybscal;
 591
 592     nri              = *p_nri;
 593     ntype            = *p_ntype;
 594     nthreads         = *p_nthreads;
 595     facel            = *p_facel;
 596     krf              = *p_krf;
 597     crf              = *p_crf;
 598     tabscale         = *p_tabscale;
 599     ii               = iinr[0];
 600     qO               = charge[ii];
 601     qH               = charge[ii+1];
 602     qqOO             = facel*qO*qO;
 603     qqOH             = facel*qO*qH;
 604     qqHH             = facel*qH*qH;
 605
 606     nouter           = 0;
 607     ninner           = 0;
 608
 609     do
 610     {
 611         #ifdef GMX_THREAD_SHM_FDECOMP
 612         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 613         nn0              = *count;
 614         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 615         *count           = nn1;
 616         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 617         if(nn1>nri) nn1=nri;
 618         #else
 619         nn0 = 0;
 620         nn1 = nri;
 621         #endif
 622
 623         for(n=nn0; (n<nn1); n++)
 624         {
 625             is3              = 3*shift[n];
 626             shX              = shiftvec[is3];
 627             shY              = shiftvec[is3+1];
 628             shZ              = shiftvec[is3+2];
 629             nj0              = jindex[n];
 630             nj1              = jindex[n+1];
 631             ii               = iinr[n];
 632             ii3              = 3*ii;
 633             ix1              = shX + pos[ii3+0];
 634             iy1              = shY + pos[ii3+1];
 635             iz1              = shZ + pos[ii3+2];
 636             ix2              = shX + pos[ii3+3];
 637             iy2              = shY + pos[ii3+4];
 638             iz2              = shZ + pos[ii3+5];
 639             ix3              = shX + pos[ii3+6];
 640             iy3              = shY + pos[ii3+7];
 641             iz3              = shZ + pos[ii3+8];
 642             weight_cg1       = wf[ii];
 643             vctot            = 0;
 644             fix1             = 0;
 645             fiy1             = 0;
 646             fiz1             = 0;
 647             fix2             = 0;
 648             fiy2             = 0;
 649             fiz2             = 0;
 650             fix3             = 0;
 651             fiy3             = 0;
 652             fiz3             = 0;
 653
 654             for(k=nj0; (k<nj1); k++)
 655             {
 656                 jnr              = jjnr[k];
 657                 weight_cg2       = wf[jnr];
 658                 weight_product   = weight_cg1*weight_cg2;
 659                 if (weight_product < ALMOST_ZERO) {
 660                 /* force is zero, skip this molecule */
 661                  continue;
 662                 }
 663                 else if (weight_product >= ALMOST_ONE)
 664                 {
 665                        hybscal = 1.0;
 666                 }
 667                 else
 668                 {
 669                    hybscal = weight_product;
 670                 }
 671                 j3               = 3*jnr;
 672                 jx1              = pos[j3+0];
 673                 jy1              = pos[j3+1];
 674                 jz1              = pos[j3+2];
 675                 jx2              = pos[j3+3];
 676                 jy2              = pos[j3+4];
 677                 jz2              = pos[j3+5];
 678                 jx3              = pos[j3+6];
 679                 jy3              = pos[j3+7];
 680                 jz3              = pos[j3+8];
 681                 dx11             = ix1 - jx1;
 682                 dy11             = iy1 - jy1;
 683                 dz11             = iz1 - jz1;
 684                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 685                 dx12             = ix1 - jx2;
 686                 dy12             = iy1 - jy2;
 687                 dz12             = iz1 - jz2;
 688                 rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 689                 dx13             = ix1 - jx3;
 690                 dy13             = iy1 - jy3;
 691                 dz13             = iz1 - jz3;
 692                 rsq13            = dx13*dx13+dy13*dy13+dz13*dz13;
 693                 dx21             = ix2 - jx1;
 694                 dy21             = iy2 - jy1;
 695                 dz21             = iz2 - jz1;
 696                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 697                 dx22             = ix2 - jx2;
 698                 dy22             = iy2 - jy2;
 699                 dz22             = iz2 - jz2;
 700                 rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 701                 dx23             = ix2 - jx3;
 702                 dy23             = iy2 - jy3;
 703                 dz23             = iz2 - jz3;
 704                 rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 705                 dx31             = ix3 - jx1;
 706                 dy31             = iy3 - jy1;
 707                 dz31             = iz3 - jz1;
 708                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 709                 dx32             = ix3 - jx2;
 710                 dy32             = iy3 - jy2;
 711                 dz32             = iz3 - jz2;
 712                 rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 713                 dx33             = ix3 - jx3;
 714                 dy33             = iy3 - jy3;
 715                 dz33             = iz3 - jz3;
 716                 rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 717                 rinv11           = 1.0/sqrt(rsq11);
 718                 rinv12           = 1.0/sqrt(rsq12);
 719                 rinv13           = 1.0/sqrt(rsq13);
 720                 rinv21           = 1.0/sqrt(rsq21);
 721                 rinv22           = 1.0/sqrt(rsq22);
 722                 rinv23           = 1.0/sqrt(rsq23);
 723                 rinv31           = 1.0/sqrt(rsq31);
 724                 rinv32           = 1.0/sqrt(rsq32);
 725                 rinv33           = 1.0/sqrt(rsq33);
 726                 qq               = qqOO;
 727                 r                = rsq11*rinv11;
 728                 rt               = r*tabscale;
 729                 n0               = rt;
 730                 eps              = rt-n0;
 731                 eps2             = eps*eps;
 732                 nnn              = 4*n0;
 733                 Y                = VFtab[nnn];
 734                 F                = VFtab[nnn+1];
 735                 Geps             = eps*VFtab[nnn+2];
 736                 Heps2            = eps2*VFtab[nnn+3];
 737                 Fp               = F+Geps+Heps2;
 738                 VV               = Y+eps*Fp;
 739                 FF               = Fp+Geps+2.0*Heps2;
 740                 vcoul            = qq*VV;
 741                 fijC             = qq*FF;
 742                 vctot            = vctot + vcoul;
 743                 fscal            = -((fijC)*tabscale)*rinv11;
 744                 fscal *= hybscal;
 745                 tx               = fscal*dx11;
 746                 ty               = fscal*dy11;
 747                 tz               = fscal*dz11;
 748                 fix1             = fix1 + tx;
 749                 fiy1             = fiy1 + ty;
 750                 fiz1             = fiz1 + tz;
 751                 fjx1             = faction[j3+0] - tx;
 752                 fjy1             = faction[j3+1] - ty;
 753                 fjz1             = faction[j3+2] - tz;
 754                 qq               = qqOH;
 755                 r                = rsq12*rinv12;
 756                 rt               = r*tabscale;
 757                 n0               = rt;
 758                 eps              = rt-n0;
 759                 eps2             = eps*eps;
 760                 nnn              = 4*n0;
 761                 Y                = VFtab[nnn];
 762                 F                = VFtab[nnn+1];
 763                 Geps             = eps*VFtab[nnn+2];
 764                 Heps2            = eps2*VFtab[nnn+3];
 765                 Fp               = F+Geps+Heps2;
 766                 VV               = Y+eps*Fp;
 767                 FF               = Fp+Geps+2.0*Heps2;
 768                 vcoul            = qq*VV;
 769                 fijC             = qq*FF;
 770                 vctot            = vctot + vcoul;
 771                 fscal            = -((fijC)*tabscale)*rinv12;
 772                 fscal *= hybscal;
 773                 tx               = fscal*dx12;
 774                 ty               = fscal*dy12;
 775                 tz               = fscal*dz12;
 776                 fix1             = fix1 + tx;
 777                 fiy1             = fiy1 + ty;
 778                 fiz1             = fiz1 + tz;
 779                 fjx2             = faction[j3+3] - tx;
 780                 fjy2             = faction[j3+4] - ty;
 781                 fjz2             = faction[j3+5] - tz;
 782                 qq               = qqOH;
 783                 r                = rsq13*rinv13;
 784                 rt               = r*tabscale;
 785                 n0               = rt;
 786                 eps              = rt-n0;
 787                 eps2             = eps*eps;
 788                 nnn              = 4*n0;
 789                 Y                = VFtab[nnn];
 790                 F                = VFtab[nnn+1];
 791                 Geps             = eps*VFtab[nnn+2];
 792                 Heps2            = eps2*VFtab[nnn+3];
 793                 Fp               = F+Geps+Heps2;
 794                 VV               = Y+eps*Fp;
 795                 FF               = Fp+Geps+2.0*Heps2;
 796                 vcoul            = qq*VV;
 797                 fijC             = qq*FF;
 798                 vctot            = vctot + vcoul;
 799                 fscal            = -((fijC)*tabscale)*rinv13;
 800                 fscal *= hybscal;
 801                 tx               = fscal*dx13;
 802                 ty               = fscal*dy13;
 803                 tz               = fscal*dz13;
 804                 fix1             = fix1 + tx;
 805                 fiy1             = fiy1 + ty;
 806                 fiz1             = fiz1 + tz;
 807                 fjx3             = faction[j3+6] - tx;
 808                 fjy3             = faction[j3+7] - ty;
 809                 fjz3             = faction[j3+8] - tz;
 810                 qq               = qqOH;
 811                 r                = rsq21*rinv21;
 812                 rt               = r*tabscale;
 813                 n0               = rt;
 814                 eps              = rt-n0;
 815                 eps2             = eps*eps;
 816                 nnn              = 4*n0;
 817                 Y                = VFtab[nnn];
 818                 F                = VFtab[nnn+1];
 819                 Geps             = eps*VFtab[nnn+2];
 820                 Heps2            = eps2*VFtab[nnn+3];
 821                 Fp               = F+Geps+Heps2;
 822                 VV               = Y+eps*Fp;
 823                 FF               = Fp+Geps+2.0*Heps2;
 824                 vcoul            = qq*VV;
 825                 fijC             = qq*FF;
 826                 vctot            = vctot + vcoul;
 827                 fscal            = -((fijC)*tabscale)*rinv21;
 828                 fscal *= hybscal;
 829                 tx               = fscal*dx21;
 830                 ty               = fscal*dy21;
 831                 tz               = fscal*dz21;
 832                 fix2             = fix2 + tx;
 833                 fiy2             = fiy2 + ty;
 834                 fiz2             = fiz2 + tz;
 835                 fjx1             = fjx1 - tx;
 836                 fjy1             = fjy1 - ty;
 837                 fjz1             = fjz1 - tz;
 838                 qq               = qqHH;
 839                 r                = rsq22*rinv22;
 840                 rt               = r*tabscale;
 841                 n0               = rt;
 842                 eps              = rt-n0;
 843                 eps2             = eps*eps;
 844                 nnn              = 4*n0;
 845                 Y                = VFtab[nnn];
 846                 F                = VFtab[nnn+1];
 847                 Geps             = eps*VFtab[nnn+2];
 848                 Heps2            = eps2*VFtab[nnn+3];
 849                 Fp               = F+Geps+Heps2;
 850                 VV               = Y+eps*Fp;
 851                 FF               = Fp+Geps+2.0*Heps2;
 852                 vcoul            = qq*VV;
 853                 fijC             = qq*FF;
 854                 vctot            = vctot + vcoul;
 855                 fscal            = -((fijC)*tabscale)*rinv22;
 856                 fscal *= hybscal;
 857                 tx               = fscal*dx22;
 858                 ty               = fscal*dy22;
 859                 tz               = fscal*dz22;
 860                 fix2             = fix2 + tx;
 861                 fiy2             = fiy2 + ty;
 862                 fiz2             = fiz2 + tz;
 863                 fjx2             = fjx2 - tx;
 864                 fjy2             = fjy2 - ty;
 865                 fjz2             = fjz2 - tz;
 866                 qq               = qqHH;
 867                 r                = rsq23*rinv23;
 868                 rt               = r*tabscale;
 869                 n0               = rt;
 870                 eps              = rt-n0;
 871                 eps2             = eps*eps;
 872                 nnn              = 4*n0;
 873                 Y                = VFtab[nnn];
 874                 F                = VFtab[nnn+1];
 875                 Geps             = eps*VFtab[nnn+2];
 876                 Heps2            = eps2*VFtab[nnn+3];
 877                 Fp               = F+Geps+Heps2;
 878                 VV               = Y+eps*Fp;
 879                 FF               = Fp+Geps+2.0*Heps2;
 880                 vcoul            = qq*VV;
 881                 fijC             = qq*FF;
 882                 vctot            = vctot + vcoul;
 883                 fscal            = -((fijC)*tabscale)*rinv23;
 884                 fscal *= hybscal;
 885                 tx               = fscal*dx23;
 886                 ty               = fscal*dy23;
 887                 tz               = fscal*dz23;
 888                 fix2             = fix2 + tx;
 889                 fiy2             = fiy2 + ty;
 890                 fiz2             = fiz2 + tz;
 891                 fjx3             = fjx3 - tx;
 892                 fjy3             = fjy3 - ty;
 893                 fjz3             = fjz3 - tz;
 894                 qq               = qqOH;
 895                 r                = rsq31*rinv31;
 896                 rt               = r*tabscale;
 897                 n0               = rt;
 898                 eps              = rt-n0;
 899                 eps2             = eps*eps;
 900                 nnn              = 4*n0;
 901                 Y                = VFtab[nnn];
 902                 F                = VFtab[nnn+1];
 903                 Geps             = eps*VFtab[nnn+2];
 904                 Heps2            = eps2*VFtab[nnn+3];
 905                 Fp               = F+Geps+Heps2;
 906                 VV               = Y+eps*Fp;
 907                 FF               = Fp+Geps+2.0*Heps2;
 908                 vcoul            = qq*VV;
 909                 fijC             = qq*FF;
 910                 vctot            = vctot + vcoul;
 911                 fscal            = -((fijC)*tabscale)*rinv31;
 912                 fscal *= hybscal;
 913                 tx               = fscal*dx31;
 914                 ty               = fscal*dy31;
 915                 tz               = fscal*dz31;
 916                 fix3             = fix3 + tx;
 917                 fiy3             = fiy3 + ty;
 918                 fiz3             = fiz3 + tz;
 919                 faction[j3+0]    = fjx1 - tx;
 920                 faction[j3+1]    = fjy1 - ty;
 921                 faction[j3+2]    = fjz1 - tz;
 922                 qq               = qqHH;
 923                 r                = rsq32*rinv32;
 924                 rt               = r*tabscale;
 925                 n0               = rt;
 926                 eps              = rt-n0;
 927                 eps2             = eps*eps;
 928                 nnn              = 4*n0;
 929                 Y                = VFtab[nnn];
 930                 F                = VFtab[nnn+1];
 931                 Geps             = eps*VFtab[nnn+2];
 932                 Heps2            = eps2*VFtab[nnn+3];
 933                 Fp               = F+Geps+Heps2;
 934                 VV               = Y+eps*Fp;
 935                 FF               = Fp+Geps+2.0*Heps2;
 936                 vcoul            = qq*VV;
 937                 fijC             = qq*FF;
 938                 vctot            = vctot + vcoul;
 939                 fscal            = -((fijC)*tabscale)*rinv32;
 940                 fscal *= hybscal;
 941                 tx               = fscal*dx32;
 942                 ty               = fscal*dy32;
 943                 tz               = fscal*dz32;
 944                 fix3             = fix3 + tx;
 945                 fiy3             = fiy3 + ty;
 946                 fiz3             = fiz3 + tz;
 947                 faction[j3+3]    = fjx2 - tx;
 948                 faction[j3+4]    = fjy2 - ty;
 949                 faction[j3+5]    = fjz2 - tz;
 950                 qq               = qqHH;
 951                 r                = rsq33*rinv33;
 952                 rt               = r*tabscale;
 953                 n0               = rt;
 954                 eps              = rt-n0;
 955                 eps2             = eps*eps;
 956                 nnn              = 4*n0;
 957                 Y                = VFtab[nnn];
 958                 F                = VFtab[nnn+1];
 959                 Geps             = eps*VFtab[nnn+2];
 960                 Heps2            = eps2*VFtab[nnn+3];
 961                 Fp               = F+Geps+Heps2;
 962                 VV               = Y+eps*Fp;
 963                 FF               = Fp+Geps+2.0*Heps2;
 964                 vcoul            = qq*VV;
 965                 fijC             = qq*FF;
 966                 vctot            = vctot + vcoul;
 967                 fscal            = -((fijC)*tabscale)*rinv33;
 968                 fscal *= hybscal;
 969                 tx               = fscal*dx33;
 970                 ty               = fscal*dy33;
 971                 tz               = fscal*dz33;
 972                 fix3             = fix3 + tx;
 973                 fiy3             = fiy3 + ty;
 974                 fiz3             = fiz3 + tz;
 975                 faction[j3+6]    = fjx3 - tx;
 976                 faction[j3+7]    = fjy3 - ty;
 977                 faction[j3+8]    = fjz3 - tz;
 978             }
 979
 980             faction[ii3+0]   = faction[ii3+0] + fix1;
 981             faction[ii3+1]   = faction[ii3+1] + fiy1;
 982             faction[ii3+2]   = faction[ii3+2] + fiz1;
 983             faction[ii3+3]   = faction[ii3+3] + fix2;
 984             faction[ii3+4]   = faction[ii3+4] + fiy2;
 985             faction[ii3+5]   = faction[ii3+5] + fiz2;
 986             faction[ii3+6]   = faction[ii3+6] + fix3;
 987             faction[ii3+7]   = faction[ii3+7] + fiy3;
 988             faction[ii3+8]   = faction[ii3+8] + fiz3;
 989             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
 990             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
 991             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
 992             ggid             = gid[n];
 993             Vc[ggid]         = Vc[ggid] + vctot;
 994             ninner           = ninner + nj1 - nj0;
 995         }
 996
 997         nouter           = nouter + nn1 - nn0;
 998     }
 999     while (nn1<nri);
1000
1001     *outeriter       = nouter;
1002     *inneriter       = ninner;
1003 }
1004
1005