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