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