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