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