src/gmxlib/nonbonded/nb_kernel_power6/pwr6kernel333.F

   1 C
   2 C                This source code is part of
   3 C
   4 C                 G   R   O   M   A   C   S
   5 C
   6 C Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   7 C Copyright (c) 2001-2009, The GROMACS Development Team
   8 C
   9 C Gromacs is a library for molecular simulation and trajectory analysis,
  10 C written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
  11 C a full list of developers and information, check out http://www.gromacs.org
  12 C
  13 C This program is free software; you can redistribute it and/or modify it under
  14 C the terms of the GNU Lesser General Public License as published by the Free
  15 C Software Foundation; either version 2 of the License, or (at your option) any
  16 C later version.
  17 C As a special exception, you may use this file as part of a free software
  18 C library without restriction.  Specifically, if other files instantiate
  19 C templates or use macros or inline functions from this file, or you compile
  20 C this file and link it with other files to produce an executable, this
  21 C file does not by itself cause the resulting executable to be covered by
  22 C the GNU Lesser General Public License.
  23 C
  24 C In plain-speak: do not worry about classes/macros/templates either - only
  25 C changes to the library have to be LGPL, not an application linking with it.
  26 C
  27 C To help fund GROMACS development, we humbly ask that you cite
  28 C the papers people have written on it - you can find them on the website!
  29 C
  30
  31 #ifdef HAVE_CONFIG_H
  32 #  include<config.h>
  33 #endif
  34
  35 #ifdef GMX_DOUBLE
  36 #  define gmxreal real*8
  37 #else
  38 #  define gmxreal real*4
  39 #endif
  40
  41
  42
  43 C
  44 C Gromacs nonbonded kernel pwr6kernel333
  45 C Coulomb interaction:     Tabulated
  46 C VdW interaction:         Tabulated
  47 C water optimization:      TIP4P - other atoms
  48 C Calculate forces:        yes
  49 C
  50       subroutine pwr6kernel333(
  51      &                          nri,
  52      &                          iinr,
  53      &                          jindex,
  54      &                          jjnr,
  55      &                          shift,
  56      &                          shiftvec,
  57      &                          fshift,
  58      &                          gid,
  59      &                          pos,
  60      &                          faction,
  61      &                          charge,
  62      &                          facel,
  63      &                          krf,
  64      &                          crf,
  65      &                          Vc,
  66      &                          type,
  67      &                          ntype,
  68      &                          vdwparam,
  69      &                          Vvdw,
  70      &                          tabscale,
  71      &                          VFtab,
  72      &                          invsqrta,
  73      &                          dvda,
  74      &                          gbtabscale,
  75      &                          GBtab,
  76      &                          nthreads,
  77      &                          count,
  78      &                          mtx,
  79      &                          outeriter,
  80      &                          inneriter,
  81      &                          work)
  82       implicit      none
  83       integer*4     nri,iinr(*),jindex(*),jjnr(*),shift(*)
  84       gmxreal       shiftvec(*),fshift(*),pos(*),faction(*)
  85       integer*4     gid(*),type(*),ntype
  86       gmxreal       charge(*),facel,krf,crf,Vc(*),vdwparam(*)
  87       gmxreal       Vvdw(*),tabscale,VFtab(*)
  88       gmxreal       invsqrta(*),dvda(*),gbtabscale,GBtab(*)
  89       integer*4     nthreads,count,mtx,outeriter,inneriter
  90       gmxreal       work(*)
  91
  92       integer*4     n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
  93       integer*4     nn0,nn1,nouter,ninner
  94       gmxreal       shX,shY,shZ
  95       gmxreal       fscal,tx,ty,tz
  96       gmxreal       jq
  97       gmxreal       qq,vcoul,vctot
  98       integer*4     nti
  99       integer*4     tj
 100       gmxreal       Vvdw6,Vvdwtot
 101       gmxreal       Vvdw12
 102       gmxreal       r,rt,eps,eps2
 103       integer*4     n0,nnn
 104       gmxreal       Y,F,Geps,Heps2,Fp,VV
 105       gmxreal       FF
 106       gmxreal       fijC
 107       gmxreal       fijD,fijR
 108       gmxreal       ix1,iy1,iz1,fix1,fiy1,fiz1
 109       gmxreal       ix2,iy2,iz2,fix2,fiy2,fiz2
 110       gmxreal       ix3,iy3,iz3,fix3,fiy3,fiz3
 111       gmxreal       ix4,iy4,iz4,fix4,fiy4,fiz4
 112       gmxreal       jx1,jy1,jz1,fjx1,fjy1,fjz1
 113       gmxreal       dx11,dy11,dz11,rsq11,rinv11
 114       gmxreal       dx21,dy21,dz21,rsq21,rinv21
 115       gmxreal       dx31,dy31,dz31,rsq31,rinv31
 116       gmxreal       dx41,dy41,dz41,rsq41,rinv41
 117       gmxreal       qH,qM
 118       gmxreal       c6,c12
 119
 120
 121 C    Initialize water data
 122       ii               = iinr(1)+1
 123       qH               = facel*charge(ii+1)
 124       qM               = facel*charge(ii+3)
 125       nti              = 2*ntype*type(ii)
 126
 127
 128 C    Reset outer and inner iteration counters
 129       nouter           = 0
 130       ninner           = 0
 131
 132 C    Loop over thread workunits
 133    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 134         if(nn1.gt.nri) nn1=nri
 135
 136 C      Start outer loop over neighborlists
 137
 138         do n=nn0+1,nn1
 139
 140 C        Load shift vector for this list
 141           is3              = 3*shift(n)+1
 142           shX              = shiftvec(is3)
 143           shY              = shiftvec(is3+1)
 144           shZ              = shiftvec(is3+2)
 145
 146 C        Load limits for loop over neighbors
 147           nj0              = jindex(n)+1
 148           nj1              = jindex(n+1)
 149
 150 C        Get outer coordinate index
 151           ii               = iinr(n)+1
 152           ii3              = 3*ii-2
 153
 154 C        Load i atom data, add shift vector
 155           ix1              = shX + pos(ii3+0)
 156           iy1              = shY + pos(ii3+1)
 157           iz1              = shZ + pos(ii3+2)
 158           ix2              = shX + pos(ii3+3)
 159           iy2              = shY + pos(ii3+4)
 160           iz2              = shZ + pos(ii3+5)
 161           ix3              = shX + pos(ii3+6)
 162           iy3              = shY + pos(ii3+7)
 163           iz3              = shZ + pos(ii3+8)
 164           ix4              = shX + pos(ii3+9)
 165           iy4              = shY + pos(ii3+10)
 166           iz4              = shZ + pos(ii3+11)
 167
 168 C        Zero the potential energy for this list
 169           vctot            = 0
 170           Vvdwtot          = 0
 171
 172 C        Clear i atom forces
 173           fix1             = 0
 174           fiy1             = 0
 175           fiz1             = 0
 176           fix2             = 0
 177           fiy2             = 0
 178           fiz2             = 0
 179           fix3             = 0
 180           fiy3             = 0
 181           fiz3             = 0
 182           fix4             = 0
 183           fiy4             = 0
 184           fiz4             = 0
 185
 186           do k=nj0,nj1
 187
 188 C          Get j neighbor index, and coordinate index
 189             jnr              = jjnr(k)+1
 190             j3               = 3*jnr-2
 191
 192 C          load j atom coordinates
 193             jx1              = pos(j3+0)
 194             jy1              = pos(j3+1)
 195             jz1              = pos(j3+2)
 196
 197 C          Calculate distance
 198             dx11             = ix1 - jx1
 199             dy11             = iy1 - jy1
 200             dz11             = iz1 - jz1
 201             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11
 202             dx21             = ix2 - jx1
 203             dy21             = iy2 - jy1
 204             dz21             = iz2 - jz1
 205             rsq21            = dx21*dx21+dy21*dy21+dz21*dz21
 206             dx31             = ix3 - jx1
 207             dy31             = iy3 - jy1
 208             dz31             = iz3 - jz1
 209             rsq31            = dx31*dx31+dy31*dy31+dz31*dz31
 210             dx41             = ix4 - jx1
 211             dy41             = iy4 - jy1
 212             dz41             = iz4 - jz1
 213             rsq41            = dx41*dx41+dy41*dy41+dz41*dz41
 214
 215 C          Calculate 1/r and 1/r2
 216
 217 C          PowerPC intrinsics 1/sqrt lookup table
 218 #ifndef GMX_BLUEGENE
 219             rinv11           = frsqrtes(rsq11)
 220 #else
 221             rinv11           = frsqrte(dble(rsq11))
 222 #endif
 223             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 224      &  *rinv11)))
 225 #ifdef GMX_DOUBLE
 226             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 227      &  *rinv11)))
 228 #endif
 229
 230 C          PowerPC intrinsics 1/sqrt lookup table
 231 #ifndef GMX_BLUEGENE
 232             rinv21           = frsqrtes(rsq21)
 233 #else
 234             rinv21           = frsqrte(dble(rsq21))
 235 #endif
 236             rinv21           = (0.5*rinv21*(3.0-((rsq21*rinv21)
 237      &  *rinv21)))
 238 #ifdef GMX_DOUBLE
 239             rinv21           = (0.5*rinv21*(3.0-((rsq21*rinv21)
 240      &  *rinv21)))
 241 #endif
 242
 243 C          PowerPC intrinsics 1/sqrt lookup table
 244 #ifndef GMX_BLUEGENE
 245             rinv31           = frsqrtes(rsq31)
 246 #else
 247             rinv31           = frsqrte(dble(rsq31))
 248 #endif
 249             rinv31           = (0.5*rinv31*(3.0-((rsq31*rinv31)
 250      &  *rinv31)))
 251 #ifdef GMX_DOUBLE
 252             rinv31           = (0.5*rinv31*(3.0-((rsq31*rinv31)
 253      &  *rinv31)))
 254 #endif
 255
 256 C          PowerPC intrinsics 1/sqrt lookup table
 257 #ifndef GMX_BLUEGENE
 258             rinv41           = frsqrtes(rsq41)
 259 #else
 260             rinv41           = frsqrte(dble(rsq41))
 261 #endif
 262             rinv41           = (0.5*rinv41*(3.0-((rsq41*rinv41)
 263      &  *rinv41)))
 264 #ifdef GMX_DOUBLE
 265             rinv41           = (0.5*rinv41*(3.0-((rsq41*rinv41)
 266      &  *rinv41)))
 267 #endif
 268
 269 C          Load parameters for j atom
 270             tj               = nti+2*type(jnr)+1
 271             c6               = vdwparam(tj)
 272             c12              = vdwparam(tj+1)
 273
 274 C          Calculate table index
 275             r                = rsq11*rinv11
 276
 277 C          Calculate table index
 278             rt               = r*tabscale
 279             n0               = rt
 280             eps              = rt-n0
 281             eps2             = eps*eps
 282             nnn              = 12*n0+1
 283
 284 C          Tabulated VdW interaction - dispersion
 285             nnn              = nnn+4
 286             Y                = VFtab(nnn)
 287             F                = VFtab(nnn+1)
 288             Geps             = eps*VFtab(nnn+2)
 289             Heps2            = eps2*VFtab(nnn+3)
 290             Fp               = F+Geps+Heps2
 291             VV               = Y+eps*Fp
 292             FF               = Fp+Geps+2.0*Heps2
 293             Vvdw6            = c6*VV
 294             fijD             = c6*FF
 295
 296 C          Tabulated VdW interaction - repulsion
 297             nnn              = nnn+4
 298             Y                = VFtab(nnn)
 299             F                = VFtab(nnn+1)
 300             Geps             = eps*VFtab(nnn+2)
 301             Heps2            = eps2*VFtab(nnn+3)
 302             Fp               = F+Geps+Heps2
 303             VV               = Y+eps*Fp
 304             FF               = Fp+Geps+2.0*Heps2
 305             Vvdw12           = c12*VV
 306             fijR             = c12*FF
 307             Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12
 308             fscal            = -((fijD+fijR)*tabscale)*rinv11
 309
 310 C          Calculate temporary vectorial force
 311             tx               = fscal*dx11
 312             ty               = fscal*dy11
 313             tz               = fscal*dz11
 314
 315 C          Increment i atom force
 316             fix1             = fix1 + tx
 317             fiy1             = fiy1 + ty
 318             fiz1             = fiz1 + tz
 319
 320 C          Decrement j atom force
 321             fjx1             = faction(j3+0) - tx
 322             fjy1             = faction(j3+1) - ty
 323             fjz1             = faction(j3+2) - tz
 324
 325 C          Load parameters for j atom
 326             jq               = charge(jnr+0)
 327             qq               = qH*jq
 328
 329 C          Calculate table index
 330             r                = rsq21*rinv21
 331
 332 C          Calculate table index
 333             rt               = r*tabscale
 334             n0               = rt
 335             eps              = rt-n0
 336             eps2             = eps*eps
 337             nnn              = 12*n0+1
 338
 339 C          Tabulated coulomb interaction
 340             Y                = VFtab(nnn)
 341             F                = VFtab(nnn+1)
 342             Geps             = eps*VFtab(nnn+2)
 343             Heps2            = eps2*VFtab(nnn+3)
 344             Fp               = F+Geps+Heps2
 345             VV               = Y+eps*Fp
 346             FF               = Fp+Geps+2.0*Heps2
 347             vcoul            = qq*VV
 348             fijC             = qq*FF
 349             vctot            = vctot + vcoul
 350             fscal            = -((fijC)*tabscale)*rinv21
 351
 352 C          Calculate temporary vectorial force
 353             tx               = fscal*dx21
 354             ty               = fscal*dy21
 355             tz               = fscal*dz21
 356
 357 C          Increment i atom force
 358             fix2             = fix2 + tx
 359             fiy2             = fiy2 + ty
 360             fiz2             = fiz2 + tz
 361
 362 C          Decrement j atom force
 363             fjx1             = fjx1 - tx
 364             fjy1             = fjy1 - ty
 365             fjz1             = fjz1 - tz
 366
 367 C          Load parameters for j atom
 368
 369 C          Calculate table index
 370             r                = rsq31*rinv31
 371
 372 C          Calculate table index
 373             rt               = r*tabscale
 374             n0               = rt
 375             eps              = rt-n0
 376             eps2             = eps*eps
 377             nnn              = 12*n0+1
 378
 379 C          Tabulated coulomb interaction
 380             Y                = VFtab(nnn)
 381             F                = VFtab(nnn+1)
 382             Geps             = eps*VFtab(nnn+2)
 383             Heps2            = eps2*VFtab(nnn+3)
 384             Fp               = F+Geps+Heps2
 385             VV               = Y+eps*Fp
 386             FF               = Fp+Geps+2.0*Heps2
 387             vcoul            = qq*VV
 388             fijC             = qq*FF
 389             vctot            = vctot + vcoul
 390             fscal            = -((fijC)*tabscale)*rinv31
 391
 392 C          Calculate temporary vectorial force
 393             tx               = fscal*dx31
 394             ty               = fscal*dy31
 395             tz               = fscal*dz31
 396
 397 C          Increment i atom force
 398             fix3             = fix3 + tx
 399             fiy3             = fiy3 + ty
 400             fiz3             = fiz3 + tz
 401
 402 C          Decrement j atom force
 403             fjx1             = fjx1 - tx
 404             fjy1             = fjy1 - ty
 405             fjz1             = fjz1 - tz
 406
 407 C          Load parameters for j atom
 408             qq               = qM*jq
 409
 410 C          Calculate table index
 411             r                = rsq41*rinv41
 412
 413 C          Calculate table index
 414             rt               = r*tabscale
 415             n0               = rt
 416             eps              = rt-n0
 417             eps2             = eps*eps
 418             nnn              = 12*n0+1
 419
 420 C          Tabulated coulomb interaction
 421             Y                = VFtab(nnn)
 422             F                = VFtab(nnn+1)
 423             Geps             = eps*VFtab(nnn+2)
 424             Heps2            = eps2*VFtab(nnn+3)
 425             Fp               = F+Geps+Heps2
 426             VV               = Y+eps*Fp
 427             FF               = Fp+Geps+2.0*Heps2
 428             vcoul            = qq*VV
 429             fijC             = qq*FF
 430             vctot            = vctot + vcoul
 431             fscal            = -((fijC)*tabscale)*rinv41
 432
 433 C          Calculate temporary vectorial force
 434             tx               = fscal*dx41
 435             ty               = fscal*dy41
 436             tz               = fscal*dz41
 437
 438 C          Increment i atom force
 439             fix4             = fix4 + tx
 440             fiy4             = fiy4 + ty
 441             fiz4             = fiz4 + tz
 442
 443 C          Decrement j atom force
 444             faction(j3+0)    = fjx1 - tx
 445             faction(j3+1)    = fjy1 - ty
 446             faction(j3+2)    = fjz1 - tz
 447
 448 C          Inner loop uses 183 flops/iteration
 449           end do
 450
 451
 452 C        Add i forces to mem and shifted force list
 453           faction(ii3+0)   = faction(ii3+0) + fix1
 454           faction(ii3+1)   = faction(ii3+1) + fiy1
 455           faction(ii3+2)   = faction(ii3+2) + fiz1
 456           faction(ii3+3)   = faction(ii3+3) + fix2
 457           faction(ii3+4)   = faction(ii3+4) + fiy2
 458           faction(ii3+5)   = faction(ii3+5) + fiz2
 459           faction(ii3+6)   = faction(ii3+6) + fix3
 460           faction(ii3+7)   = faction(ii3+7) + fiy3
 461           faction(ii3+8)   = faction(ii3+8) + fiz3
 462           faction(ii3+9)   = faction(ii3+9) + fix4
 463           faction(ii3+10)  = faction(ii3+10) + fiy4
 464           faction(ii3+11)  = faction(ii3+11) + fiz4
 465           fshift(is3)      = fshift(is3)+fix1+fix2+fix3+fix4
 466           fshift(is3+1)    = fshift(is3+1)+fiy1+fiy2+fiy3+fiy4
 467           fshift(is3+2)    = fshift(is3+2)+fiz1+fiz2+fiz3+fiz4
 468
 469 C        Add potential energies to the group for this list
 470           ggid             = gid(n)+1
 471           Vc(ggid)         = Vc(ggid) + vctot
 472           Vvdw(ggid)       = Vvdw(ggid) + Vvdwtot
 473
 474 C        Increment number of inner iterations
 475           ninner           = ninner + nj1 - nj0
 476
 477 C        Outer loop uses 38 flops/iteration
 478         end do
 479
 480
 481 C      Increment number of outer iterations
 482         nouter           = nouter + nn1 - nn0
 483       if(nn1.lt.nri) goto 10
 484
 485 C    Write outer/inner iteration count to pointers
 486       outeriter        = nouter
 487       inneriter        = ninner
 488       return
 489       end
 490
 491
 492
 493
 494
 495
 496 C
 497 C Gromacs nonbonded kernel pwr6kernel333nf
 498 C Coulomb interaction:     Tabulated
 499 C VdW interaction:         Tabulated
 500 C water optimization:      TIP4P - other atoms
 501 C Calculate forces:        no
 502 C
 503       subroutine pwr6kernel333nf(
 504      &                          nri,
 505      &                          iinr,
 506      &                          jindex,
 507      &                          jjnr,
 508      &                          shift,
 509      &                          shiftvec,
 510      &                          fshift,
 511      &                          gid,
 512      &                          pos,
 513      &                          faction,
 514      &                          charge,
 515      &                          facel,
 516      &                          krf,
 517      &                          crf,
 518      &                          Vc,
 519      &                          type,
 520      &                          ntype,
 521      &                          vdwparam,
 522      &                          Vvdw,
 523      &                          tabscale,
 524      &                          VFtab,
 525      &                          invsqrta,
 526      &                          dvda,
 527      &                          gbtabscale,
 528      &                          GBtab,
 529      &                          nthreads,
 530      &                          count,
 531      &                          mtx,
 532      &                          outeriter,
 533      &                          inneriter,
 534      &                          work)
 535       implicit      none
 536       integer*4     nri,iinr(*),jindex(*),jjnr(*),shift(*)
 537       gmxreal       shiftvec(*),fshift(*),pos(*),faction(*)
 538       integer*4     gid(*),type(*),ntype
 539       gmxreal       charge(*),facel,krf,crf,Vc(*),vdwparam(*)
 540       gmxreal       Vvdw(*),tabscale,VFtab(*)
 541       gmxreal       invsqrta(*),dvda(*),gbtabscale,GBtab(*)
 542       integer*4     nthreads,count,mtx,outeriter,inneriter
 543       gmxreal       work(*)
 544
 545       integer*4     n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
 546       integer*4     nn0,nn1,nouter,ninner
 547       gmxreal       shX,shY,shZ
 548       gmxreal       jq
 549       gmxreal       qq,vcoul,vctot
 550       integer*4     nti
 551       integer*4     tj
 552       gmxreal       Vvdw6,Vvdwtot
 553       gmxreal       Vvdw12
 554       gmxreal       r,rt,eps,eps2
 555       integer*4     n0,nnn
 556       gmxreal       Y,F,Geps,Heps2,Fp,VV
 557       gmxreal       ix1,iy1,iz1
 558       gmxreal       ix2,iy2,iz2
 559       gmxreal       ix3,iy3,iz3
 560       gmxreal       ix4,iy4,iz4
 561       gmxreal       jx1,jy1,jz1
 562       gmxreal       dx11,dy11,dz11,rsq11,rinv11
 563       gmxreal       dx21,dy21,dz21,rsq21,rinv21
 564       gmxreal       dx31,dy31,dz31,rsq31,rinv31
 565       gmxreal       dx41,dy41,dz41,rsq41,rinv41
 566       gmxreal       qH,qM
 567       gmxreal       c6,c12
 568
 569
 570 C    Initialize water data
 571       ii               = iinr(1)+1
 572       qH               = facel*charge(ii+1)
 573       qM               = facel*charge(ii+3)
 574       nti              = 2*ntype*type(ii)
 575
 576
 577 C    Reset outer and inner iteration counters
 578       nouter           = 0
 579       ninner           = 0
 580
 581 C    Loop over thread workunits
 582    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 583         if(nn1.gt.nri) nn1=nri
 584
 585 C      Start outer loop over neighborlists
 586
 587         do n=nn0+1,nn1
 588
 589 C        Load shift vector for this list
 590           is3              = 3*shift(n)+1
 591           shX              = shiftvec(is3)
 592           shY              = shiftvec(is3+1)
 593           shZ              = shiftvec(is3+2)
 594
 595 C        Load limits for loop over neighbors
 596           nj0              = jindex(n)+1
 597           nj1              = jindex(n+1)
 598
 599 C        Get outer coordinate index
 600           ii               = iinr(n)+1
 601           ii3              = 3*ii-2
 602
 603 C        Load i atom data, add shift vector
 604           ix1              = shX + pos(ii3+0)
 605           iy1              = shY + pos(ii3+1)
 606           iz1              = shZ + pos(ii3+2)
 607           ix2              = shX + pos(ii3+3)
 608           iy2              = shY + pos(ii3+4)
 609           iz2              = shZ + pos(ii3+5)
 610           ix3              = shX + pos(ii3+6)
 611           iy3              = shY + pos(ii3+7)
 612           iz3              = shZ + pos(ii3+8)
 613           ix4              = shX + pos(ii3+9)
 614           iy4              = shY + pos(ii3+10)
 615           iz4              = shZ + pos(ii3+11)
 616
 617 C        Zero the potential energy for this list
 618           vctot            = 0
 619           Vvdwtot          = 0
 620
 621 C        Clear i atom forces
 622
 623           do k=nj0,nj1
 624
 625 C          Get j neighbor index, and coordinate index
 626             jnr              = jjnr(k)+1
 627             j3               = 3*jnr-2
 628
 629 C          load j atom coordinates
 630             jx1              = pos(j3+0)
 631             jy1              = pos(j3+1)
 632             jz1              = pos(j3+2)
 633
 634 C          Calculate distance
 635             dx11             = ix1 - jx1
 636             dy11             = iy1 - jy1
 637             dz11             = iz1 - jz1
 638             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11
 639             dx21             = ix2 - jx1
 640             dy21             = iy2 - jy1
 641             dz21             = iz2 - jz1
 642             rsq21            = dx21*dx21+dy21*dy21+dz21*dz21
 643             dx31             = ix3 - jx1
 644             dy31             = iy3 - jy1
 645             dz31             = iz3 - jz1
 646             rsq31            = dx31*dx31+dy31*dy31+dz31*dz31
 647             dx41             = ix4 - jx1
 648             dy41             = iy4 - jy1
 649             dz41             = iz4 - jz1
 650             rsq41            = dx41*dx41+dy41*dy41+dz41*dz41
 651
 652 C          Calculate 1/r and 1/r2
 653
 654 C          PowerPC intrinsics 1/sqrt lookup table
 655 #ifndef GMX_BLUEGENE
 656             rinv11           = frsqrtes(rsq11)
 657 #else
 658             rinv11           = frsqrte(dble(rsq11))
 659 #endif
 660             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 661      &  *rinv11)))
 662 #ifdef GMX_DOUBLE
 663             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 664      &  *rinv11)))
 665 #endif
 666
 667 C          PowerPC intrinsics 1/sqrt lookup table
 668 #ifndef GMX_BLUEGENE
 669             rinv21           = frsqrtes(rsq21)
 670 #else
 671             rinv21           = frsqrte(dble(rsq21))
 672 #endif
 673             rinv21           = (0.5*rinv21*(3.0-((rsq21*rinv21)
 674      &  *rinv21)))
 675 #ifdef GMX_DOUBLE
 676             rinv21           = (0.5*rinv21*(3.0-((rsq21*rinv21)
 677      &  *rinv21)))
 678 #endif
 679
 680 C          PowerPC intrinsics 1/sqrt lookup table
 681 #ifndef GMX_BLUEGENE
 682             rinv31           = frsqrtes(rsq31)
 683 #else
 684             rinv31           = frsqrte(dble(rsq31))
 685 #endif
 686             rinv31           = (0.5*rinv31*(3.0-((rsq31*rinv31)
 687      &  *rinv31)))
 688 #ifdef GMX_DOUBLE
 689             rinv31           = (0.5*rinv31*(3.0-((rsq31*rinv31)
 690      &  *rinv31)))
 691 #endif
 692
 693 C          PowerPC intrinsics 1/sqrt lookup table
 694 #ifndef GMX_BLUEGENE
 695             rinv41           = frsqrtes(rsq41)
 696 #else
 697             rinv41           = frsqrte(dble(rsq41))
 698 #endif
 699             rinv41           = (0.5*rinv41*(3.0-((rsq41*rinv41)
 700      &  *rinv41)))
 701 #ifdef GMX_DOUBLE
 702             rinv41           = (0.5*rinv41*(3.0-((rsq41*rinv41)
 703      &  *rinv41)))
 704 #endif
 705
 706 C          Load parameters for j atom
 707             tj               = nti+2*type(jnr)+1
 708             c6               = vdwparam(tj)
 709             c12              = vdwparam(tj+1)
 710
 711 C          Calculate table index
 712             r                = rsq11*rinv11
 713
 714 C          Calculate table index
 715             rt               = r*tabscale
 716             n0               = rt
 717             eps              = rt-n0
 718             eps2             = eps*eps
 719             nnn              = 12*n0+1
 720
 721 C          Tabulated VdW interaction - dispersion
 722             nnn              = nnn+4
 723             Y                = VFtab(nnn)
 724             F                = VFtab(nnn+1)
 725             Geps             = eps*VFtab(nnn+2)
 726             Heps2            = eps2*VFtab(nnn+3)
 727             Fp               = F+Geps+Heps2
 728             VV               = Y+eps*Fp
 729             Vvdw6            = c6*VV
 730
 731 C          Tabulated VdW interaction - repulsion
 732             nnn              = nnn+4
 733             Y                = VFtab(nnn)
 734             F                = VFtab(nnn+1)
 735             Geps             = eps*VFtab(nnn+2)
 736             Heps2            = eps2*VFtab(nnn+3)
 737             Fp               = F+Geps+Heps2
 738             VV               = Y+eps*Fp
 739             Vvdw12           = c12*VV
 740             Vvdwtot          = Vvdwtot+ Vvdw6 + Vvdw12
 741
 742 C          Load parameters for j atom
 743             jq               = charge(jnr+0)
 744             qq               = qH*jq
 745
 746 C          Calculate table index
 747             r                = rsq21*rinv21
 748
 749 C          Calculate table index
 750             rt               = r*tabscale
 751             n0               = rt
 752             eps              = rt-n0
 753             eps2             = eps*eps
 754             nnn              = 12*n0+1
 755
 756 C          Tabulated coulomb interaction
 757             Y                = VFtab(nnn)
 758             F                = VFtab(nnn+1)
 759             Geps             = eps*VFtab(nnn+2)
 760             Heps2            = eps2*VFtab(nnn+3)
 761             Fp               = F+Geps+Heps2
 762             VV               = Y+eps*Fp
 763             vcoul            = qq*VV
 764             vctot            = vctot + vcoul
 765
 766 C          Load parameters for j atom
 767
 768 C          Calculate table index
 769             r                = rsq31*rinv31
 770
 771 C          Calculate table index
 772             rt               = r*tabscale
 773             n0               = rt
 774             eps              = rt-n0
 775             eps2             = eps*eps
 776             nnn              = 12*n0+1
 777
 778 C          Tabulated coulomb interaction
 779             Y                = VFtab(nnn)
 780             F                = VFtab(nnn+1)
 781             Geps             = eps*VFtab(nnn+2)
 782             Heps2            = eps2*VFtab(nnn+3)
 783             Fp               = F+Geps+Heps2
 784             VV               = Y+eps*Fp
 785             vcoul            = qq*VV
 786             vctot            = vctot + vcoul
 787
 788 C          Load parameters for j atom
 789             qq               = qM*jq
 790
 791 C          Calculate table index
 792             r                = rsq41*rinv41
 793
 794 C          Calculate table index
 795             rt               = r*tabscale
 796             n0               = rt
 797             eps              = rt-n0
 798             eps2             = eps*eps
 799             nnn              = 12*n0+1
 800
 801 C          Tabulated coulomb interaction
 802             Y                = VFtab(nnn)
 803             F                = VFtab(nnn+1)
 804             Geps             = eps*VFtab(nnn+2)
 805             Heps2            = eps2*VFtab(nnn+3)
 806             Fp               = F+Geps+Heps2
 807             VV               = Y+eps*Fp
 808             vcoul            = qq*VV
 809             vctot            = vctot + vcoul
 810
 811 C          Inner loop uses 114 flops/iteration
 812           end do
 813
 814
 815 C        Add i forces to mem and shifted force list
 816
 817 C        Add potential energies to the group for this list
 818           ggid             = gid(n)+1
 819           Vc(ggid)         = Vc(ggid) + vctot
 820           Vvdw(ggid)       = Vvdw(ggid) + Vvdwtot
 821
 822 C        Increment number of inner iterations
 823           ninner           = ninner + nj1 - nj0
 824
 825 C        Outer loop uses 14 flops/iteration
 826         end do
 827
 828
 829 C      Increment number of outer iterations
 830         nouter           = nouter + nn1 - nn0
 831       if(nn1.lt.nri) goto 10
 832
 833 C    Write outer/inner iteration count to pointers
 834       outeriter        = nouter
 835       inneriter        = ninner
 836       return
 837       end
 838
 839
 840