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