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