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