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