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