src/gmxlib/nonbonded/nb_kernel_power6/pwr6kernel104.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 pwr6kernel104
  45 C Coulomb interaction:     Normal Coulomb
  46 C VdW interaction:         Not calculated
  47 C water optimization:      pairs of TIP4P interactions
  48 C Calculate forces:        yes
  49 C
  50       subroutine pwr6kernel104(
  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       gmxreal       ix2,iy2,iz2,fix2,fiy2,fiz2
  99       gmxreal       ix3,iy3,iz3,fix3,fiy3,fiz3
 100       gmxreal       ix4,iy4,iz4,fix4,fiy4,fiz4
 101       gmxreal       jx2,jy2,jz2,fjx2,fjy2,fjz2
 102       gmxreal       jx3,jy3,jz3,fjx3,fjy3,fjz3
 103       gmxreal       jx4,jy4,jz4,fjx4,fjy4,fjz4
 104       gmxreal       dx22,dy22,dz22,rsq22,rinv22
 105       gmxreal       dx23,dy23,dz23,rsq23,rinv23
 106       gmxreal       dx24,dy24,dz24,rsq24,rinv24
 107       gmxreal       dx32,dy32,dz32,rsq32,rinv32
 108       gmxreal       dx33,dy33,dz33,rsq33,rinv33
 109       gmxreal       dx34,dy34,dz34,rsq34,rinv34
 110       gmxreal       dx42,dy42,dz42,rsq42,rinv42
 111       gmxreal       dx43,dy43,dz43,rsq43,rinv43
 112       gmxreal       dx44,dy44,dz44,rsq44,rinv44
 113       gmxreal       qH,qM,qqMM,qqMH,qqHH
 114
 115
 116 C    Initialize water data
 117       ii               = iinr(1)+1
 118       qH               = charge(ii+1)
 119       qM               = charge(ii+3)
 120       qqMM             = facel*qM*qM
 121       qqMH             = facel*qM*qH
 122       qqHH             = facel*qH*qH
 123
 124
 125 C    Reset outer and inner iteration counters
 126       nouter           = 0
 127       ninner           = 0
 128
 129 C    Loop over thread workunits
 130    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 131         if(nn1.gt.nri) nn1=nri
 132
 133 C      Start outer loop over neighborlists
 134
 135         do n=nn0+1,nn1
 136
 137 C        Load shift vector for this list
 138           is3              = 3*shift(n)+1
 139           shX              = shiftvec(is3)
 140           shY              = shiftvec(is3+1)
 141           shZ              = shiftvec(is3+2)
 142
 143 C        Load limits for loop over neighbors
 144           nj0              = jindex(n)+1
 145           nj1              = jindex(n+1)
 146
 147 C        Get outer coordinate index
 148           ii               = iinr(n)+1
 149           ii3              = 3*ii-2
 150
 151 C        Load i atom data, add shift vector
 152           ix2              = shX + pos(ii3+3)
 153           iy2              = shY + pos(ii3+4)
 154           iz2              = shZ + pos(ii3+5)
 155           ix3              = shX + pos(ii3+6)
 156           iy3              = shY + pos(ii3+7)
 157           iz3              = shZ + pos(ii3+8)
 158           ix4              = shX + pos(ii3+9)
 159           iy4              = shY + pos(ii3+10)
 160           iz4              = shZ + pos(ii3+11)
 161
 162 C        Zero the potential energy for this list
 163           vctot            = 0
 164
 165 C        Clear i atom forces
 166           fix2             = 0
 167           fiy2             = 0
 168           fiz2             = 0
 169           fix3             = 0
 170           fiy3             = 0
 171           fiz3             = 0
 172           fix4             = 0
 173           fiy4             = 0
 174           fiz4             = 0
 175
 176           do k=nj0,nj1
 177
 178 C          Get j neighbor index, and coordinate index
 179             jnr              = jjnr(k)+1
 180             j3               = 3*jnr-2
 181
 182 C          load j atom coordinates
 183             jx2              = pos(j3+3)
 184             jy2              = pos(j3+4)
 185             jz2              = pos(j3+5)
 186             jx3              = pos(j3+6)
 187             jy3              = pos(j3+7)
 188             jz3              = pos(j3+8)
 189             jx4              = pos(j3+9)
 190             jy4              = pos(j3+10)
 191             jz4              = pos(j3+11)
 192
 193 C          Calculate distance
 194             dx22             = ix2 - jx2
 195             dy22             = iy2 - jy2
 196             dz22             = iz2 - jz2
 197             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22
 198             dx23             = ix2 - jx3
 199             dy23             = iy2 - jy3
 200             dz23             = iz2 - jz3
 201             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23
 202             dx24             = ix2 - jx4
 203             dy24             = iy2 - jy4
 204             dz24             = iz2 - jz4
 205             rsq24            = dx24*dx24+dy24*dy24+dz24*dz24
 206             dx32             = ix3 - jx2
 207             dy32             = iy3 - jy2
 208             dz32             = iz3 - jz2
 209             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32
 210             dx33             = ix3 - jx3
 211             dy33             = iy3 - jy3
 212             dz33             = iz3 - jz3
 213             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33
 214             dx34             = ix3 - jx4
 215             dy34             = iy3 - jy4
 216             dz34             = iz3 - jz4
 217             rsq34            = dx34*dx34+dy34*dy34+dz34*dz34
 218             dx42             = ix4 - jx2
 219             dy42             = iy4 - jy2
 220             dz42             = iz4 - jz2
 221             rsq42            = dx42*dx42+dy42*dy42+dz42*dz42
 222             dx43             = ix4 - jx3
 223             dy43             = iy4 - jy3
 224             dz43             = iz4 - jz3
 225             rsq43            = dx43*dx43+dy43*dy43+dz43*dz43
 226             dx44             = ix4 - jx4
 227             dy44             = iy4 - jy4
 228             dz44             = iz4 - jz4
 229             rsq44            = dx44*dx44+dy44*dy44+dz44*dz44
 230
 231 C          Calculate 1/r and 1/r2
 232
 233 C          PowerPC intrinsics 1/sqrt lookup table
 234 #ifndef GMX_BLUEGENE
 235             rinv22           = frsqrtes(rsq22)
 236 #else
 237             rinv22           = frsqrte(dble(rsq22))
 238 #endif
 239             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 240      &  *rinv22)))
 241 #ifdef GMX_DOUBLE
 242             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 243      &  *rinv22)))
 244 #endif
 245
 246 C          PowerPC intrinsics 1/sqrt lookup table
 247 #ifndef GMX_BLUEGENE
 248             rinv23           = frsqrtes(rsq23)
 249 #else
 250             rinv23           = frsqrte(dble(rsq23))
 251 #endif
 252             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 253      &  *rinv23)))
 254 #ifdef GMX_DOUBLE
 255             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 256      &  *rinv23)))
 257 #endif
 258
 259 C          PowerPC intrinsics 1/sqrt lookup table
 260 #ifndef GMX_BLUEGENE
 261             rinv24           = frsqrtes(rsq24)
 262 #else
 263             rinv24           = frsqrte(dble(rsq24))
 264 #endif
 265             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 266      &  *rinv24)))
 267 #ifdef GMX_DOUBLE
 268             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 269      &  *rinv24)))
 270 #endif
 271
 272 C          PowerPC intrinsics 1/sqrt lookup table
 273 #ifndef GMX_BLUEGENE
 274             rinv32           = frsqrtes(rsq32)
 275 #else
 276             rinv32           = frsqrte(dble(rsq32))
 277 #endif
 278             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 279      &  *rinv32)))
 280 #ifdef GMX_DOUBLE
 281             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 282      &  *rinv32)))
 283 #endif
 284
 285 C          PowerPC intrinsics 1/sqrt lookup table
 286 #ifndef GMX_BLUEGENE
 287             rinv33           = frsqrtes(rsq33)
 288 #else
 289             rinv33           = frsqrte(dble(rsq33))
 290 #endif
 291             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 292      &  *rinv33)))
 293 #ifdef GMX_DOUBLE
 294             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 295      &  *rinv33)))
 296 #endif
 297
 298 C          PowerPC intrinsics 1/sqrt lookup table
 299 #ifndef GMX_BLUEGENE
 300             rinv34           = frsqrtes(rsq34)
 301 #else
 302             rinv34           = frsqrte(dble(rsq34))
 303 #endif
 304             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 305      &  *rinv34)))
 306 #ifdef GMX_DOUBLE
 307             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 308      &  *rinv34)))
 309 #endif
 310
 311 C          PowerPC intrinsics 1/sqrt lookup table
 312 #ifndef GMX_BLUEGENE
 313             rinv42           = frsqrtes(rsq42)
 314 #else
 315             rinv42           = frsqrte(dble(rsq42))
 316 #endif
 317             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 318      &  *rinv42)))
 319 #ifdef GMX_DOUBLE
 320             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 321      &  *rinv42)))
 322 #endif
 323
 324 C          PowerPC intrinsics 1/sqrt lookup table
 325 #ifndef GMX_BLUEGENE
 326             rinv43           = frsqrtes(rsq43)
 327 #else
 328             rinv43           = frsqrte(dble(rsq43))
 329 #endif
 330             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 331      &  *rinv43)))
 332 #ifdef GMX_DOUBLE
 333             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 334      &  *rinv43)))
 335 #endif
 336
 337 C          PowerPC intrinsics 1/sqrt lookup table
 338 #ifndef GMX_BLUEGENE
 339             rinv44           = frsqrtes(rsq44)
 340 #else
 341             rinv44           = frsqrte(dble(rsq44))
 342 #endif
 343             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 344      &  *rinv44)))
 345 #ifdef GMX_DOUBLE
 346             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 347      &  *rinv44)))
 348 #endif
 349
 350 C          Load parameters for j atom
 351             qq               = qqHH
 352             rinvsq           = rinv22*rinv22
 353
 354 C          Coulomb interaction
 355             vcoul            = qq*rinv22
 356             vctot            = vctot+vcoul
 357             fscal            = (vcoul)*rinvsq
 358
 359 C          Calculate temporary vectorial force
 360             tx               = fscal*dx22
 361             ty               = fscal*dy22
 362             tz               = fscal*dz22
 363
 364 C          Increment i atom force
 365             fix2             = fix2 + tx
 366             fiy2             = fiy2 + ty
 367             fiz2             = fiz2 + tz
 368
 369 C          Decrement j atom force
 370             fjx2             = faction(j3+3) - tx
 371             fjy2             = faction(j3+4) - ty
 372             fjz2             = faction(j3+5) - tz
 373
 374 C          Load parameters for j atom
 375             qq               = qqHH
 376             rinvsq           = rinv23*rinv23
 377
 378 C          Coulomb interaction
 379             vcoul            = qq*rinv23
 380             vctot            = vctot+vcoul
 381             fscal            = (vcoul)*rinvsq
 382
 383 C          Calculate temporary vectorial force
 384             tx               = fscal*dx23
 385             ty               = fscal*dy23
 386             tz               = fscal*dz23
 387
 388 C          Increment i atom force
 389             fix2             = fix2 + tx
 390             fiy2             = fiy2 + ty
 391             fiz2             = fiz2 + tz
 392
 393 C          Decrement j atom force
 394             fjx3             = faction(j3+6) - tx
 395             fjy3             = faction(j3+7) - ty
 396             fjz3             = faction(j3+8) - tz
 397
 398 C          Load parameters for j atom
 399             qq               = qqMH
 400             rinvsq           = rinv24*rinv24
 401
 402 C          Coulomb interaction
 403             vcoul            = qq*rinv24
 404             vctot            = vctot+vcoul
 405             fscal            = (vcoul)*rinvsq
 406
 407 C          Calculate temporary vectorial force
 408             tx               = fscal*dx24
 409             ty               = fscal*dy24
 410             tz               = fscal*dz24
 411
 412 C          Increment i atom force
 413             fix2             = fix2 + tx
 414             fiy2             = fiy2 + ty
 415             fiz2             = fiz2 + tz
 416
 417 C          Decrement j atom force
 418             fjx4             = faction(j3+9) - tx
 419             fjy4             = faction(j3+10) - ty
 420             fjz4             = faction(j3+11) - tz
 421
 422 C          Load parameters for j atom
 423             qq               = qqHH
 424             rinvsq           = rinv32*rinv32
 425
 426 C          Coulomb interaction
 427             vcoul            = qq*rinv32
 428             vctot            = vctot+vcoul
 429             fscal            = (vcoul)*rinvsq
 430
 431 C          Calculate temporary vectorial force
 432             tx               = fscal*dx32
 433             ty               = fscal*dy32
 434             tz               = fscal*dz32
 435
 436 C          Increment i atom force
 437             fix3             = fix3 + tx
 438             fiy3             = fiy3 + ty
 439             fiz3             = fiz3 + tz
 440
 441 C          Decrement j atom force
 442             fjx2             = fjx2 - tx
 443             fjy2             = fjy2 - ty
 444             fjz2             = fjz2 - tz
 445
 446 C          Load parameters for j atom
 447             qq               = qqHH
 448             rinvsq           = rinv33*rinv33
 449
 450 C          Coulomb interaction
 451             vcoul            = qq*rinv33
 452             vctot            = vctot+vcoul
 453             fscal            = (vcoul)*rinvsq
 454
 455 C          Calculate temporary vectorial force
 456             tx               = fscal*dx33
 457             ty               = fscal*dy33
 458             tz               = fscal*dz33
 459
 460 C          Increment i atom force
 461             fix3             = fix3 + tx
 462             fiy3             = fiy3 + ty
 463             fiz3             = fiz3 + tz
 464
 465 C          Decrement j atom force
 466             fjx3             = fjx3 - tx
 467             fjy3             = fjy3 - ty
 468             fjz3             = fjz3 - tz
 469
 470 C          Load parameters for j atom
 471             qq               = qqMH
 472             rinvsq           = rinv34*rinv34
 473
 474 C          Coulomb interaction
 475             vcoul            = qq*rinv34
 476             vctot            = vctot+vcoul
 477             fscal            = (vcoul)*rinvsq
 478
 479 C          Calculate temporary vectorial force
 480             tx               = fscal*dx34
 481             ty               = fscal*dy34
 482             tz               = fscal*dz34
 483
 484 C          Increment i atom force
 485             fix3             = fix3 + tx
 486             fiy3             = fiy3 + ty
 487             fiz3             = fiz3 + tz
 488
 489 C          Decrement j atom force
 490             fjx4             = fjx4 - tx
 491             fjy4             = fjy4 - ty
 492             fjz4             = fjz4 - tz
 493
 494 C          Load parameters for j atom
 495             qq               = qqMH
 496             rinvsq           = rinv42*rinv42
 497
 498 C          Coulomb interaction
 499             vcoul            = qq*rinv42
 500             vctot            = vctot+vcoul
 501             fscal            = (vcoul)*rinvsq
 502
 503 C          Calculate temporary vectorial force
 504             tx               = fscal*dx42
 505             ty               = fscal*dy42
 506             tz               = fscal*dz42
 507
 508 C          Increment i atom force
 509             fix4             = fix4 + tx
 510             fiy4             = fiy4 + ty
 511             fiz4             = fiz4 + tz
 512
 513 C          Decrement j atom force
 514             faction(j3+3)    = fjx2 - tx
 515             faction(j3+4)    = fjy2 - ty
 516             faction(j3+5)    = fjz2 - tz
 517
 518 C          Load parameters for j atom
 519             qq               = qqMH
 520             rinvsq           = rinv43*rinv43
 521
 522 C          Coulomb interaction
 523             vcoul            = qq*rinv43
 524             vctot            = vctot+vcoul
 525             fscal            = (vcoul)*rinvsq
 526
 527 C          Calculate temporary vectorial force
 528             tx               = fscal*dx43
 529             ty               = fscal*dy43
 530             tz               = fscal*dz43
 531
 532 C          Increment i atom force
 533             fix4             = fix4 + tx
 534             fiy4             = fiy4 + ty
 535             fiz4             = fiz4 + tz
 536
 537 C          Decrement j atom force
 538             faction(j3+6)    = fjx3 - tx
 539             faction(j3+7)    = fjy3 - ty
 540             faction(j3+8)    = fjz3 - tz
 541
 542 C          Load parameters for j atom
 543             qq               = qqMM
 544             rinvsq           = rinv44*rinv44
 545
 546 C          Coulomb interaction
 547             vcoul            = qq*rinv44
 548             vctot            = vctot+vcoul
 549             fscal            = (vcoul)*rinvsq
 550
 551 C          Calculate temporary vectorial force
 552             tx               = fscal*dx44
 553             ty               = fscal*dy44
 554             tz               = fscal*dz44
 555
 556 C          Increment i atom force
 557             fix4             = fix4 + tx
 558             fiy4             = fiy4 + ty
 559             fiz4             = fiz4 + tz
 560
 561 C          Decrement j atom force
 562             faction(j3+9)    = fjx4 - tx
 563             faction(j3+10)   = fjy4 - ty
 564             faction(j3+11)   = fjz4 - tz
 565
 566 C          Inner loop uses 243 flops/iteration
 567           end do
 568
 569
 570 C        Add i forces to mem and shifted force list
 571           faction(ii3+3)   = faction(ii3+3) + fix2
 572           faction(ii3+4)   = faction(ii3+4) + fiy2
 573           faction(ii3+5)   = faction(ii3+5) + fiz2
 574           faction(ii3+6)   = faction(ii3+6) + fix3
 575           faction(ii3+7)   = faction(ii3+7) + fiy3
 576           faction(ii3+8)   = faction(ii3+8) + fiz3
 577           faction(ii3+9)   = faction(ii3+9) + fix4
 578           faction(ii3+10)  = faction(ii3+10) + fiy4
 579           faction(ii3+11)  = faction(ii3+11) + fiz4
 580           fshift(is3)      = fshift(is3)+fix2+fix3+fix4
 581           fshift(is3+1)    = fshift(is3+1)+fiy2+fiy3+fiy4
 582           fshift(is3+2)    = fshift(is3+2)+fiz2+fiz3+fiz4
 583
 584 C        Add potential energies to the group for this list
 585           ggid             = gid(n)+1
 586           Vc(ggid)         = Vc(ggid) + vctot
 587
 588 C        Increment number of inner iterations
 589           ninner           = ninner + nj1 - nj0
 590
 591 C        Outer loop uses 28 flops/iteration
 592         end do
 593
 594
 595 C      Increment number of outer iterations
 596         nouter           = nouter + nn1 - nn0
 597       if(nn1.lt.nri) goto 10
 598
 599 C    Write outer/inner iteration count to pointers
 600       outeriter        = nouter
 601       inneriter        = ninner
 602       return
 603       end
 604
 605
 606
 607
 608
 609
 610 C
 611 C Gromacs nonbonded kernel pwr6kernel104nf
 612 C Coulomb interaction:     Normal Coulomb
 613 C VdW interaction:         Not calculated
 614 C water optimization:      pairs of TIP4P interactions
 615 C Calculate forces:        no
 616 C
 617       subroutine pwr6kernel104nf(
 618      &                          nri,
 619      &                          iinr,
 620      &                          jindex,
 621      &                          jjnr,
 622      &                          shift,
 623      &                          shiftvec,
 624      &                          fshift,
 625      &                          gid,
 626      &                          pos,
 627      &                          faction,
 628      &                          charge,
 629      &                          facel,
 630      &                          krf,
 631      &                          crf,
 632      &                          Vc,
 633      &                          type,
 634      &                          ntype,
 635      &                          vdwparam,
 636      &                          Vvdw,
 637      &                          tabscale,
 638      &                          VFtab,
 639      &                          invsqrta,
 640      &                          dvda,
 641      &                          gbtabscale,
 642      &                          GBtab,
 643      &                          nthreads,
 644      &                          count,
 645      &                          mtx,
 646      &                          outeriter,
 647      &                          inneriter,
 648      &                          work)
 649       implicit      none
 650       integer*4     nri,iinr(*),jindex(*),jjnr(*),shift(*)
 651       gmxreal       shiftvec(*),fshift(*),pos(*),faction(*)
 652       integer*4     gid(*),type(*),ntype
 653       gmxreal       charge(*),facel,krf,crf,Vc(*),vdwparam(*)
 654       gmxreal       Vvdw(*),tabscale,VFtab(*)
 655       gmxreal       invsqrta(*),dvda(*),gbtabscale,GBtab(*)
 656       integer*4     nthreads,count,mtx,outeriter,inneriter
 657       gmxreal       work(*)
 658
 659       integer*4     n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
 660       integer*4     nn0,nn1,nouter,ninner
 661       gmxreal       shX,shY,shZ
 662       gmxreal       qq,vcoul,vctot
 663       gmxreal       ix2,iy2,iz2
 664       gmxreal       ix3,iy3,iz3
 665       gmxreal       ix4,iy4,iz4
 666       gmxreal       jx2,jy2,jz2
 667       gmxreal       jx3,jy3,jz3
 668       gmxreal       jx4,jy4,jz4
 669       gmxreal       dx22,dy22,dz22,rsq22,rinv22
 670       gmxreal       dx23,dy23,dz23,rsq23,rinv23
 671       gmxreal       dx24,dy24,dz24,rsq24,rinv24
 672       gmxreal       dx32,dy32,dz32,rsq32,rinv32
 673       gmxreal       dx33,dy33,dz33,rsq33,rinv33
 674       gmxreal       dx34,dy34,dz34,rsq34,rinv34
 675       gmxreal       dx42,dy42,dz42,rsq42,rinv42
 676       gmxreal       dx43,dy43,dz43,rsq43,rinv43
 677       gmxreal       dx44,dy44,dz44,rsq44,rinv44
 678       gmxreal       qH,qM,qqMM,qqMH,qqHH
 679
 680
 681 C    Initialize water data
 682       ii               = iinr(1)+1
 683       qH               = charge(ii+1)
 684       qM               = charge(ii+3)
 685       qqMM             = facel*qM*qM
 686       qqMH             = facel*qM*qH
 687       qqHH             = facel*qH*qH
 688
 689
 690 C    Reset outer and inner iteration counters
 691       nouter           = 0
 692       ninner           = 0
 693
 694 C    Loop over thread workunits
 695    10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
 696         if(nn1.gt.nri) nn1=nri
 697
 698 C      Start outer loop over neighborlists
 699
 700         do n=nn0+1,nn1
 701
 702 C        Load shift vector for this list
 703           is3              = 3*shift(n)+1
 704           shX              = shiftvec(is3)
 705           shY              = shiftvec(is3+1)
 706           shZ              = shiftvec(is3+2)
 707
 708 C        Load limits for loop over neighbors
 709           nj0              = jindex(n)+1
 710           nj1              = jindex(n+1)
 711
 712 C        Get outer coordinate index
 713           ii               = iinr(n)+1
 714           ii3              = 3*ii-2
 715
 716 C        Load i atom data, add shift vector
 717           ix2              = shX + pos(ii3+3)
 718           iy2              = shY + pos(ii3+4)
 719           iz2              = shZ + pos(ii3+5)
 720           ix3              = shX + pos(ii3+6)
 721           iy3              = shY + pos(ii3+7)
 722           iz3              = shZ + pos(ii3+8)
 723           ix4              = shX + pos(ii3+9)
 724           iy4              = shY + pos(ii3+10)
 725           iz4              = shZ + pos(ii3+11)
 726
 727 C        Zero the potential energy for this list
 728           vctot            = 0
 729
 730 C        Clear i atom forces
 731
 732           do k=nj0,nj1
 733
 734 C          Get j neighbor index, and coordinate index
 735             jnr              = jjnr(k)+1
 736             j3               = 3*jnr-2
 737
 738 C          load j atom coordinates
 739             jx2              = pos(j3+3)
 740             jy2              = pos(j3+4)
 741             jz2              = pos(j3+5)
 742             jx3              = pos(j3+6)
 743             jy3              = pos(j3+7)
 744             jz3              = pos(j3+8)
 745             jx4              = pos(j3+9)
 746             jy4              = pos(j3+10)
 747             jz4              = pos(j3+11)
 748
 749 C          Calculate distance
 750             dx22             = ix2 - jx2
 751             dy22             = iy2 - jy2
 752             dz22             = iz2 - jz2
 753             rsq22            = dx22*dx22+dy22*dy22+dz22*dz22
 754             dx23             = ix2 - jx3
 755             dy23             = iy2 - jy3
 756             dz23             = iz2 - jz3
 757             rsq23            = dx23*dx23+dy23*dy23+dz23*dz23
 758             dx24             = ix2 - jx4
 759             dy24             = iy2 - jy4
 760             dz24             = iz2 - jz4
 761             rsq24            = dx24*dx24+dy24*dy24+dz24*dz24
 762             dx32             = ix3 - jx2
 763             dy32             = iy3 - jy2
 764             dz32             = iz3 - jz2
 765             rsq32            = dx32*dx32+dy32*dy32+dz32*dz32
 766             dx33             = ix3 - jx3
 767             dy33             = iy3 - jy3
 768             dz33             = iz3 - jz3
 769             rsq33            = dx33*dx33+dy33*dy33+dz33*dz33
 770             dx34             = ix3 - jx4
 771             dy34             = iy3 - jy4
 772             dz34             = iz3 - jz4
 773             rsq34            = dx34*dx34+dy34*dy34+dz34*dz34
 774             dx42             = ix4 - jx2
 775             dy42             = iy4 - jy2
 776             dz42             = iz4 - jz2
 777             rsq42            = dx42*dx42+dy42*dy42+dz42*dz42
 778             dx43             = ix4 - jx3
 779             dy43             = iy4 - jy3
 780             dz43             = iz4 - jz3
 781             rsq43            = dx43*dx43+dy43*dy43+dz43*dz43
 782             dx44             = ix4 - jx4
 783             dy44             = iy4 - jy4
 784             dz44             = iz4 - jz4
 785             rsq44            = dx44*dx44+dy44*dy44+dz44*dz44
 786
 787 C          Calculate 1/r and 1/r2
 788
 789 C          PowerPC intrinsics 1/sqrt lookup table
 790 #ifndef GMX_BLUEGENE
 791             rinv22           = frsqrtes(rsq22)
 792 #else
 793             rinv22           = frsqrte(dble(rsq22))
 794 #endif
 795             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 796      &  *rinv22)))
 797 #ifdef GMX_DOUBLE
 798             rinv22           = (0.5*rinv22*(3.0-((rsq22*rinv22)
 799      &  *rinv22)))
 800 #endif
 801
 802 C          PowerPC intrinsics 1/sqrt lookup table
 803 #ifndef GMX_BLUEGENE
 804             rinv23           = frsqrtes(rsq23)
 805 #else
 806             rinv23           = frsqrte(dble(rsq23))
 807 #endif
 808             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 809      &  *rinv23)))
 810 #ifdef GMX_DOUBLE
 811             rinv23           = (0.5*rinv23*(3.0-((rsq23*rinv23)
 812      &  *rinv23)))
 813 #endif
 814
 815 C          PowerPC intrinsics 1/sqrt lookup table
 816 #ifndef GMX_BLUEGENE
 817             rinv24           = frsqrtes(rsq24)
 818 #else
 819             rinv24           = frsqrte(dble(rsq24))
 820 #endif
 821             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 822      &  *rinv24)))
 823 #ifdef GMX_DOUBLE
 824             rinv24           = (0.5*rinv24*(3.0-((rsq24*rinv24)
 825      &  *rinv24)))
 826 #endif
 827
 828 C          PowerPC intrinsics 1/sqrt lookup table
 829 #ifndef GMX_BLUEGENE
 830             rinv32           = frsqrtes(rsq32)
 831 #else
 832             rinv32           = frsqrte(dble(rsq32))
 833 #endif
 834             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 835      &  *rinv32)))
 836 #ifdef GMX_DOUBLE
 837             rinv32           = (0.5*rinv32*(3.0-((rsq32*rinv32)
 838      &  *rinv32)))
 839 #endif
 840
 841 C          PowerPC intrinsics 1/sqrt lookup table
 842 #ifndef GMX_BLUEGENE
 843             rinv33           = frsqrtes(rsq33)
 844 #else
 845             rinv33           = frsqrte(dble(rsq33))
 846 #endif
 847             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 848      &  *rinv33)))
 849 #ifdef GMX_DOUBLE
 850             rinv33           = (0.5*rinv33*(3.0-((rsq33*rinv33)
 851      &  *rinv33)))
 852 #endif
 853
 854 C          PowerPC intrinsics 1/sqrt lookup table
 855 #ifndef GMX_BLUEGENE
 856             rinv34           = frsqrtes(rsq34)
 857 #else
 858             rinv34           = frsqrte(dble(rsq34))
 859 #endif
 860             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 861      &  *rinv34)))
 862 #ifdef GMX_DOUBLE
 863             rinv34           = (0.5*rinv34*(3.0-((rsq34*rinv34)
 864      &  *rinv34)))
 865 #endif
 866
 867 C          PowerPC intrinsics 1/sqrt lookup table
 868 #ifndef GMX_BLUEGENE
 869             rinv42           = frsqrtes(rsq42)
 870 #else
 871             rinv42           = frsqrte(dble(rsq42))
 872 #endif
 873             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 874      &  *rinv42)))
 875 #ifdef GMX_DOUBLE
 876             rinv42           = (0.5*rinv42*(3.0-((rsq42*rinv42)
 877      &  *rinv42)))
 878 #endif
 879
 880 C          PowerPC intrinsics 1/sqrt lookup table
 881 #ifndef GMX_BLUEGENE
 882             rinv43           = frsqrtes(rsq43)
 883 #else
 884             rinv43           = frsqrte(dble(rsq43))
 885 #endif
 886             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 887      &  *rinv43)))
 888 #ifdef GMX_DOUBLE
 889             rinv43           = (0.5*rinv43*(3.0-((rsq43*rinv43)
 890      &  *rinv43)))
 891 #endif
 892
 893 C          PowerPC intrinsics 1/sqrt lookup table
 894 #ifndef GMX_BLUEGENE
 895             rinv44           = frsqrtes(rsq44)
 896 #else
 897             rinv44           = frsqrte(dble(rsq44))
 898 #endif
 899             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 900      &  *rinv44)))
 901 #ifdef GMX_DOUBLE
 902             rinv44           = (0.5*rinv44*(3.0-((rsq44*rinv44)
 903      &  *rinv44)))
 904 #endif
 905
 906 C          Load parameters for j atom
 907             qq               = qqHH
 908
 909 C          Coulomb interaction
 910             vcoul            = qq*rinv22
 911             vctot            = vctot+vcoul
 912
 913 C          Load parameters for j atom
 914             qq               = qqHH
 915
 916 C          Coulomb interaction
 917             vcoul            = qq*rinv23
 918             vctot            = vctot+vcoul
 919
 920 C          Load parameters for j atom
 921             qq               = qqMH
 922
 923 C          Coulomb interaction
 924             vcoul            = qq*rinv24
 925             vctot            = vctot+vcoul
 926
 927 C          Load parameters for j atom
 928             qq               = qqHH
 929
 930 C          Coulomb interaction
 931             vcoul            = qq*rinv32
 932             vctot            = vctot+vcoul
 933
 934 C          Load parameters for j atom
 935             qq               = qqHH
 936
 937 C          Coulomb interaction
 938             vcoul            = qq*rinv33
 939             vctot            = vctot+vcoul
 940
 941 C          Load parameters for j atom
 942             qq               = qqMH
 943
 944 C          Coulomb interaction
 945             vcoul            = qq*rinv34
 946             vctot            = vctot+vcoul
 947
 948 C          Load parameters for j atom
 949             qq               = qqMH
 950
 951 C          Coulomb interaction
 952             vcoul            = qq*rinv42
 953             vctot            = vctot+vcoul
 954
 955 C          Load parameters for j atom
 956             qq               = qqMH
 957
 958 C          Coulomb interaction
 959             vcoul            = qq*rinv43
 960             vctot            = vctot+vcoul
 961
 962 C          Load parameters for j atom
 963             qq               = qqMM
 964
 965 C          Coulomb interaction
 966             vcoul            = qq*rinv44
 967             vctot            = vctot+vcoul
 968
 969 C          Inner loop uses 144 flops/iteration
 970           end do
 971
 972
 973 C        Add i forces to mem and shifted force list
 974
 975 C        Add potential energies to the group for this list
 976           ggid             = gid(n)+1
 977           Vc(ggid)         = Vc(ggid) + vctot
 978
 979 C        Increment number of inner iterations
 980           ninner           = ninner + nj1 - nj0
 981
 982 C        Outer loop uses 10 flops/iteration
 983         end do
 984
 985
 986 C      Increment number of outer iterations
 987         nouter           = nouter + nn1 - nn0
 988       if(nn1.lt.nri) goto 10
 989
 990 C    Write outer/inner iteration count to pointers
 991       outeriter        = nouter
 992       inneriter        = ninner
 993       return
 994       end
 995
 996
 997