src/gromacs/gmxlib/nonbonded/nb_kernel_adress_c/nb_kernel232_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_kernel232_adress.h"
  28
  29
  30
  31 /*
  32  * Gromacs nonbonded kernel nb_kernel232_adress_cg
  33  * Coulomb interaction:     Reaction field
  34  * VdW interaction:         Tabulated
  35  * water optimization:      pairs of SPC/TIP3P interactions
  36  * Calculate forces:        yes
  37  */
  38 void nb_kernel232_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         krsq;
  89     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
  90     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
  91     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
  92     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
  93     real         jx2,jy2,jz2,fjx2,fjy2,fjz2;
  94     real         jx3,jy3,jz3,fjx3,fjy3,fjz3;
  95     real         dx11,dy11,dz11,rsq11,rinv11;
  96     real         dx12,dy12,dz12,rsq12,rinv12;
  97     real         dx13,dy13,dz13,rsq13,rinv13;
  98     real         dx21,dy21,dz21,rsq21,rinv21;
  99     real         dx22,dy22,dz22,rsq22,rinv22;
 100     real         dx23,dy23,dz23,rsq23,rinv23;
 101     real         dx31,dy31,dz31,rsq31,rinv31;
 102     real         dx32,dy32,dz32,rsq32,rinv32;
 103     real         dx33,dy33,dz33,rsq33,rinv33;
 104     real         qO,qH,qqOO,qqOH,qqHH;
 105     real         c6,c12;
 106     real         weight_cg1, weight_cg2, weight_product;
 107     real         hybscal;
 108
 109     nri              = *p_nri;
 110     ntype            = *p_ntype;
 111     nthreads         = *p_nthreads;
 112     facel            = *p_facel;
 113     krf              = *p_krf;
 114     crf              = *p_crf;
 115     tabscale         = *p_tabscale;
 116     ii               = iinr[0];
 117     qO               = charge[ii];
 118     qH               = charge[ii+1];
 119     qqOO             = facel*qO*qO;
 120     qqOH             = facel*qO*qH;
 121     qqHH             = facel*qH*qH;
 122     tj               = 2*(ntype+1)*type[ii];
 123     c6               = vdwparam[tj];
 124     c12              = vdwparam[tj+1];
 125
 126     nouter           = 0;
 127     ninner           = 0;
 128
 129     do
 130     {
 131         #ifdef GMX_THREAD_SHM_FDECOMP
 132         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 133         nn0              = *count;
 134         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 135         *count           = nn1;
 136         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 137         if(nn1>nri) nn1=nri;
 138         #else
 139         nn0 = 0;
 140         nn1 = nri;
 141         #endif
 142
 143         for(n=nn0; (n<nn1); n++)
 144         {
 145             is3              = 3*shift[n];
 146             shX              = shiftvec[is3];
 147             shY              = shiftvec[is3+1];
 148             shZ              = shiftvec[is3+2];
 149             nj0              = jindex[n];
 150             nj1              = jindex[n+1];
 151             ii               = iinr[n];
 152             ii3              = 3*ii;
 153             ix1              = shX + pos[ii3+0];
 154             iy1              = shY + pos[ii3+1];
 155             iz1              = shZ + pos[ii3+2];
 156             ix2              = shX + pos[ii3+3];
 157             iy2              = shY + pos[ii3+4];
 158             iz2              = shZ + pos[ii3+5];
 159             ix3              = shX + pos[ii3+6];
 160             iy3              = shY + pos[ii3+7];
 161             iz3              = shZ + pos[ii3+8];
 162             weight_cg1       = wf[ii];
 163             vctot            = 0;
 164             Vvdwtot          = 0;
 165             fix1             = 0;
 166             fiy1             = 0;
 167             fiz1             = 0;
 168             fix2             = 0;
 169             fiy2             = 0;
 170             fiz2             = 0;
 171             fix3             = 0;
 172             fiy3             = 0;
 173             fiz3             = 0;
 174
 175             for(k=nj0; (k<nj1); k++)
 176             {
 177                 jnr              = jjnr[k];
 178                 weight_cg2       = wf[jnr];
 179                 weight_product   = weight_cg1*weight_cg2;
 180                 if (weight_product < ALMOST_ZERO) {
 181                        hybscal = 1.0;
 182                 }
 183                 else if (weight_product >= ALMOST_ONE)
 184                 {
 185                   /* force is zero, skip this molecule */
 186                        continue;
 187                 }
 188                 else
 189                 {
 190                    hybscal = 1.0 - weight_product;
 191                 }
 192                 j3               = 3*jnr;
 193                 jx1              = pos[j3+0];
 194                 jy1              = pos[j3+1];
 195                 jz1              = pos[j3+2];
 196                 jx2              = pos[j3+3];
 197                 jy2              = pos[j3+4];
 198                 jz2              = pos[j3+5];
 199                 jx3              = pos[j3+6];
 200                 jy3              = pos[j3+7];
 201                 jz3              = pos[j3+8];
 202                 dx11             = ix1 - jx1;
 203                 dy11             = iy1 - jy1;
 204                 dz11             = iz1 - jz1;
 205                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 206                 dx12             = ix1 - jx2;
 207                 dy12             = iy1 - jy2;
 208                 dz12             = iz1 - jz2;
 209                 rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 210                 dx13             = ix1 - jx3;
 211                 dy13             = iy1 - jy3;
 212                 dz13             = iz1 - jz3;
 213                 rsq13            = dx13*dx13+dy13*dy13+dz13*dz13;
 214                 dx21             = ix2 - jx1;
 215                 dy21             = iy2 - jy1;
 216                 dz21             = iz2 - jz1;
 217                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 218                 dx22             = ix2 - jx2;
 219                 dy22             = iy2 - jy2;
 220                 dz22             = iz2 - jz2;
 221                 rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 222                 dx23             = ix2 - jx3;
 223                 dy23             = iy2 - jy3;
 224                 dz23             = iz2 - jz3;
 225                 rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 226                 dx31             = ix3 - jx1;
 227                 dy31             = iy3 - jy1;
 228                 dz31             = iz3 - jz1;
 229                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 230                 dx32             = ix3 - jx2;
 231                 dy32             = iy3 - jy2;
 232                 dz32             = iz3 - jz2;
 233                 rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 234                 dx33             = ix3 - jx3;
 235                 dy33             = iy3 - jy3;
 236                 dz33             = iz3 - jz3;
 237                 rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 238                 rinv11           = 1.0/sqrt(rsq11);
 239                 rinv12           = 1.0/sqrt(rsq12);
 240                 rinv13           = 1.0/sqrt(rsq13);
 241                 rinv21           = 1.0/sqrt(rsq21);
 242                 rinv22           = 1.0/sqrt(rsq22);
 243                 rinv23           = 1.0/sqrt(rsq23);
 244                 rinv31           = 1.0/sqrt(rsq31);
 245                 rinv32           = 1.0/sqrt(rsq32);
 246                 rinv33           = 1.0/sqrt(rsq33);
 247                 qq               = qqOO;
 248                 rinvsq           = rinv11*rinv11;
 249                 krsq             = krf*rsq11;
 250                 vcoul            = qq*(rinv11+krsq-crf);
 251                 vctot            = vctot+vcoul;
 252                 r                = rsq11*rinv11;
 253                 rt               = r*tabscale;
 254                 n0               = rt;
 255                 eps              = rt-n0;
 256                 eps2             = eps*eps;
 257                 nnn              = 8*n0;
 258                 Y                = VFtab[nnn];
 259                 F                = VFtab[nnn+1];
 260                 Geps             = eps*VFtab[nnn+2];
 261                 Heps2            = eps2*VFtab[nnn+3];
 262                 Fp               = F+Geps+Heps2;
 263                 VV               = Y+eps*Fp;
 264                 FF               = Fp+Geps+2.0*Heps2;
 265                 Vvdw6            = c6*VV;
 266                 fijD             = c6*FF;
 267                 nnn              = nnn+4;
 268                 Y                = VFtab[nnn];
 269                 F                = VFtab[nnn+1];
 270                 Geps             = eps*VFtab[nnn+2];
 271                 Heps2            = eps2*VFtab[nnn+3];
 272                 Fp               = F+Geps+Heps2;
 273                 VV               = Y+eps*Fp;
 274                 FF               = Fp+Geps+2.0*Heps2;
 275                 Vvdw12           = c12*VV;
 276                 fijR             = c12*FF;
 277                 Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12;
 278                 fscal            = (qq*(rinv11-2.0*krsq))*rinvsq-((fijD+fijR)*tabscale)*rinv11;
 279                 fscal *= hybscal;
 280                 tx               = fscal*dx11;
 281                 ty               = fscal*dy11;
 282                 tz               = fscal*dz11;
 283                 fix1             = fix1 + tx;
 284                 fiy1             = fiy1 + ty;
 285                 fiz1             = fiz1 + tz;
 286                 fjx1             = faction[j3+0] - tx;
 287                 fjy1             = faction[j3+1] - ty;
 288                 fjz1             = faction[j3+2] - tz;
 289                 qq               = qqOH;
 290                 rinvsq           = rinv12*rinv12;
 291                 krsq             = krf*rsq12;
 292                 vcoul            = qq*(rinv12+krsq-crf);
 293                 vctot            = vctot+vcoul;
 294                 fscal            = (qq*(rinv12-2.0*krsq))*rinvsq;
 295                 fscal *= hybscal;
 296                 tx               = fscal*dx12;
 297                 ty               = fscal*dy12;
 298                 tz               = fscal*dz12;
 299                 fix1             = fix1 + tx;
 300                 fiy1             = fiy1 + ty;
 301                 fiz1             = fiz1 + tz;
 302                 fjx2             = faction[j3+3] - tx;
 303                 fjy2             = faction[j3+4] - ty;
 304                 fjz2             = faction[j3+5] - tz;
 305                 qq               = qqOH;
 306                 rinvsq           = rinv13*rinv13;
 307                 krsq             = krf*rsq13;
 308                 vcoul            = qq*(rinv13+krsq-crf);
 309                 vctot            = vctot+vcoul;
 310                 fscal            = (qq*(rinv13-2.0*krsq))*rinvsq;
 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                 rinvsq           = rinv21*rinv21;
 323                 krsq             = krf*rsq21;
 324                 vcoul            = qq*(rinv21+krsq-crf);
 325                 vctot            = vctot+vcoul;
 326                 fscal            = (qq*(rinv21-2.0*krsq))*rinvsq;
 327                 fscal *= hybscal;
 328                 tx               = fscal*dx21;
 329                 ty               = fscal*dy21;
 330                 tz               = fscal*dz21;
 331                 fix2             = fix2 + tx;
 332                 fiy2             = fiy2 + ty;
 333                 fiz2             = fiz2 + tz;
 334                 fjx1             = fjx1 - tx;
 335                 fjy1             = fjy1 - ty;
 336                 fjz1             = fjz1 - tz;
 337                 qq               = qqHH;
 338                 rinvsq           = rinv22*rinv22;
 339                 krsq             = krf*rsq22;
 340                 vcoul            = qq*(rinv22+krsq-crf);
 341                 vctot            = vctot+vcoul;
 342                 fscal            = (qq*(rinv22-2.0*krsq))*rinvsq;
 343                 fscal *= hybscal;
 344                 tx               = fscal*dx22;
 345                 ty               = fscal*dy22;
 346                 tz               = fscal*dz22;
 347                 fix2             = fix2 + tx;
 348                 fiy2             = fiy2 + ty;
 349                 fiz2             = fiz2 + tz;
 350                 fjx2             = fjx2 - tx;
 351                 fjy2             = fjy2 - ty;
 352                 fjz2             = fjz2 - tz;
 353                 qq               = qqHH;
 354                 rinvsq           = rinv23*rinv23;
 355                 krsq             = krf*rsq23;
 356                 vcoul            = qq*(rinv23+krsq-crf);
 357                 vctot            = vctot+vcoul;
 358                 fscal            = (qq*(rinv23-2.0*krsq))*rinvsq;
 359                 fscal *= hybscal;
 360                 tx               = fscal*dx23;
 361                 ty               = fscal*dy23;
 362                 tz               = fscal*dz23;
 363                 fix2             = fix2 + tx;
 364                 fiy2             = fiy2 + ty;
 365                 fiz2             = fiz2 + tz;
 366                 fjx3             = fjx3 - tx;
 367                 fjy3             = fjy3 - ty;
 368                 fjz3             = fjz3 - tz;
 369                 qq               = qqOH;
 370                 rinvsq           = rinv31*rinv31;
 371                 krsq             = krf*rsq31;
 372                 vcoul            = qq*(rinv31+krsq-crf);
 373                 vctot            = vctot+vcoul;
 374                 fscal            = (qq*(rinv31-2.0*krsq))*rinvsq;
 375                 fscal *= hybscal;
 376                 tx               = fscal*dx31;
 377                 ty               = fscal*dy31;
 378                 tz               = fscal*dz31;
 379                 fix3             = fix3 + tx;
 380                 fiy3             = fiy3 + ty;
 381                 fiz3             = fiz3 + tz;
 382                 faction[j3+0]    = fjx1 - tx;
 383                 faction[j3+1]    = fjy1 - ty;
 384                 faction[j3+2]    = fjz1 - tz;
 385                 qq               = qqHH;
 386                 rinvsq           = rinv32*rinv32;
 387                 krsq             = krf*rsq32;
 388                 vcoul            = qq*(rinv32+krsq-crf);
 389                 vctot            = vctot+vcoul;
 390                 fscal            = (qq*(rinv32-2.0*krsq))*rinvsq;
 391                 fscal *= hybscal;
 392                 tx               = fscal*dx32;
 393                 ty               = fscal*dy32;
 394                 tz               = fscal*dz32;
 395                 fix3             = fix3 + tx;
 396                 fiy3             = fiy3 + ty;
 397                 fiz3             = fiz3 + tz;
 398                 faction[j3+3]    = fjx2 - tx;
 399                 faction[j3+4]    = fjy2 - ty;
 400                 faction[j3+5]    = fjz2 - tz;
 401                 qq               = qqHH;
 402                 rinvsq           = rinv33*rinv33;
 403                 krsq             = krf*rsq33;
 404                 vcoul            = qq*(rinv33+krsq-crf);
 405                 vctot            = vctot+vcoul;
 406                 fscal            = (qq*(rinv33-2.0*krsq))*rinvsq;
 407                 fscal *= hybscal;
 408                 tx               = fscal*dx33;
 409                 ty               = fscal*dy33;
 410                 tz               = fscal*dz33;
 411                 fix3             = fix3 + tx;
 412                 fiy3             = fiy3 + ty;
 413                 fiz3             = fiz3 + tz;
 414                 faction[j3+6]    = fjx3 - tx;
 415                 faction[j3+7]    = fjy3 - ty;
 416                 faction[j3+8]    = fjz3 - tz;
 417             }
 418
 419             faction[ii3+0]   = faction[ii3+0] + fix1;
 420             faction[ii3+1]   = faction[ii3+1] + fiy1;
 421             faction[ii3+2]   = faction[ii3+2] + fiz1;
 422             faction[ii3+3]   = faction[ii3+3] + fix2;
 423             faction[ii3+4]   = faction[ii3+4] + fiy2;
 424             faction[ii3+5]   = faction[ii3+5] + fiz2;
 425             faction[ii3+6]   = faction[ii3+6] + fix3;
 426             faction[ii3+7]   = faction[ii3+7] + fiy3;
 427             faction[ii3+8]   = faction[ii3+8] + fiz3;
 428             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
 429             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
 430             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
 431             ggid             = gid[n];
 432             Vc[ggid]         = Vc[ggid] + vctot;
 433             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 434             ninner           = ninner + nj1 - nj0;
 435         }
 436
 437         nouter           = nouter + nn1 - nn0;
 438     }
 439     while (nn1<nri);
 440
 441     *outeriter       = nouter;
 442     *inneriter       = ninner;
 443 }
 444
 445
 446
 447
 448
 449 /*
 450  * Gromacs nonbonded kernel nb_kernel232_adress_ex
 451  * Coulomb interaction:     Reaction field
 452  * VdW interaction:         Tabulated
 453  * water optimization:      pairs of SPC/TIP3P interactions
 454  * Calculate forces:        yes
 455  */
 456 void nb_kernel232_adress_ex(
 457                     int *           p_nri,
 458                     int *           iinr,
 459                     int *           jindex,
 460                     int *           jjnr,
 461                     int *           shift,
 462                     real *         shiftvec,
 463                     real *         fshift,
 464                     int *           gid,
 465                     real *         pos,
 466                     real *         faction,
 467                     real *         charge,
 468                     real *         p_facel,
 469                     real *         p_krf,
 470                     real *         p_crf,
 471                     real *         Vc,
 472                     int *           type,
 473                     int *           p_ntype,
 474                     real *         vdwparam,
 475                     real *         Vvdw,
 476                     real *         p_tabscale,
 477                     real *         VFtab,
 478                     real *         invsqrta,
 479                     real *         dvda,
 480                     real *         p_gbtabscale,
 481                     real *         GBtab,
 482                     int *           p_nthreads,
 483                     int *           count,
 484                     void *          mtx,
 485                     int *           outeriter,
 486                     int *           inneriter,
 487                     real           force_cap,
 488                     real *         wf)
 489 {
 490     int           nri,ntype,nthreads;
 491     real         facel,krf,crf,tabscale,gbtabscale;
 492     int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
 493     int           nn0,nn1,nouter,ninner;
 494     real         shX,shY,shZ;
 495     real         fscal,tx,ty,tz;
 496     real         rinvsq;
 497     real         qq,vcoul,vctot;
 498     int           tj;
 499     real         Vvdw6,Vvdwtot;
 500     real         Vvdw12;
 501     real         r,rt,eps,eps2;
 502     int           n0,nnn;
 503     real         Y,F,Geps,Heps2,Fp,VV;
 504     real         FF;
 505     real         fijD,fijR;
 506     real         krsq;
 507     real         ix1,iy1,iz1,fix1,fiy1,fiz1;
 508     real         ix2,iy2,iz2,fix2,fiy2,fiz2;
 509     real         ix3,iy3,iz3,fix3,fiy3,fiz3;
 510     real         jx1,jy1,jz1,fjx1,fjy1,fjz1;
 511     real         jx2,jy2,jz2,fjx2,fjy2,fjz2;
 512     real         jx3,jy3,jz3,fjx3,fjy3,fjz3;
 513     real         dx11,dy11,dz11,rsq11,rinv11;
 514     real         dx12,dy12,dz12,rsq12,rinv12;
 515     real         dx13,dy13,dz13,rsq13,rinv13;
 516     real         dx21,dy21,dz21,rsq21,rinv21;
 517     real         dx22,dy22,dz22,rsq22,rinv22;
 518     real         dx23,dy23,dz23,rsq23,rinv23;
 519     real         dx31,dy31,dz31,rsq31,rinv31;
 520     real         dx32,dy32,dz32,rsq32,rinv32;
 521     real         dx33,dy33,dz33,rsq33,rinv33;
 522     real         qO,qH,qqOO,qqOH,qqHH;
 523     real         c6,c12;
 524     real         weight_cg1, weight_cg2, weight_product;
 525     real         hybscal;
 526
 527     nri              = *p_nri;
 528     ntype            = *p_ntype;
 529     nthreads         = *p_nthreads;
 530     facel            = *p_facel;
 531     krf              = *p_krf;
 532     crf              = *p_crf;
 533     tabscale         = *p_tabscale;
 534     ii               = iinr[0];
 535     qO               = charge[ii];
 536     qH               = charge[ii+1];
 537     qqOO             = facel*qO*qO;
 538     qqOH             = facel*qO*qH;
 539     qqHH             = facel*qH*qH;
 540     tj               = 2*(ntype+1)*type[ii];
 541     c6               = vdwparam[tj];
 542     c12              = vdwparam[tj+1];
 543
 544     nouter           = 0;
 545     ninner           = 0;
 546
 547     do
 548     {
 549         #ifdef GMX_THREAD_SHM_FDECOMP
 550         tMPI_Thread_mutex_lock((tMPI_Thread_mutex_t *)mtx);
 551         nn0              = *count;
 552         nn1              = nn0+(nri-nn0)/(2*nthreads)+10;
 553         *count           = nn1;
 554         tMPI_Thread_mutex_unlock((tMPI_Thread_mutex_t *)mtx);
 555         if(nn1>nri) nn1=nri;
 556         #else
 557         nn0 = 0;
 558         nn1 = nri;
 559         #endif
 560
 561         for(n=nn0; (n<nn1); n++)
 562         {
 563             is3              = 3*shift[n];
 564             shX              = shiftvec[is3];
 565             shY              = shiftvec[is3+1];
 566             shZ              = shiftvec[is3+2];
 567             nj0              = jindex[n];
 568             nj1              = jindex[n+1];
 569             ii               = iinr[n];
 570             ii3              = 3*ii;
 571             ix1              = shX + pos[ii3+0];
 572             iy1              = shY + pos[ii3+1];
 573             iz1              = shZ + pos[ii3+2];
 574             ix2              = shX + pos[ii3+3];
 575             iy2              = shY + pos[ii3+4];
 576             iz2              = shZ + pos[ii3+5];
 577             ix3              = shX + pos[ii3+6];
 578             iy3              = shY + pos[ii3+7];
 579             iz3              = shZ + pos[ii3+8];
 580             weight_cg1       = wf[ii];
 581             vctot            = 0;
 582             Vvdwtot          = 0;
 583             fix1             = 0;
 584             fiy1             = 0;
 585             fiz1             = 0;
 586             fix2             = 0;
 587             fiy2             = 0;
 588             fiz2             = 0;
 589             fix3             = 0;
 590             fiy3             = 0;
 591             fiz3             = 0;
 592
 593             for(k=nj0; (k<nj1); k++)
 594             {
 595                 jnr              = jjnr[k];
 596                 weight_cg2       = wf[jnr];
 597                 weight_product   = weight_cg1*weight_cg2;
 598                 if (weight_product < ALMOST_ZERO) {
 599                 /* force is zero, skip this molecule */
 600                  continue;
 601                 }
 602                 else if (weight_product >= ALMOST_ONE)
 603                 {
 604                        hybscal = 1.0;
 605                 }
 606                 else
 607                 {
 608                    hybscal = weight_product;
 609                 }
 610                 j3               = 3*jnr;
 611                 jx1              = pos[j3+0];
 612                 jy1              = pos[j3+1];
 613                 jz1              = pos[j3+2];
 614                 jx2              = pos[j3+3];
 615                 jy2              = pos[j3+4];
 616                 jz2              = pos[j3+5];
 617                 jx3              = pos[j3+6];
 618                 jy3              = pos[j3+7];
 619                 jz3              = pos[j3+8];
 620                 dx11             = ix1 - jx1;
 621                 dy11             = iy1 - jy1;
 622                 dz11             = iz1 - jz1;
 623                 rsq11            = dx11*dx11+dy11*dy11+dz11*dz11;
 624                 dx12             = ix1 - jx2;
 625                 dy12             = iy1 - jy2;
 626                 dz12             = iz1 - jz2;
 627                 rsq12            = dx12*dx12+dy12*dy12+dz12*dz12;
 628                 dx13             = ix1 - jx3;
 629                 dy13             = iy1 - jy3;
 630                 dz13             = iz1 - jz3;
 631                 rsq13            = dx13*dx13+dy13*dy13+dz13*dz13;
 632                 dx21             = ix2 - jx1;
 633                 dy21             = iy2 - jy1;
 634                 dz21             = iz2 - jz1;
 635                 rsq21            = dx21*dx21+dy21*dy21+dz21*dz21;
 636                 dx22             = ix2 - jx2;
 637                 dy22             = iy2 - jy2;
 638                 dz22             = iz2 - jz2;
 639                 rsq22            = dx22*dx22+dy22*dy22+dz22*dz22;
 640                 dx23             = ix2 - jx3;
 641                 dy23             = iy2 - jy3;
 642                 dz23             = iz2 - jz3;
 643                 rsq23            = dx23*dx23+dy23*dy23+dz23*dz23;
 644                 dx31             = ix3 - jx1;
 645                 dy31             = iy3 - jy1;
 646                 dz31             = iz3 - jz1;
 647                 rsq31            = dx31*dx31+dy31*dy31+dz31*dz31;
 648                 dx32             = ix3 - jx2;
 649                 dy32             = iy3 - jy2;
 650                 dz32             = iz3 - jz2;
 651                 rsq32            = dx32*dx32+dy32*dy32+dz32*dz32;
 652                 dx33             = ix3 - jx3;
 653                 dy33             = iy3 - jy3;
 654                 dz33             = iz3 - jz3;
 655                 rsq33            = dx33*dx33+dy33*dy33+dz33*dz33;
 656                 rinv11           = 1.0/sqrt(rsq11);
 657                 rinv12           = 1.0/sqrt(rsq12);
 658                 rinv13           = 1.0/sqrt(rsq13);
 659                 rinv21           = 1.0/sqrt(rsq21);
 660                 rinv22           = 1.0/sqrt(rsq22);
 661                 rinv23           = 1.0/sqrt(rsq23);
 662                 rinv31           = 1.0/sqrt(rsq31);
 663                 rinv32           = 1.0/sqrt(rsq32);
 664                 rinv33           = 1.0/sqrt(rsq33);
 665                 qq               = qqOO;
 666                 rinvsq           = rinv11*rinv11;
 667                 krsq             = krf*rsq11;
 668                 vcoul            = qq*(rinv11+krsq-crf);
 669                 vctot            = vctot+vcoul;
 670                 r                = rsq11*rinv11;
 671                 rt               = r*tabscale;
 672                 n0               = rt;
 673                 eps              = rt-n0;
 674                 eps2             = eps*eps;
 675                 nnn              = 8*n0;
 676                 Y                = VFtab[nnn];
 677                 F                = VFtab[nnn+1];
 678                 Geps             = eps*VFtab[nnn+2];
 679                 Heps2            = eps2*VFtab[nnn+3];
 680                 Fp               = F+Geps+Heps2;
 681                 VV               = Y+eps*Fp;
 682                 FF               = Fp+Geps+2.0*Heps2;
 683                 Vvdw6            = c6*VV;
 684                 fijD             = c6*FF;
 685                 nnn              = nnn+4;
 686                 Y                = VFtab[nnn];
 687                 F                = VFtab[nnn+1];
 688                 Geps             = eps*VFtab[nnn+2];
 689                 Heps2            = eps2*VFtab[nnn+3];
 690                 Fp               = F+Geps+Heps2;
 691                 VV               = Y+eps*Fp;
 692                 FF               = Fp+Geps+2.0*Heps2;
 693                 Vvdw12           = c12*VV;
 694                 fijR             = c12*FF;
 695                 Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12;
 696                 fscal            = (qq*(rinv11-2.0*krsq))*rinvsq-((fijD+fijR)*tabscale)*rinv11;
 697                 fscal *= hybscal;
 698                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 699                 fscal=force_cap*fscal/fabs(fscal);
 700                 }
 701                 tx               = fscal*dx11;
 702                 ty               = fscal*dy11;
 703                 tz               = fscal*dz11;
 704                 fix1             = fix1 + tx;
 705                 fiy1             = fiy1 + ty;
 706                 fiz1             = fiz1 + tz;
 707                 fjx1             = faction[j3+0] - tx;
 708                 fjy1             = faction[j3+1] - ty;
 709                 fjz1             = faction[j3+2] - tz;
 710                 qq               = qqOH;
 711                 rinvsq           = rinv12*rinv12;
 712                 krsq             = krf*rsq12;
 713                 vcoul            = qq*(rinv12+krsq-crf);
 714                 vctot            = vctot+vcoul;
 715                 fscal            = (qq*(rinv12-2.0*krsq))*rinvsq;
 716                 fscal *= hybscal;
 717                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 718                 fscal=force_cap*fscal/fabs(fscal);
 719                 }
 720                 tx               = fscal*dx12;
 721                 ty               = fscal*dy12;
 722                 tz               = fscal*dz12;
 723                 fix1             = fix1 + tx;
 724                 fiy1             = fiy1 + ty;
 725                 fiz1             = fiz1 + tz;
 726                 fjx2             = faction[j3+3] - tx;
 727                 fjy2             = faction[j3+4] - ty;
 728                 fjz2             = faction[j3+5] - tz;
 729                 qq               = qqOH;
 730                 rinvsq           = rinv13*rinv13;
 731                 krsq             = krf*rsq13;
 732                 vcoul            = qq*(rinv13+krsq-crf);
 733                 vctot            = vctot+vcoul;
 734                 fscal            = (qq*(rinv13-2.0*krsq))*rinvsq;
 735                 fscal *= hybscal;
 736                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 737                 fscal=force_cap*fscal/fabs(fscal);
 738                 }
 739                 tx               = fscal*dx13;
 740                 ty               = fscal*dy13;
 741                 tz               = fscal*dz13;
 742                 fix1             = fix1 + tx;
 743                 fiy1             = fiy1 + ty;
 744                 fiz1             = fiz1 + tz;
 745                 fjx3             = faction[j3+6] - tx;
 746                 fjy3             = faction[j3+7] - ty;
 747                 fjz3             = faction[j3+8] - tz;
 748                 qq               = qqOH;
 749                 rinvsq           = rinv21*rinv21;
 750                 krsq             = krf*rsq21;
 751                 vcoul            = qq*(rinv21+krsq-crf);
 752                 vctot            = vctot+vcoul;
 753                 fscal            = (qq*(rinv21-2.0*krsq))*rinvsq;
 754                 fscal *= hybscal;
 755                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 756                 fscal=force_cap*fscal/fabs(fscal);
 757                 }
 758                 tx               = fscal*dx21;
 759                 ty               = fscal*dy21;
 760                 tz               = fscal*dz21;
 761                 fix2             = fix2 + tx;
 762                 fiy2             = fiy2 + ty;
 763                 fiz2             = fiz2 + tz;
 764                 fjx1             = fjx1 - tx;
 765                 fjy1             = fjy1 - ty;
 766                 fjz1             = fjz1 - tz;
 767                 qq               = qqHH;
 768                 rinvsq           = rinv22*rinv22;
 769                 krsq             = krf*rsq22;
 770                 vcoul            = qq*(rinv22+krsq-crf);
 771                 vctot            = vctot+vcoul;
 772                 fscal            = (qq*(rinv22-2.0*krsq))*rinvsq;
 773                 fscal *= hybscal;
 774                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 775                 fscal=force_cap*fscal/fabs(fscal);
 776                 }
 777                 tx               = fscal*dx22;
 778                 ty               = fscal*dy22;
 779                 tz               = fscal*dz22;
 780                 fix2             = fix2 + tx;
 781                 fiy2             = fiy2 + ty;
 782                 fiz2             = fiz2 + tz;
 783                 fjx2             = fjx2 - tx;
 784                 fjy2             = fjy2 - ty;
 785                 fjz2             = fjz2 - tz;
 786                 qq               = qqHH;
 787                 rinvsq           = rinv23*rinv23;
 788                 krsq             = krf*rsq23;
 789                 vcoul            = qq*(rinv23+krsq-crf);
 790                 vctot            = vctot+vcoul;
 791                 fscal            = (qq*(rinv23-2.0*krsq))*rinvsq;
 792                 fscal *= hybscal;
 793                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 794                 fscal=force_cap*fscal/fabs(fscal);
 795                 }
 796                 tx               = fscal*dx23;
 797                 ty               = fscal*dy23;
 798                 tz               = fscal*dz23;
 799                 fix2             = fix2 + tx;
 800                 fiy2             = fiy2 + ty;
 801                 fiz2             = fiz2 + tz;
 802                 fjx3             = fjx3 - tx;
 803                 fjy3             = fjy3 - ty;
 804                 fjz3             = fjz3 - tz;
 805                 qq               = qqOH;
 806                 rinvsq           = rinv31*rinv31;
 807                 krsq             = krf*rsq31;
 808                 vcoul            = qq*(rinv31+krsq-crf);
 809                 vctot            = vctot+vcoul;
 810                 fscal            = (qq*(rinv31-2.0*krsq))*rinvsq;
 811                 fscal *= hybscal;
 812                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 813                 fscal=force_cap*fscal/fabs(fscal);
 814                 }
 815                 tx               = fscal*dx31;
 816                 ty               = fscal*dy31;
 817                 tz               = fscal*dz31;
 818                 fix3             = fix3 + tx;
 819                 fiy3             = fiy3 + ty;
 820                 fiz3             = fiz3 + tz;
 821                 faction[j3+0]    = fjx1 - tx;
 822                 faction[j3+1]    = fjy1 - ty;
 823                 faction[j3+2]    = fjz1 - tz;
 824                 qq               = qqHH;
 825                 rinvsq           = rinv32*rinv32;
 826                 krsq             = krf*rsq32;
 827                 vcoul            = qq*(rinv32+krsq-crf);
 828                 vctot            = vctot+vcoul;
 829                 fscal            = (qq*(rinv32-2.0*krsq))*rinvsq;
 830                 fscal *= hybscal;
 831                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 832                 fscal=force_cap*fscal/fabs(fscal);
 833                 }
 834                 tx               = fscal*dx32;
 835                 ty               = fscal*dy32;
 836                 tz               = fscal*dz32;
 837                 fix3             = fix3 + tx;
 838                 fiy3             = fiy3 + ty;
 839                 fiz3             = fiz3 + tz;
 840                 faction[j3+3]    = fjx2 - tx;
 841                 faction[j3+4]    = fjy2 - ty;
 842                 faction[j3+5]    = fjz2 - tz;
 843                 qq               = qqHH;
 844                 rinvsq           = rinv33*rinv33;
 845                 krsq             = krf*rsq33;
 846                 vcoul            = qq*(rinv33+krsq-crf);
 847                 vctot            = vctot+vcoul;
 848                 fscal            = (qq*(rinv33-2.0*krsq))*rinvsq;
 849                 fscal *= hybscal;
 850                 if(force_cap>0 && (fabs(fscal)> force_cap)){
 851                 fscal=force_cap*fscal/fabs(fscal);
 852                 }
 853                 tx               = fscal*dx33;
 854                 ty               = fscal*dy33;
 855                 tz               = fscal*dz33;
 856                 fix3             = fix3 + tx;
 857                 fiy3             = fiy3 + ty;
 858                 fiz3             = fiz3 + tz;
 859                 faction[j3+6]    = fjx3 - tx;
 860                 faction[j3+7]    = fjy3 - ty;
 861                 faction[j3+8]    = fjz3 - tz;
 862             }
 863
 864             faction[ii3+0]   = faction[ii3+0] + fix1;
 865             faction[ii3+1]   = faction[ii3+1] + fiy1;
 866             faction[ii3+2]   = faction[ii3+2] + fiz1;
 867             faction[ii3+3]   = faction[ii3+3] + fix2;
 868             faction[ii3+4]   = faction[ii3+4] + fiy2;
 869             faction[ii3+5]   = faction[ii3+5] + fiz2;
 870             faction[ii3+6]   = faction[ii3+6] + fix3;
 871             faction[ii3+7]   = faction[ii3+7] + fiy3;
 872             faction[ii3+8]   = faction[ii3+8] + fiz3;
 873             fshift[is3]      = fshift[is3]+fix1+fix2+fix3;
 874             fshift[is3+1]    = fshift[is3+1]+fiy1+fiy2+fiy3;
 875             fshift[is3+2]    = fshift[is3+2]+fiz1+fiz2+fiz3;
 876             ggid             = gid[n];
 877             Vc[ggid]         = Vc[ggid] + vctot;
 878             Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
 879             ninner           = ninner + nj1 - nj0;
 880         }
 881
 882         nouter           = nouter + nn1 - nn0;
 883     }
 884     while (nn1<nri);
 885
 886     *outeriter       = nouter;
 887     *inneriter       = ninner;
 888 }
 889
 890