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