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