src/gmxlib/nonbonded/nb_kernel_power6/pwr6kernel324.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 pwr6kernel324
  45 C Coulomb interaction:     Tabulated
  46 C VdW interaction:         Buckingham
  47 C water optimization:      pairs of TIP4P interactions
  48 C Calculate forces:        yes
  49 C
  50       subroutine pwr6kernel324(
  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       rinvsq
  97       gmxreal       qq,vcoul,vctot
  98       integer*4     tj
  99       gmxreal       rinvsix
 100       gmxreal       Vvdw6,Vvdwtot
 101       gmxreal       r,rt,eps,eps2
 102       integer*4     n0,nnn
 103       gmxreal       Y,F,Geps,Heps2,Fp,VV
 104       gmxreal       FF
 105       gmxreal       fijC
 106       gmxreal       Vvdwexp,br
 107       gmxreal       ix1,iy1,iz1,fix1,fiy1,fiz1
 108       gmxreal       ix2,iy2,iz2,fix2,fiy2,fiz2
 109       gmxreal       ix3,iy3,iz3,fix3,fiy3,fiz3
 110       gmxreal       ix4,iy4,iz4,fix4,fiy4,fiz4
 111       gmxreal       jx1,jy1,jz1
 112       gmxreal       jx2,jy2,jz2,fjx2,fjy2,fjz2
 113       gmxreal       jx3,jy3,jz3,fjx3,fjy3,fjz3
 114       gmxreal       jx4,jy4,jz4,fjx4,fjy4,fjz4
 115       gmxreal       dx11,dy11,dz11,rsq11,rinv11
 116       gmxreal       dx22,dy22,dz22,rsq22,rinv22
 117       gmxreal       dx23,dy23,dz23,rsq23,rinv23
 118       gmxreal       dx24,dy24,dz24,rsq24,rinv24
 119       gmxreal       dx32,dy32,dz32,rsq32,rinv32
 120       gmxreal       dx33,dy33,dz33,rsq33,rinv33
 121       gmxreal       dx34,dy34,dz34,rsq34,rinv34
 122       gmxreal       dx42,dy42,dz42,rsq42,rinv42
 123       gmxreal       dx43,dy43,dz43,rsq43,rinv43
 124       gmxreal       dx44,dy44,dz44,rsq44,rinv44
 125       gmxreal       qH,qM,qqMM,qqMH,qqHH
 126       gmxreal       c6,cexp1,cexp2
 127
 128
 129 C    Initialize water data
 130       ii               = iinr(1)+1
 131       qH               = charge(ii+1)
 132       qM               = charge(ii+3)
 133       qqMM             = facel*qM*qM
 134       qqMH             = facel*qM*qH
 135       qqHH             = facel*qH*qH
 136       tj               = 3*(ntype+1)*type(ii)+1
 137       c6               = vdwparam(tj)
 138       cexp1            = vdwparam(tj+1)
 139       cexp2            = vdwparam(tj+2)
 140
 141
 142 C    Reset outer and inner iteration counters
 143       nouter           = 0
 144       ninner           = 0
 145
 146 C    Loop over thread workunits
 147    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 148         if(nn1.gt.nri) nn1=nri
 149
 150 C      Start outer loop over neighborlists
 151
 152         do n=nn0+1,nn1
 153
 154 C        Load shift vector for this list
 155           is3              = 3*shift(n)+1
 156           shX              = shiftvec(is3)
 157           shY              = shiftvec(is3+1)
 158           shZ              = shiftvec(is3+2)
 159
 160 C        Load limits for loop over neighbors
 161           nj0              = jindex(n)+1
 162           nj1              = jindex(n+1)
 163
 164 C        Get outer coordinate index
 165           ii               = iinr(n)+1
 166           ii3              = 3*ii-2
 167
 168 C        Load i atom data, add shift vector
 169           ix1              = shX + pos(ii3+0)
 170           iy1              = shY + pos(ii3+1)
 171           iz1              = shZ + pos(ii3+2)
 172           ix2              = shX + pos(ii3+3)
 173           iy2              = shY + pos(ii3+4)
 174           iz2              = shZ + pos(ii3+5)
 175           ix3              = shX + pos(ii3+6)
 176           iy3              = shY + pos(ii3+7)
 177           iz3              = shZ + pos(ii3+8)
 178           ix4              = shX + pos(ii3+9)
 179           iy4              = shY + pos(ii3+10)
 180           iz4              = shZ + pos(ii3+11)
 181
 182 C        Zero the potential energy for this list
 183           vctot            = 0
 184           Vvdwtot          = 0
 185
 186 C        Clear i atom forces
 187           fix1             = 0
 188           fiy1             = 0
 189           fiz1             = 0
 190           fix2             = 0
 191           fiy2             = 0
 192           fiz2             = 0
 193           fix3             = 0
 194           fiy3             = 0
 195           fiz3             = 0
 196           fix4             = 0
 197           fiy4             = 0
 198           fiz4             = 0
 199
 200           do k=nj0,nj1
 201
 202 C          Get j neighbor index, and coordinate index
 203             jnr              = jjnr(k)+1
 204             j3               = 3*jnr-2
 205
 206 C          load j atom coordinates
 207             jx1              = pos(j3+0)
 208             jy1              = pos(j3+1)
 209             jz1              = pos(j3+2)
 210             jx2              = pos(j3+3)
 211             jy2              = pos(j3+4)
 212             jz2              = pos(j3+5)
 213             jx3              = pos(j3+6)
 214             jy3              = pos(j3+7)
 215             jz3              = pos(j3+8)
 216             jx4              = pos(j3+9)
 217             jy4              = pos(j3+10)
 218             jz4              = pos(j3+11)
 219
 220 C          Calculate distance
 221             dx11             = ix1 - jx1
 222             dy11             = iy1 - jy1
 223             dz11             = iz1 - jz1
 224             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11
 225             dx22             = ix2 - jx2
 226             dy22             = iy2 - jy2
 227             dz22             = iz2 - jz2
 228             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22
 229             dx23             = ix2 - jx3
 230             dy23             = iy2 - jy3
 231             dz23             = iz2 - jz3
 232             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23
 233             dx24             = ix2 - jx4
 234             dy24             = iy2 - jy4
 235             dz24             = iz2 - jz4
 236             rsq24            = dx24*dx24+dy24*dy24+dz24*dz24
 237             dx32             = ix3 - jx2
 238             dy32             = iy3 - jy2
 239             dz32             = iz3 - jz2
 240             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32
 241             dx33             = ix3 - jx3
 242             dy33             = iy3 - jy3
 243             dz33             = iz3 - jz3
 244             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33
 245             dx34             = ix3 - jx4
 246             dy34             = iy3 - jy4
 247             dz34             = iz3 - jz4
 248             rsq34            = dx34*dx34+dy34*dy34+dz34*dz34
 249             dx42             = ix4 - jx2
 250             dy42             = iy4 - jy2
 251             dz42             = iz4 - jz2
 252             rsq42            = dx42*dx42+dy42*dy42+dz42*dz42
 253             dx43             = ix4 - jx3
 254             dy43             = iy4 - jy3
 255             dz43             = iz4 - jz3
 256             rsq43            = dx43*dx43+dy43*dy43+dz43*dz43
 257             dx44             = ix4 - jx4
 258             dy44             = iy4 - jy4
 259             dz44             = iz4 - jz4
 260             rsq44            = dx44*dx44+dy44*dy44+dz44*dz44
 261
 262 C          Calculate 1/r and 1/r2
 263
 264 C          PowerPC intrinsics 1/sqrt lookup table
 265 #ifndef GMX_BLUEGENE
 266             rinv11           = frsqrtes(rsq11)
 267 #else
 268             rinv11           = frsqrte(dble(rsq11))
 269 #endif
 270             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 271      &  *rinv11)))
 272 #ifdef GMX_DOUBLE
 273             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
 274      &  *rinv11)))
 275 #endif
 276
 277 C          PowerPC intrinsics 1/sqrt lookup table
 278 #ifndef GMX_BLUEGENE
 279             rinv22           = frsqrtes(rsq22)
 280 #else
 281             rinv22           = frsqrte(dble(rsq22))
 282 #endif
 283             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 284      &  *rinv22)))
 285 #ifdef GMX_DOUBLE
 286             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 287      &  *rinv22)))
 288 #endif
 289
 290 C          PowerPC intrinsics 1/sqrt lookup table
 291 #ifndef GMX_BLUEGENE
 292             rinv23           = frsqrtes(rsq23)
 293 #else
 294             rinv23           = frsqrte(dble(rsq23))
 295 #endif
 296             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 297      &  *rinv23)))
 298 #ifdef GMX_DOUBLE
 299             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 300      &  *rinv23)))
 301 #endif
 302
 303 C          PowerPC intrinsics 1/sqrt lookup table
 304 #ifndef GMX_BLUEGENE
 305             rinv24           = frsqrtes(rsq24)
 306 #else
 307             rinv24           = frsqrte(dble(rsq24))
 308 #endif
 309             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 310      &  *rinv24)))
 311 #ifdef GMX_DOUBLE
 312             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 313      &  *rinv24)))
 314 #endif
 315
 316 C          PowerPC intrinsics 1/sqrt lookup table
 317 #ifndef GMX_BLUEGENE
 318             rinv32           = frsqrtes(rsq32)
 319 #else
 320             rinv32           = frsqrte(dble(rsq32))
 321 #endif
 322             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 323      &  *rinv32)))
 324 #ifdef GMX_DOUBLE
 325             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 326      &  *rinv32)))
 327 #endif
 328
 329 C          PowerPC intrinsics 1/sqrt lookup table
 330 #ifndef GMX_BLUEGENE
 331             rinv33           = frsqrtes(rsq33)
 332 #else
 333             rinv33           = frsqrte(dble(rsq33))
 334 #endif
 335             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 336      &  *rinv33)))
 337 #ifdef GMX_DOUBLE
 338             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 339      &  *rinv33)))
 340 #endif
 341
 342 C          PowerPC intrinsics 1/sqrt lookup table
 343 #ifndef GMX_BLUEGENE
 344             rinv34           = frsqrtes(rsq34)
 345 #else
 346             rinv34           = frsqrte(dble(rsq34))
 347 #endif
 348             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 349      &  *rinv34)))
 350 #ifdef GMX_DOUBLE
 351             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 352      &  *rinv34)))
 353 #endif
 354
 355 C          PowerPC intrinsics 1/sqrt lookup table
 356 #ifndef GMX_BLUEGENE
 357             rinv42           = frsqrtes(rsq42)
 358 #else
 359             rinv42           = frsqrte(dble(rsq42))
 360 #endif
 361             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 362      &  *rinv42)))
 363 #ifdef GMX_DOUBLE
 364             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 365      &  *rinv42)))
 366 #endif
 367
 368 C          PowerPC intrinsics 1/sqrt lookup table
 369 #ifndef GMX_BLUEGENE
 370             rinv43           = frsqrtes(rsq43)
 371 #else
 372             rinv43           = frsqrte(dble(rsq43))
 373 #endif
 374             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 375      &  *rinv43)))
 376 #ifdef GMX_DOUBLE
 377             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 378      &  *rinv43)))
 379 #endif
 380
 381 C          PowerPC intrinsics 1/sqrt lookup table
 382 #ifndef GMX_BLUEGENE
 383             rinv44           = frsqrtes(rsq44)
 384 #else
 385             rinv44           = frsqrte(dble(rsq44))
 386 #endif
 387             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 388      &  *rinv44)))
 389 #ifdef GMX_DOUBLE
 390             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 391      &  *rinv44)))
 392 #endif
 393
 394 C          Load parameters for j atom
 395             rinvsq           = rinv11*rinv11
 396
 397 C          Buckingham interaction
 398             rinvsix          = rinvsq*rinvsq*rinvsq
 399             Vvdw6            = c6*rinvsix
 400             br               = cexp2*rsq11*rinv11
 401             Vvdwexp          = cexp1*exp(-br)
 402             Vvdwtot          = Vvdwtot+Vvdwexp-Vvdw6
 403             fscal            = (br*Vvdwexp-6.0*Vvdw6)*rinvsq
 404
 405 C          Calculate temporary vectorial force
 406             tx               = fscal*dx11
 407             ty               = fscal*dy11
 408             tz               = fscal*dz11
 409
 410 C          Increment i atom force
 411             fix1             = fix1 + tx
 412             fiy1             = fiy1 + ty
 413             fiz1             = fiz1 + tz
 414
 415 C          Decrement j atom force
 416             faction(j3+0)    = faction(j3+0) - tx
 417             faction(j3+1)    = faction(j3+1) - ty
 418             faction(j3+2)    = faction(j3+2) - tz
 419
 420 C          Load parameters for j atom
 421             qq               = qqHH
 422
 423 C          Calculate table index
 424             r                = rsq22*rinv22
 425
 426 C          Calculate table index
 427             rt               = r*tabscale
 428             n0               = rt
 429             eps              = rt-n0
 430             eps2             = eps*eps
 431             nnn              = 4*n0+1
 432
 433 C          Tabulated coulomb interaction
 434             Y                = VFtab(nnn)
 435             F                = VFtab(nnn+1)
 436             Geps             = eps*VFtab(nnn+2)
 437             Heps2            = eps2*VFtab(nnn+3)
 438             Fp               = F+Geps+Heps2
 439             VV               = Y+eps*Fp
 440             FF               = Fp+Geps+2.0*Heps2
 441             vcoul            = qq*VV
 442             fijC             = qq*FF
 443             vctot            = vctot + vcoul
 444             fscal            = -((fijC)*tabscale)*rinv22
 445
 446 C          Calculate temporary vectorial force
 447             tx               = fscal*dx22
 448             ty               = fscal*dy22
 449             tz               = fscal*dz22
 450
 451 C          Increment i atom force
 452             fix2             = fix2 + tx
 453             fiy2             = fiy2 + ty
 454             fiz2             = fiz2 + tz
 455
 456 C          Decrement j atom force
 457             fjx2             = faction(j3+3) - tx
 458             fjy2             = faction(j3+4) - ty
 459             fjz2             = faction(j3+5) - tz
 460
 461 C          Load parameters for j atom
 462             qq               = qqHH
 463
 464 C          Calculate table index
 465             r                = rsq23*rinv23
 466
 467 C          Calculate table index
 468             rt               = r*tabscale
 469             n0               = rt
 470             eps              = rt-n0
 471             eps2             = eps*eps
 472             nnn              = 4*n0+1
 473
 474 C          Tabulated coulomb interaction
 475             Y                = VFtab(nnn)
 476             F                = VFtab(nnn+1)
 477             Geps             = eps*VFtab(nnn+2)
 478             Heps2            = eps2*VFtab(nnn+3)
 479             Fp               = F+Geps+Heps2
 480             VV               = Y+eps*Fp
 481             FF               = Fp+Geps+2.0*Heps2
 482             vcoul            = qq*VV
 483             fijC             = qq*FF
 484             vctot            = vctot + vcoul
 485             fscal            = -((fijC)*tabscale)*rinv23
 486
 487 C          Calculate temporary vectorial force
 488             tx               = fscal*dx23
 489             ty               = fscal*dy23
 490             tz               = fscal*dz23
 491
 492 C          Increment i atom force
 493             fix2             = fix2 + tx
 494             fiy2             = fiy2 + ty
 495             fiz2             = fiz2 + tz
 496
 497 C          Decrement j atom force
 498             fjx3             = faction(j3+6) - tx
 499             fjy3             = faction(j3+7) - ty
 500             fjz3             = faction(j3+8) - tz
 501
 502 C          Load parameters for j atom
 503             qq               = qqMH
 504
 505 C          Calculate table index
 506             r                = rsq24*rinv24
 507
 508 C          Calculate table index
 509             rt               = r*tabscale
 510             n0               = rt
 511             eps              = rt-n0
 512             eps2             = eps*eps
 513             nnn              = 4*n0+1
 514
 515 C          Tabulated coulomb interaction
 516             Y                = VFtab(nnn)
 517             F                = VFtab(nnn+1)
 518             Geps             = eps*VFtab(nnn+2)
 519             Heps2            = eps2*VFtab(nnn+3)
 520             Fp               = F+Geps+Heps2
 521             VV               = Y+eps*Fp
 522             FF               = Fp+Geps+2.0*Heps2
 523             vcoul            = qq*VV
 524             fijC             = qq*FF
 525             vctot            = vctot + vcoul
 526             fscal            = -((fijC)*tabscale)*rinv24
 527
 528 C          Calculate temporary vectorial force
 529             tx               = fscal*dx24
 530             ty               = fscal*dy24
 531             tz               = fscal*dz24
 532
 533 C          Increment i atom force
 534             fix2             = fix2 + tx
 535             fiy2             = fiy2 + ty
 536             fiz2             = fiz2 + tz
 537
 538 C          Decrement j atom force
 539             fjx4             = faction(j3+9) - tx
 540             fjy4             = faction(j3+10) - ty
 541             fjz4             = faction(j3+11) - tz
 542
 543 C          Load parameters for j atom
 544             qq               = qqHH
 545
 546 C          Calculate table index
 547             r                = rsq32*rinv32
 548
 549 C          Calculate table index
 550             rt               = r*tabscale
 551             n0               = rt
 552             eps              = rt-n0
 553             eps2             = eps*eps
 554             nnn              = 4*n0+1
 555
 556 C          Tabulated coulomb interaction
 557             Y                = VFtab(nnn)
 558             F                = VFtab(nnn+1)
 559             Geps             = eps*VFtab(nnn+2)
 560             Heps2            = eps2*VFtab(nnn+3)
 561             Fp               = F+Geps+Heps2
 562             VV               = Y+eps*Fp
 563             FF               = Fp+Geps+2.0*Heps2
 564             vcoul            = qq*VV
 565             fijC             = qq*FF
 566             vctot            = vctot + vcoul
 567             fscal            = -((fijC)*tabscale)*rinv32
 568
 569 C          Calculate temporary vectorial force
 570             tx               = fscal*dx32
 571             ty               = fscal*dy32
 572             tz               = fscal*dz32
 573
 574 C          Increment i atom force
 575             fix3             = fix3 + tx
 576             fiy3             = fiy3 + ty
 577             fiz3             = fiz3 + tz
 578
 579 C          Decrement j atom force
 580             fjx2             = fjx2 - tx
 581             fjy2             = fjy2 - ty
 582             fjz2             = fjz2 - tz
 583
 584 C          Load parameters for j atom
 585             qq               = qqHH
 586
 587 C          Calculate table index
 588             r                = rsq33*rinv33
 589
 590 C          Calculate table index
 591             rt               = r*tabscale
 592             n0               = rt
 593             eps              = rt-n0
 594             eps2             = eps*eps
 595             nnn              = 4*n0+1
 596
 597 C          Tabulated coulomb interaction
 598             Y                = VFtab(nnn)
 599             F                = VFtab(nnn+1)
 600             Geps             = eps*VFtab(nnn+2)
 601             Heps2            = eps2*VFtab(nnn+3)
 602             Fp               = F+Geps+Heps2
 603             VV               = Y+eps*Fp
 604             FF               = Fp+Geps+2.0*Heps2
 605             vcoul            = qq*VV
 606             fijC             = qq*FF
 607             vctot            = vctot + vcoul
 608             fscal            = -((fijC)*tabscale)*rinv33
 609
 610 C          Calculate temporary vectorial force
 611             tx               = fscal*dx33
 612             ty               = fscal*dy33
 613             tz               = fscal*dz33
 614
 615 C          Increment i atom force
 616             fix3             = fix3 + tx
 617             fiy3             = fiy3 + ty
 618             fiz3             = fiz3 + tz
 619
 620 C          Decrement j atom force
 621             fjx3             = fjx3 - tx
 622             fjy3             = fjy3 - ty
 623             fjz3             = fjz3 - tz
 624
 625 C          Load parameters for j atom
 626             qq               = qqMH
 627
 628 C          Calculate table index
 629             r                = rsq34*rinv34
 630
 631 C          Calculate table index
 632             rt               = r*tabscale
 633             n0               = rt
 634             eps              = rt-n0
 635             eps2             = eps*eps
 636             nnn              = 4*n0+1
 637
 638 C          Tabulated coulomb interaction
 639             Y                = VFtab(nnn)
 640             F                = VFtab(nnn+1)
 641             Geps             = eps*VFtab(nnn+2)
 642             Heps2            = eps2*VFtab(nnn+3)
 643             Fp               = F+Geps+Heps2
 644             VV               = Y+eps*Fp
 645             FF               = Fp+Geps+2.0*Heps2
 646             vcoul            = qq*VV
 647             fijC             = qq*FF
 648             vctot            = vctot + vcoul
 649             fscal            = -((fijC)*tabscale)*rinv34
 650
 651 C          Calculate temporary vectorial force
 652             tx               = fscal*dx34
 653             ty               = fscal*dy34
 654             tz               = fscal*dz34
 655
 656 C          Increment i atom force
 657             fix3             = fix3 + tx
 658             fiy3             = fiy3 + ty
 659             fiz3             = fiz3 + tz
 660
 661 C          Decrement j atom force
 662             fjx4             = fjx4 - tx
 663             fjy4             = fjy4 - ty
 664             fjz4             = fjz4 - tz
 665
 666 C          Load parameters for j atom
 667             qq               = qqMH
 668
 669 C          Calculate table index
 670             r                = rsq42*rinv42
 671
 672 C          Calculate table index
 673             rt               = r*tabscale
 674             n0               = rt
 675             eps              = rt-n0
 676             eps2             = eps*eps
 677             nnn              = 4*n0+1
 678
 679 C          Tabulated coulomb interaction
 680             Y                = VFtab(nnn)
 681             F                = VFtab(nnn+1)
 682             Geps             = eps*VFtab(nnn+2)
 683             Heps2            = eps2*VFtab(nnn+3)
 684             Fp               = F+Geps+Heps2
 685             VV               = Y+eps*Fp
 686             FF               = Fp+Geps+2.0*Heps2
 687             vcoul            = qq*VV
 688             fijC             = qq*FF
 689             vctot            = vctot + vcoul
 690             fscal            = -((fijC)*tabscale)*rinv42
 691
 692 C          Calculate temporary vectorial force
 693             tx               = fscal*dx42
 694             ty               = fscal*dy42
 695             tz               = fscal*dz42
 696
 697 C          Increment i atom force
 698             fix4             = fix4 + tx
 699             fiy4             = fiy4 + ty
 700             fiz4             = fiz4 + tz
 701
 702 C          Decrement j atom force
 703             faction(j3+3)    = fjx2 - tx
 704             faction(j3+4)    = fjy2 - ty
 705             faction(j3+5)    = fjz2 - tz
 706
 707 C          Load parameters for j atom
 708             qq               = qqMH
 709
 710 C          Calculate table index
 711             r                = rsq43*rinv43
 712
 713 C          Calculate table index
 714             rt               = r*tabscale
 715             n0               = rt
 716             eps              = rt-n0
 717             eps2             = eps*eps
 718             nnn              = 4*n0+1
 719
 720 C          Tabulated coulomb interaction
 721             Y                = VFtab(nnn)
 722             F                = VFtab(nnn+1)
 723             Geps             = eps*VFtab(nnn+2)
 724             Heps2            = eps2*VFtab(nnn+3)
 725             Fp               = F+Geps+Heps2
 726             VV               = Y+eps*Fp
 727             FF               = Fp+Geps+2.0*Heps2
 728             vcoul            = qq*VV
 729             fijC             = qq*FF
 730             vctot            = vctot + vcoul
 731             fscal            = -((fijC)*tabscale)*rinv43
 732
 733 C          Calculate temporary vectorial force
 734             tx               = fscal*dx43
 735             ty               = fscal*dy43
 736             tz               = fscal*dz43
 737
 738 C          Increment i atom force
 739             fix4             = fix4 + tx
 740             fiy4             = fiy4 + ty
 741             fiz4             = fiz4 + tz
 742
 743 C          Decrement j atom force
 744             faction(j3+6)    = fjx3 - tx
 745             faction(j3+7)    = fjy3 - ty
 746             faction(j3+8)    = fjz3 - tz
 747
 748 C          Load parameters for j atom
 749             qq               = qqMM
 750
 751 C          Calculate table index
 752             r                = rsq44*rinv44
 753
 754 C          Calculate table index
 755             rt               = r*tabscale
 756             n0               = rt
 757             eps              = rt-n0
 758             eps2             = eps*eps
 759             nnn              = 4*n0+1
 760
 761 C          Tabulated coulomb interaction
 762             Y                = VFtab(nnn)
 763             F                = VFtab(nnn+1)
 764             Geps             = eps*VFtab(nnn+2)
 765             Heps2            = eps2*VFtab(nnn+3)
 766             Fp               = F+Geps+Heps2
 767             VV               = Y+eps*Fp
 768             FF               = Fp+Geps+2.0*Heps2
 769             vcoul            = qq*VV
 770             fijC             = qq*FF
 771             vctot            = vctot + vcoul
 772             fscal            = -((fijC)*tabscale)*rinv44
 773
 774 C          Calculate temporary vectorial force
 775             tx               = fscal*dx44
 776             ty               = fscal*dy44
 777             tz               = fscal*dz44
 778
 779 C          Increment i atom force
 780             fix4             = fix4 + tx
 781             fiy4             = fiy4 + ty
 782             fiz4             = fiz4 + tz
 783
 784 C          Decrement j atom force
 785             faction(j3+9)    = fjx4 - tx
 786             faction(j3+10)   = fjy4 - ty
 787             faction(j3+11)   = fjz4 - tz
 788
 789 C          Inner loop uses 440 flops/iteration
 790           end do
 791
 792
 793 C        Add i forces to mem and shifted force list
 794           faction(ii3+0)   = faction(ii3+0) + fix1
 795           faction(ii3+1)   = faction(ii3+1) + fiy1
 796           faction(ii3+2)   = faction(ii3+2) + fiz1
 797           faction(ii3+3)   = faction(ii3+3) + fix2
 798           faction(ii3+4)   = faction(ii3+4) + fiy2
 799           faction(ii3+5)   = faction(ii3+5) + fiz2
 800           faction(ii3+6)   = faction(ii3+6) + fix3
 801           faction(ii3+7)   = faction(ii3+7) + fiy3
 802           faction(ii3+8)   = faction(ii3+8) + fiz3
 803           faction(ii3+9)   = faction(ii3+9) + fix4
 804           faction(ii3+10)  = faction(ii3+10) + fiy4
 805           faction(ii3+11)  = faction(ii3+11) + fiz4
 806           fshift(is3)      = fshift(is3)+fix1+fix2+fix3+fix4
 807           fshift(is3+1)    = fshift(is3+1)+fiy1+fiy2+fiy3+fiy4
 808           fshift(is3+2)    = fshift(is3+2)+fiz1+fiz2+fiz3+fiz4
 809
 810 C        Add potential energies to the group for this list
 811           ggid             = gid(n)+1
 812           Vc(ggid)         = Vc(ggid) + vctot
 813           Vvdw(ggid)       = Vvdw(ggid) + Vvdwtot
 814
 815 C        Increment number of inner iterations
 816           ninner           = ninner + nj1 - nj0
 817
 818 C        Outer loop uses 38 flops/iteration
 819         end do
 820
 821
 822 C      Increment number of outer iterations
 823         nouter           = nouter + nn1 - nn0
 824       if(nn1.lt.nri) goto 10
 825
 826 C    Write outer/inner iteration count to pointers
 827       outeriter        = nouter
 828       inneriter        = ninner
 829       return
 830       end
 831
 832
 833
 834
 835
 836
 837 C
 838 C Gromacs nonbonded kernel pwr6kernel324nf
 839 C Coulomb interaction:     Tabulated
 840 C VdW interaction:         Buckingham
 841 C water optimization:      pairs of TIP4P interactions
 842 C Calculate forces:        no
 843 C
 844       subroutine pwr6kernel324nf(
 845      &                          nri,
 846      &                          iinr,
 847      &                          jindex,
 848      &                          jjnr,
 849      &                          shift,
 850      &                          shiftvec,
 851      &                          fshift,
 852      &                          gid,
 853      &                          pos,
 854      &                          faction,
 855      &                          charge,
 856      &                          facel,
 857      &                          krf,
 858      &                          crf,
 859      &                          Vc,
 860      &                          type,
 861      &                          ntype,
 862      &                          vdwparam,
 863      &                          Vvdw,
 864      &                          tabscale,
 865      &                          VFtab,
 866      &                          invsqrta,
 867      &                          dvda,
 868      &                          gbtabscale,
 869      &                          GBtab,
 870      &                          nthreads,
 871      &                          count,
 872      &                          mtx,
 873      &                          outeriter,
 874      &                          inneriter,
 875      &                          work)
 876       implicit      none
 877       integer*4     nri,iinr(*),jindex(*),jjnr(*),shift(*)
 878       gmxreal       shiftvec(*),fshift(*),pos(*),faction(*)
 879       integer*4     gid(*),type(*),ntype
 880       gmxreal       charge(*),facel,krf,crf,Vc(*),vdwparam(*)
 881       gmxreal       Vvdw(*),tabscale,VFtab(*)
 882       gmxreal       invsqrta(*),dvda(*),gbtabscale,GBtab(*)
 883       integer*4     nthreads,count,mtx,outeriter,inneriter
 884       gmxreal       work(*)
 885
 886       integer*4     n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
 887       integer*4     nn0,nn1,nouter,ninner
 888       gmxreal       shX,shY,shZ
 889       gmxreal       rinvsq
 890       gmxreal       qq,vcoul,vctot
 891       integer*4     tj
 892       gmxreal       rinvsix
 893       gmxreal       Vvdw6,Vvdwtot
 894       gmxreal       r,rt,eps,eps2
 895       integer*4     n0,nnn
 896       gmxreal       Y,F,Geps,Heps2,Fp,VV
 897       gmxreal       Vvdwexp,br
 898       gmxreal       ix1,iy1,iz1
 899       gmxreal       ix2,iy2,iz2
 900       gmxreal       ix3,iy3,iz3
 901       gmxreal       ix4,iy4,iz4
 902       gmxreal       jx1,jy1,jz1
 903       gmxreal       jx2,jy2,jz2
 904       gmxreal       jx3,jy3,jz3
 905       gmxreal       jx4,jy4,jz4
 906       gmxreal       dx11,dy11,dz11,rsq11,rinv11
 907       gmxreal       dx22,dy22,dz22,rsq22,rinv22
 908       gmxreal       dx23,dy23,dz23,rsq23,rinv23
 909       gmxreal       dx24,dy24,dz24,rsq24,rinv24
 910       gmxreal       dx32,dy32,dz32,rsq32,rinv32
 911       gmxreal       dx33,dy33,dz33,rsq33,rinv33
 912       gmxreal       dx34,dy34,dz34,rsq34,rinv34
 913       gmxreal       dx42,dy42,dz42,rsq42,rinv42
 914       gmxreal       dx43,dy43,dz43,rsq43,rinv43
 915       gmxreal       dx44,dy44,dz44,rsq44,rinv44
 916       gmxreal       qH,qM,qqMM,qqMH,qqHH
 917       gmxreal       c6,cexp1,cexp2
 918
 919
 920 C    Initialize water data
 921       ii               = iinr(1)+1
 922       qH               = charge(ii+1)
 923       qM               = charge(ii+3)
 924       qqMM             = facel*qM*qM
 925       qqMH             = facel*qM*qH
 926       qqHH             = facel*qH*qH
 927       tj               = 3*(ntype+1)*type(ii)+1
 928       c6               = vdwparam(tj)
 929       cexp1            = vdwparam(tj+1)
 930       cexp2            = vdwparam(tj+2)
 931
 932
 933 C    Reset outer and inner iteration counters
 934       nouter           = 0
 935       ninner           = 0
 936
 937 C    Loop over thread workunits
 938    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 939         if(nn1.gt.nri) nn1=nri
 940
 941 C      Start outer loop over neighborlists
 942
 943         do n=nn0+1,nn1
 944
 945 C        Load shift vector for this list
 946           is3              = 3*shift(n)+1
 947           shX              = shiftvec(is3)
 948           shY              = shiftvec(is3+1)
 949           shZ              = shiftvec(is3+2)
 950
 951 C        Load limits for loop over neighbors
 952           nj0              = jindex(n)+1
 953           nj1              = jindex(n+1)
 954
 955 C        Get outer coordinate index
 956           ii               = iinr(n)+1
 957           ii3              = 3*ii-2
 958
 959 C        Load i atom data, add shift vector
 960           ix1              = shX + pos(ii3+0)
 961           iy1              = shY + pos(ii3+1)
 962           iz1              = shZ + pos(ii3+2)
 963           ix2              = shX + pos(ii3+3)
 964           iy2              = shY + pos(ii3+4)
 965           iz2              = shZ + pos(ii3+5)
 966           ix3              = shX + pos(ii3+6)
 967           iy3              = shY + pos(ii3+7)
 968           iz3              = shZ + pos(ii3+8)
 969           ix4              = shX + pos(ii3+9)
 970           iy4              = shY + pos(ii3+10)
 971           iz4              = shZ + pos(ii3+11)
 972
 973 C        Zero the potential energy for this list
 974           vctot            = 0
 975           Vvdwtot          = 0
 976
 977 C        Clear i atom forces
 978
 979           do k=nj0,nj1
 980
 981 C          Get j neighbor index, and coordinate index
 982             jnr              = jjnr(k)+1
 983             j3               = 3*jnr-2
 984
 985 C          load j atom coordinates
 986             jx1              = pos(j3+0)
 987             jy1              = pos(j3+1)
 988             jz1              = pos(j3+2)
 989             jx2              = pos(j3+3)
 990             jy2              = pos(j3+4)
 991             jz2              = pos(j3+5)
 992             jx3              = pos(j3+6)
 993             jy3              = pos(j3+7)
 994             jz3              = pos(j3+8)
 995             jx4              = pos(j3+9)
 996             jy4              = pos(j3+10)
 997             jz4              = pos(j3+11)
 998
 999 C          Calculate distance
1000             dx11             = ix1 - jx1
1001             dy11             = iy1 - jy1
1002             dz11             = iz1 - jz1
1003             rsq11            = dx11*dx11+dy11*dy11+dz11*dz11
1004             dx22             = ix2 - jx2
1005             dy22             = iy2 - jy2
1006             dz22             = iz2 - jz2
1007             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22
1008             dx23             = ix2 - jx3
1009             dy23             = iy2 - jy3
1010             dz23             = iz2 - jz3
1011             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23
1012             dx24             = ix2 - jx4
1013             dy24             = iy2 - jy4
1014             dz24             = iz2 - jz4
1015             rsq24            = dx24*dx24+dy24*dy24+dz24*dz24
1016             dx32             = ix3 - jx2
1017             dy32             = iy3 - jy2
1018             dz32             = iz3 - jz2
1019             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32
1020             dx33             = ix3 - jx3
1021             dy33             = iy3 - jy3
1022             dz33             = iz3 - jz3
1023             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33
1024             dx34             = ix3 - jx4
1025             dy34             = iy3 - jy4
1026             dz34             = iz3 - jz4
1027             rsq34            = dx34*dx34+dy34*dy34+dz34*dz34
1028             dx42             = ix4 - jx2
1029             dy42             = iy4 - jy2
1030             dz42             = iz4 - jz2
1031             rsq42            = dx42*dx42+dy42*dy42+dz42*dz42
1032             dx43             = ix4 - jx3
1033             dy43             = iy4 - jy3
1034             dz43             = iz4 - jz3
1035             rsq43            = dx43*dx43+dy43*dy43+dz43*dz43
1036             dx44             = ix4 - jx4
1037             dy44             = iy4 - jy4
1038             dz44             = iz4 - jz4
1039             rsq44            = dx44*dx44+dy44*dy44+dz44*dz44
1040
1041 C          Calculate 1/r and 1/r2
1042
1043 C          PowerPC intrinsics 1/sqrt lookup table
1044 #ifndef GMX_BLUEGENE
1045             rinv11           = frsqrtes(rsq11)
1046 #else
1047             rinv11           = frsqrte(dble(rsq11))
1048 #endif
1049             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
1050      &  *rinv11)))
1051 #ifdef GMX_DOUBLE
1052             rinv11           = (0.5*rinv11*(3.0-((rsq11*rinv11)
1053      &  *rinv11)))
1054 #endif
1055
1056 C          PowerPC intrinsics 1/sqrt lookup table
1057 #ifndef GMX_BLUEGENE
1058             rinv22           = frsqrtes(rsq22)
1059 #else
1060             rinv22           = frsqrte(dble(rsq22))
1061 #endif
1062             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
1063      &  *rinv22)))
1064 #ifdef GMX_DOUBLE
1065             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
1066      &  *rinv22)))
1067 #endif
1068
1069 C          PowerPC intrinsics 1/sqrt lookup table
1070 #ifndef GMX_BLUEGENE
1071             rinv23           = frsqrtes(rsq23)
1072 #else
1073             rinv23           = frsqrte(dble(rsq23))
1074 #endif
1075             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
1076      &  *rinv23)))
1077 #ifdef GMX_DOUBLE
1078             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
1079      &  *rinv23)))
1080 #endif
1081
1082 C          PowerPC intrinsics 1/sqrt lookup table
1083 #ifndef GMX_BLUEGENE
1084             rinv24           = frsqrtes(rsq24)
1085 #else
1086             rinv24           = frsqrte(dble(rsq24))
1087 #endif
1088             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
1089      &  *rinv24)))
1090 #ifdef GMX_DOUBLE
1091             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
1092      &  *rinv24)))
1093 #endif
1094
1095 C          PowerPC intrinsics 1/sqrt lookup table
1096 #ifndef GMX_BLUEGENE
1097             rinv32           = frsqrtes(rsq32)
1098 #else
1099             rinv32           = frsqrte(dble(rsq32))
1100 #endif
1101             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
1102      &  *rinv32)))
1103 #ifdef GMX_DOUBLE
1104             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
1105      &  *rinv32)))
1106 #endif
1107
1108 C          PowerPC intrinsics 1/sqrt lookup table
1109 #ifndef GMX_BLUEGENE
1110             rinv33           = frsqrtes(rsq33)
1111 #else
1112             rinv33           = frsqrte(dble(rsq33))
1113 #endif
1114             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
1115      &  *rinv33)))
1116 #ifdef GMX_DOUBLE
1117             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
1118      &  *rinv33)))
1119 #endif
1120
1121 C          PowerPC intrinsics 1/sqrt lookup table
1122 #ifndef GMX_BLUEGENE
1123             rinv34           = frsqrtes(rsq34)
1124 #else
1125             rinv34           = frsqrte(dble(rsq34))
1126 #endif
1127             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
1128      &  *rinv34)))
1129 #ifdef GMX_DOUBLE
1130             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
1131      &  *rinv34)))
1132 #endif
1133
1134 C          PowerPC intrinsics 1/sqrt lookup table
1135 #ifndef GMX_BLUEGENE
1136             rinv42           = frsqrtes(rsq42)
1137 #else
1138             rinv42           = frsqrte(dble(rsq42))
1139 #endif
1140             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
1141      &  *rinv42)))
1142 #ifdef GMX_DOUBLE
1143             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
1144      &  *rinv42)))
1145 #endif
1146
1147 C          PowerPC intrinsics 1/sqrt lookup table
1148 #ifndef GMX_BLUEGENE
1149             rinv43           = frsqrtes(rsq43)
1150 #else
1151             rinv43           = frsqrte(dble(rsq43))
1152 #endif
1153             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
1154      &  *rinv43)))
1155 #ifdef GMX_DOUBLE
1156             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
1157      &  *rinv43)))
1158 #endif
1159
1160 C          PowerPC intrinsics 1/sqrt lookup table
1161 #ifndef GMX_BLUEGENE
1162             rinv44           = frsqrtes(rsq44)
1163 #else
1164             rinv44           = frsqrte(dble(rsq44))
1165 #endif
1166             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
1167      &  *rinv44)))
1168 #ifdef GMX_DOUBLE
1169             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
1170      &  *rinv44)))
1171 #endif
1172
1173 C          Load parameters for j atom
1174             rinvsq           = rinv11*rinv11
1175
1176 C          Buckingham interaction
1177             rinvsix          = rinvsq*rinvsq*rinvsq
1178             Vvdw6            = c6*rinvsix
1179             br               = cexp2*rsq11*rinv11
1180             Vvdwexp          = cexp1*exp(-br)
1181             Vvdwtot          = Vvdwtot+Vvdwexp-Vvdw6
1182
1183 C          Load parameters for j atom
1184             qq               = qqHH
1185
1186 C          Calculate table index
1187             r                = rsq22*rinv22
1188
1189 C          Calculate table index
1190             rt               = r*tabscale
1191             n0               = rt
1192             eps              = rt-n0
1193             eps2             = eps*eps
1194             nnn              = 4*n0+1
1195
1196 C          Tabulated coulomb interaction
1197             Y                = VFtab(nnn)
1198             F                = VFtab(nnn+1)
1199             Geps             = eps*VFtab(nnn+2)
1200             Heps2            = eps2*VFtab(nnn+3)
1201             Fp               = F+Geps+Heps2
1202             VV               = Y+eps*Fp
1203             vcoul            = qq*VV
1204             vctot            = vctot + vcoul
1205
1206 C          Load parameters for j atom
1207             qq               = qqHH
1208
1209 C          Calculate table index
1210             r                = rsq23*rinv23
1211
1212 C          Calculate table index
1213             rt               = r*tabscale
1214             n0               = rt
1215             eps              = rt-n0
1216             eps2             = eps*eps
1217             nnn              = 4*n0+1
1218
1219 C          Tabulated coulomb interaction
1220             Y                = VFtab(nnn)
1221             F                = VFtab(nnn+1)
1222             Geps             = eps*VFtab(nnn+2)
1223             Heps2            = eps2*VFtab(nnn+3)
1224             Fp               = F+Geps+Heps2
1225             VV               = Y+eps*Fp
1226             vcoul            = qq*VV
1227             vctot            = vctot + vcoul
1228
1229 C          Load parameters for j atom
1230             qq               = qqMH
1231
1232 C          Calculate table index
1233             r                = rsq24*rinv24
1234
1235 C          Calculate table index
1236             rt               = r*tabscale
1237             n0               = rt
1238             eps              = rt-n0
1239             eps2             = eps*eps
1240             nnn              = 4*n0+1
1241
1242 C          Tabulated coulomb interaction
1243             Y                = VFtab(nnn)
1244             F                = VFtab(nnn+1)
1245             Geps             = eps*VFtab(nnn+2)
1246             Heps2            = eps2*VFtab(nnn+3)
1247             Fp               = F+Geps+Heps2
1248             VV               = Y+eps*Fp
1249             vcoul            = qq*VV
1250             vctot            = vctot + vcoul
1251
1252 C          Load parameters for j atom
1253             qq               = qqHH
1254
1255 C          Calculate table index
1256             r                = rsq32*rinv32
1257
1258 C          Calculate table index
1259             rt               = r*tabscale
1260             n0               = rt
1261             eps              = rt-n0
1262             eps2             = eps*eps
1263             nnn              = 4*n0+1
1264
1265 C          Tabulated coulomb interaction
1266             Y                = VFtab(nnn)
1267             F                = VFtab(nnn+1)
1268             Geps             = eps*VFtab(nnn+2)
1269             Heps2            = eps2*VFtab(nnn+3)
1270             Fp               = F+Geps+Heps2
1271             VV               = Y+eps*Fp
1272             vcoul            = qq*VV
1273             vctot            = vctot + vcoul
1274
1275 C          Load parameters for j atom
1276             qq               = qqHH
1277
1278 C          Calculate table index
1279             r                = rsq33*rinv33
1280
1281 C          Calculate table index
1282             rt               = r*tabscale
1283             n0               = rt
1284             eps              = rt-n0
1285             eps2             = eps*eps
1286             nnn              = 4*n0+1
1287
1288 C          Tabulated coulomb interaction
1289             Y                = VFtab(nnn)
1290             F                = VFtab(nnn+1)
1291             Geps             = eps*VFtab(nnn+2)
1292             Heps2            = eps2*VFtab(nnn+3)
1293             Fp               = F+Geps+Heps2
1294             VV               = Y+eps*Fp
1295             vcoul            = qq*VV
1296             vctot            = vctot + vcoul
1297
1298 C          Load parameters for j atom
1299             qq               = qqMH
1300
1301 C          Calculate table index
1302             r                = rsq34*rinv34
1303
1304 C          Calculate table index
1305             rt               = r*tabscale
1306             n0               = rt
1307             eps              = rt-n0
1308             eps2             = eps*eps
1309             nnn              = 4*n0+1
1310
1311 C          Tabulated coulomb interaction
1312             Y                = VFtab(nnn)
1313             F                = VFtab(nnn+1)
1314             Geps             = eps*VFtab(nnn+2)
1315             Heps2            = eps2*VFtab(nnn+3)
1316             Fp               = F+Geps+Heps2
1317             VV               = Y+eps*Fp
1318             vcoul            = qq*VV
1319             vctot            = vctot + vcoul
1320
1321 C          Load parameters for j atom
1322             qq               = qqMH
1323
1324 C          Calculate table index
1325             r                = rsq42*rinv42
1326
1327 C          Calculate table index
1328             rt               = r*tabscale
1329             n0               = rt
1330             eps              = rt-n0
1331             eps2             = eps*eps
1332             nnn              = 4*n0+1
1333
1334 C          Tabulated coulomb interaction
1335             Y                = VFtab(nnn)
1336             F                = VFtab(nnn+1)
1337             Geps             = eps*VFtab(nnn+2)
1338             Heps2            = eps2*VFtab(nnn+3)
1339             Fp               = F+Geps+Heps2
1340             VV               = Y+eps*Fp
1341             vcoul            = qq*VV
1342             vctot            = vctot + vcoul
1343
1344 C          Load parameters for j atom
1345             qq               = qqMH
1346
1347 C          Calculate table index
1348             r                = rsq43*rinv43
1349
1350 C          Calculate table index
1351             rt               = r*tabscale
1352             n0               = rt
1353             eps              = rt-n0
1354             eps2             = eps*eps
1355             nnn              = 4*n0+1
1356
1357 C          Tabulated coulomb interaction
1358             Y                = VFtab(nnn)
1359             F                = VFtab(nnn+1)
1360             Geps             = eps*VFtab(nnn+2)
1361             Heps2            = eps2*VFtab(nnn+3)
1362             Fp               = F+Geps+Heps2
1363             VV               = Y+eps*Fp
1364             vcoul            = qq*VV
1365             vctot            = vctot + vcoul
1366
1367 C          Load parameters for j atom
1368             qq               = qqMM
1369
1370 C          Calculate table index
1371             r                = rsq44*rinv44
1372
1373 C          Calculate table index
1374             rt               = r*tabscale
1375             n0               = rt
1376             eps              = rt-n0
1377             eps2             = eps*eps
1378             nnn              = 4*n0+1
1379
1380 C          Tabulated coulomb interaction
1381             Y                = VFtab(nnn)
1382             F                = VFtab(nnn+1)
1383             Geps             = eps*VFtab(nnn+2)
1384             Heps2            = eps2*VFtab(nnn+3)
1385             Fp               = F+Geps+Heps2
1386             VV               = Y+eps*Fp
1387             vcoul            = qq*VV
1388             vctot            = vctot + vcoul
1389
1390 C          Inner loop uses 283 flops/iteration
1391           end do
1392
1393
1394 C        Add i forces to mem and shifted force list
1395
1396 C        Add potential energies to the group for this list
1397           ggid             = gid(n)+1
1398           Vc(ggid)         = Vc(ggid) + vctot
1399           Vvdw(ggid)       = Vvdw(ggid) + Vvdwtot
1400
1401 C        Increment number of inner iterations
1402           ninner           = ninner + nj1 - nj0
1403
1404 C        Outer loop uses 14 flops/iteration
1405         end do
1406
1407
1408 C      Increment number of outer iterations
1409         nouter           = nouter + nn1 - nn0
1410       if(nn1.lt.nri) goto 10
1411
1412 C    Write outer/inner iteration count to pointers
1413       outeriter        = nouter
1414       inneriter        = ninner
1415       return
1416       end
1417
1418
1419