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