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