2 C This source code is part of
6 C Copyright (c) 1991-2000, University of Groningen, The Netherlands.
7 C Copyright (c) 2001-2009, The GROMACS Development Team
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
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
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.
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.
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!
36 # define gmxreal real*8
38 # define gmxreal real*4
44 C Gromacs nonbonded kernel pwr6kernel333
45 C Coulomb interaction: Tabulated
46 C VdW interaction: Tabulated
47 C water optimization: TIP4P - other atoms
48 C Calculate forces: yes
50 subroutine pwr6kernel333(
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
92 integer*4 n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
93 integer*4 nn0,nn1,nouter,ninner
95 gmxreal fscal,tx,ty,tz
97 gmxreal qq,vcoul,vctot
100 gmxreal Vvdw6,Vvdwtot
102 gmxreal r,rt,eps,eps2
104 gmxreal Y,F,Geps,Heps2,Fp,VV
108 gmxreal ix1,iy1,iz1,fix1,fiy1,fiz1
109 gmxreal ix2,iy2,iz2,fix2,fiy2,fiz2
110 gmxreal ix3,iy3,iz3,fix3,fiy3,fiz3
111 gmxreal ix4,iy4,iz4,fix4,fiy4,fiz4
112 gmxreal jx1,jy1,jz1,fjx1,fjy1,fjz1
113 gmxreal dx11,dy11,dz11,rsq11,rinv11
114 gmxreal dx21,dy21,dz21,rsq21,rinv21
115 gmxreal dx31,dy31,dz31,rsq31,rinv31
116 gmxreal dx41,dy41,dz41,rsq41,rinv41
121 C Initialize water data
123 qH = facel*charge(ii+1)
124 qM = facel*charge(ii+3)
125 nti = 2*ntype*type(ii)
128 C Reset outer and inner iteration counters
132 C Loop over thread workunits
133 10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
134 if(nn1.gt.nri) nn1=nri
136 C Start outer loop over neighborlists
140 C Load shift vector for this list
143 shY = shiftvec(is3+1)
144 shZ = shiftvec(is3+2)
146 C Load limits for loop over neighbors
150 C Get outer coordinate index
154 C Load i atom data, add shift vector
155 ix1 = shX + pos(ii3+0)
156 iy1 = shY + pos(ii3+1)
157 iz1 = shZ + pos(ii3+2)
158 ix2 = shX + pos(ii3+3)
159 iy2 = shY + pos(ii3+4)
160 iz2 = shZ + pos(ii3+5)
161 ix3 = shX + pos(ii3+6)
162 iy3 = shY + pos(ii3+7)
163 iz3 = shZ + pos(ii3+8)
164 ix4 = shX + pos(ii3+9)
165 iy4 = shY + pos(ii3+10)
166 iz4 = shZ + pos(ii3+11)
168 C Zero the potential energy for this list
172 C Clear i atom forces
188 C Get j neighbor index, and coordinate index
192 C load j atom coordinates
201 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11
205 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21
209 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31
213 rsq41 = dx41*dx41+dy41*dy41+dz41*dz41
215 C Calculate 1/r and 1/r2
217 C PowerPC intrinsics 1/sqrt lookup table
219 rinv11 = frsqrtes(rsq11)
221 rinv11 = frsqrte(dble(rsq11))
223 rinv11 = (0.5*rinv11*(3.0-((rsq11*rinv11)
226 rinv11 = (0.5*rinv11*(3.0-((rsq11*rinv11)
230 C PowerPC intrinsics 1/sqrt lookup table
232 rinv21 = frsqrtes(rsq21)
234 rinv21 = frsqrte(dble(rsq21))
236 rinv21 = (0.5*rinv21*(3.0-((rsq21*rinv21)
239 rinv21 = (0.5*rinv21*(3.0-((rsq21*rinv21)
243 C PowerPC intrinsics 1/sqrt lookup table
245 rinv31 = frsqrtes(rsq31)
247 rinv31 = frsqrte(dble(rsq31))
249 rinv31 = (0.5*rinv31*(3.0-((rsq31*rinv31)
252 rinv31 = (0.5*rinv31*(3.0-((rsq31*rinv31)
256 C PowerPC intrinsics 1/sqrt lookup table
258 rinv41 = frsqrtes(rsq41)
260 rinv41 = frsqrte(dble(rsq41))
262 rinv41 = (0.5*rinv41*(3.0-((rsq41*rinv41)
265 rinv41 = (0.5*rinv41*(3.0-((rsq41*rinv41)
269 C Load parameters for j atom
270 tj = nti+2*type(jnr)+1
274 C Calculate table index
277 C Calculate table index
284 C Tabulated VdW interaction - dispersion
288 Geps = eps*VFtab(nnn+2)
289 Heps2 = eps2*VFtab(nnn+3)
292 FF = Fp+Geps+2.0*Heps2
296 C Tabulated VdW interaction - repulsion
300 Geps = eps*VFtab(nnn+2)
301 Heps2 = eps2*VFtab(nnn+3)
304 FF = Fp+Geps+2.0*Heps2
307 Vvdwtot = Vvdwtot+ Vvdw6 + Vvdw12
308 fscal = -((fijD+fijR)*tabscale)*rinv11
310 C Calculate temporary vectorial force
315 C Increment i atom force
320 C Decrement j atom force
321 fjx1 = faction(j3+0) - tx
322 fjy1 = faction(j3+1) - ty
323 fjz1 = faction(j3+2) - tz
325 C Load parameters for j atom
329 C Calculate table index
332 C Calculate table index
339 C Tabulated coulomb interaction
342 Geps = eps*VFtab(nnn+2)
343 Heps2 = eps2*VFtab(nnn+3)
346 FF = Fp+Geps+2.0*Heps2
349 vctot = vctot + vcoul
350 fscal = -((fijC)*tabscale)*rinv21
352 C Calculate temporary vectorial force
357 C Increment i atom force
362 C Decrement j atom force
367 C Load parameters for j atom
369 C Calculate table index
372 C Calculate table index
379 C Tabulated coulomb interaction
382 Geps = eps*VFtab(nnn+2)
383 Heps2 = eps2*VFtab(nnn+3)
386 FF = Fp+Geps+2.0*Heps2
389 vctot = vctot + vcoul
390 fscal = -((fijC)*tabscale)*rinv31
392 C Calculate temporary vectorial force
397 C Increment i atom force
402 C Decrement j atom force
407 C Load parameters for j atom
410 C Calculate table index
413 C Calculate table index
420 C Tabulated coulomb interaction
423 Geps = eps*VFtab(nnn+2)
424 Heps2 = eps2*VFtab(nnn+3)
427 FF = Fp+Geps+2.0*Heps2
430 vctot = vctot + vcoul
431 fscal = -((fijC)*tabscale)*rinv41
433 C Calculate temporary vectorial force
438 C Increment i atom force
443 C Decrement j atom force
444 faction(j3+0) = fjx1 - tx
445 faction(j3+1) = fjy1 - ty
446 faction(j3+2) = fjz1 - tz
448 C Inner loop uses 183 flops/iteration
452 C Add i forces to mem and shifted force list
453 faction(ii3+0) = faction(ii3+0) + fix1
454 faction(ii3+1) = faction(ii3+1) + fiy1
455 faction(ii3+2) = faction(ii3+2) + fiz1
456 faction(ii3+3) = faction(ii3+3) + fix2
457 faction(ii3+4) = faction(ii3+4) + fiy2
458 faction(ii3+5) = faction(ii3+5) + fiz2
459 faction(ii3+6) = faction(ii3+6) + fix3
460 faction(ii3+7) = faction(ii3+7) + fiy3
461 faction(ii3+8) = faction(ii3+8) + fiz3
462 faction(ii3+9) = faction(ii3+9) + fix4
463 faction(ii3+10) = faction(ii3+10) + fiy4
464 faction(ii3+11) = faction(ii3+11) + fiz4
465 fshift(is3) = fshift(is3)+fix1+fix2+fix3+fix4
466 fshift(is3+1) = fshift(is3+1)+fiy1+fiy2+fiy3+fiy4
467 fshift(is3+2) = fshift(is3+2)+fiz1+fiz2+fiz3+fiz4
469 C Add potential energies to the group for this list
471 Vc(ggid) = Vc(ggid) + vctot
472 Vvdw(ggid) = Vvdw(ggid) + Vvdwtot
474 C Increment number of inner iterations
475 ninner = ninner + nj1 - nj0
477 C Outer loop uses 38 flops/iteration
481 C Increment number of outer iterations
482 nouter = nouter + nn1 - nn0
483 if(nn1.lt.nri) goto 10
485 C Write outer/inner iteration count to pointers
497 C Gromacs nonbonded kernel pwr6kernel333nf
498 C Coulomb interaction: Tabulated
499 C VdW interaction: Tabulated
500 C water optimization: TIP4P - other atoms
501 C Calculate forces: no
503 subroutine pwr6kernel333nf(
536 integer*4 nri,iinr(*),jindex(*),jjnr(*),shift(*)
537 gmxreal shiftvec(*),fshift(*),pos(*),faction(*)
538 integer*4 gid(*),type(*),ntype
539 gmxreal charge(*),facel,krf,crf,Vc(*),vdwparam(*)
540 gmxreal Vvdw(*),tabscale,VFtab(*)
541 gmxreal invsqrta(*),dvda(*),gbtabscale,GBtab(*)
542 integer*4 nthreads,count,mtx,outeriter,inneriter
545 integer*4 n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid
546 integer*4 nn0,nn1,nouter,ninner
549 gmxreal qq,vcoul,vctot
552 gmxreal Vvdw6,Vvdwtot
554 gmxreal r,rt,eps,eps2
556 gmxreal Y,F,Geps,Heps2,Fp,VV
562 gmxreal dx11,dy11,dz11,rsq11,rinv11
563 gmxreal dx21,dy21,dz21,rsq21,rinv21
564 gmxreal dx31,dy31,dz31,rsq31,rinv31
565 gmxreal dx41,dy41,dz41,rsq41,rinv41
570 C Initialize water data
572 qH = facel*charge(ii+1)
573 qM = facel*charge(ii+3)
574 nti = 2*ntype*type(ii)
577 C Reset outer and inner iteration counters
581 C Loop over thread workunits
582 10 call pwr6kernelsync(mtx,count,nri,nthreads,nn0,nn1)
583 if(nn1.gt.nri) nn1=nri
585 C Start outer loop over neighborlists
589 C Load shift vector for this list
592 shY = shiftvec(is3+1)
593 shZ = shiftvec(is3+2)
595 C Load limits for loop over neighbors
599 C Get outer coordinate index
603 C Load i atom data, add shift vector
604 ix1 = shX + pos(ii3+0)
605 iy1 = shY + pos(ii3+1)
606 iz1 = shZ + pos(ii3+2)
607 ix2 = shX + pos(ii3+3)
608 iy2 = shY + pos(ii3+4)
609 iz2 = shZ + pos(ii3+5)
610 ix3 = shX + pos(ii3+6)
611 iy3 = shY + pos(ii3+7)
612 iz3 = shZ + pos(ii3+8)
613 ix4 = shX + pos(ii3+9)
614 iy4 = shY + pos(ii3+10)
615 iz4 = shZ + pos(ii3+11)
617 C Zero the potential energy for this list
621 C Clear i atom forces
625 C Get j neighbor index, and coordinate index
629 C load j atom coordinates
638 rsq11 = dx11*dx11+dy11*dy11+dz11*dz11
642 rsq21 = dx21*dx21+dy21*dy21+dz21*dz21
646 rsq31 = dx31*dx31+dy31*dy31+dz31*dz31
650 rsq41 = dx41*dx41+dy41*dy41+dz41*dz41
652 C Calculate 1/r and 1/r2
654 C PowerPC intrinsics 1/sqrt lookup table
656 rinv11 = frsqrtes(rsq11)
658 rinv11 = frsqrte(dble(rsq11))
660 rinv11 = (0.5*rinv11*(3.0-((rsq11*rinv11)
663 rinv11 = (0.5*rinv11*(3.0-((rsq11*rinv11)
667 C PowerPC intrinsics 1/sqrt lookup table
669 rinv21 = frsqrtes(rsq21)
671 rinv21 = frsqrte(dble(rsq21))
673 rinv21 = (0.5*rinv21*(3.0-((rsq21*rinv21)
676 rinv21 = (0.5*rinv21*(3.0-((rsq21*rinv21)
680 C PowerPC intrinsics 1/sqrt lookup table
682 rinv31 = frsqrtes(rsq31)
684 rinv31 = frsqrte(dble(rsq31))
686 rinv31 = (0.5*rinv31*(3.0-((rsq31*rinv31)
689 rinv31 = (0.5*rinv31*(3.0-((rsq31*rinv31)
693 C PowerPC intrinsics 1/sqrt lookup table
695 rinv41 = frsqrtes(rsq41)
697 rinv41 = frsqrte(dble(rsq41))
699 rinv41 = (0.5*rinv41*(3.0-((rsq41*rinv41)
702 rinv41 = (0.5*rinv41*(3.0-((rsq41*rinv41)
706 C Load parameters for j atom
707 tj = nti+2*type(jnr)+1
711 C Calculate table index
714 C Calculate table index
721 C Tabulated VdW interaction - dispersion
725 Geps = eps*VFtab(nnn+2)
726 Heps2 = eps2*VFtab(nnn+3)
731 C Tabulated VdW interaction - repulsion
735 Geps = eps*VFtab(nnn+2)
736 Heps2 = eps2*VFtab(nnn+3)
740 Vvdwtot = Vvdwtot+ Vvdw6 + Vvdw12
742 C Load parameters for j atom
746 C Calculate table index
749 C Calculate table index
756 C Tabulated coulomb interaction
759 Geps = eps*VFtab(nnn+2)
760 Heps2 = eps2*VFtab(nnn+3)
764 vctot = vctot + vcoul
766 C Load parameters for j atom
768 C Calculate table index
771 C Calculate table index
778 C Tabulated coulomb interaction
781 Geps = eps*VFtab(nnn+2)
782 Heps2 = eps2*VFtab(nnn+3)
786 vctot = vctot + vcoul
788 C Load parameters for j atom
791 C Calculate table index
794 C Calculate table index
801 C Tabulated coulomb interaction
804 Geps = eps*VFtab(nnn+2)
805 Heps2 = eps2*VFtab(nnn+3)
809 vctot = vctot + vcoul
811 C Inner loop uses 114 flops/iteration
815 C Add i forces to mem and shifted force list
817 C Add potential energies to the group for this list
819 Vc(ggid) = Vc(ggid) + vctot
820 Vvdw(ggid) = Vvdw(ggid) + Vvdwtot
822 C Increment number of inner iterations
823 ninner = ninner + nj1 - nj0
825 C Outer loop uses 14 flops/iteration
829 C Increment number of outer iterations
830 nouter = nouter + nn1 - nn0
831 if(nn1.lt.nri) goto 10
833 C Write outer/inner iteration count to pointers