← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/tools/ga-5-1/armci/tcgmsg/examples/md.f

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/ga-5-1/armci/tcgmsg/examples/md.f
1
 
c     
2
 
c     md example program due to dieter heermann
3
 
c     restructured and pdf added by rjh december 1990
4
 
c
5
 
c     the message passing version distributes the computation
6
 
c     of the forces (order npart**2) over the processes, assuming
7
 
c     that each process has all of the co-ordinates.
8
 
c     A global add then gives each process all the information
9
 
c     needed to compute the next update.
10
 
c     None of the order npart work has been parallelized so that
11
 
c     will begin to dominate on many processors.
12
 
c     
13
 
      implicit double precision (a-h, o-z)
14
 
c      parameter (mm = 13, lenpdf = 256)
15
 
c      parameter (mm = 8, lenpdf = 256)
16
 
c      parameter (mm = 6, lenpdf = 256)
17
 
C$Id: md.f,v 1.2 1995-02-02 23:24:18 d3g681 Exp $
18
 
      parameter (mm = 4, lenpdf = 256)
19
 
c      parameter (mm = 3, lenpdf = 256)
20
 
      parameter (npart = 4*mm*mm*mm)
21
 
      parameter (maxint = npart*150)
22
 
      include 'msgtypesf.h'
23
 
c     
24
 
      dimension  x(1:npart,1:3), vh(1:npart,1:3),f(1:npart,1:3),
25
 
     &     pdf(1:lenpdf+1), times(8), inter(2,maxint)
26
 
      data times/8*0.0d0/
27
 
c
28
 
c     initalize message passing environment
29
 
c
30
 
      call pbeginf
31
 
      me = nodeid()
32
 
c     
33
 
c     parameter definition (density, volume, temperature ...)
34
 
c     
35
 
      den = 0.83134d0
36
 
      side = (dble(npart) / den)**0.3333333d0
37
 
      tref = 0.722d0
38
 
      rcoff = min(3.5d0, side/2.0d0)
39
 
c     islow = 1 to match the rather large original timestep
40
 
      islow = 4
41
 
      h = 0.064d0 / islow
42
 
      irep = 50
43
 
      istop = 400
44
 
      iprint = 10
45
 
      ineigh = 10 * (1 + islow/4)
46
 
      movemx = 800
47
 
      delpdf = 0.5d0*side/lenpdf
48
 
      rdelp = 1.0d0 / delpdf
49
 
c
50
 
      if (me.eq.0)write(6,1) npart,side,rcoff,tref,h,delpdf,irep,istop,
51
 
     &     iprint,ineigh,movemx
52
 
 1    format(' molecular dynamics simulation example program'/
53
 
     &     ' ---------------------------------------------'//
54
 
     &     ' number of particles is ............ ',i6/
55
 
     &     ' side length of the box is ......... ',f13.6/
56
 
     &     ' cut off is ........................ ',f13.6/
57
 
     &     ' reduced temperature is ............ ',f13.6/
58
 
     &     ' basic time step is ................ ',f13.6/
59
 
     &     ' pdf sampling interval ............. ',f13.6/
60
 
     &     ' temperature scale interval ........ ',i6/
61
 
     &     ' stop scaling at move .............. ',i6/
62
 
     &     ' print interval .................... ',i6/
63
 
     &     ' update neighbor list every ........ ',i6, ' steps'/
64
 
     &     ' total no. of steps ................ ',i6)
65
 
c     call flush(6)
66
 
c     
67
 
      a = side / dble(mm)
68
 
      sideh = side * 0.5d0
69
 
      hsq = h * h
70
 
      hsq2 = hsq * 0.5d0
71
 
      npartm = npart - 1
72
 
      rcoffs = rcoff * rcoff
73
 
      tscale = 16.0d0 / (1.0d0 * npart - 1.0d0)
74
 
      vaver = 1.13d0 * sqrt(tref / 24.0d0)
75
 
      ekinavg = 0.0d0
76
 
c     
77
 
c     generate fcc lattice for atoms inside box
78
 
c     
79
 
      rjunk = timer()
80
 
      call fcc(x, npart, mm, a)
81
 
      times(1) = times(1) + timer()
82
 
c     
83
 
c     initialise velocites and forces (which are zero in fcc positions)
84
 
c     
85
 
      call mxwell(vh,3*npart,h,tref)
86
 
      call dfill(3*npart, 0.0d0, f, 1)
87
 
      times(2) = times(2) + timer()
88
 
c     
89
 
c     start of md. 
90
 
c     
91
 
      if (me.eq.0) write(6,3)
92
 
 3    format(//1x,'   i  ','     ke    ','      pe     ','      e     ',
93
 
     &     '    temp   ','   pres   ','   vel    ','  rp'/
94
 
     &     1x,' -----','  ----------','  ----------','  ----------',
95
 
     &        '  --------','  --------','  --------','  ----')
96
 
c     call flush(6)
97
 
c     
98
 
      do 200 move = 1,movemx
99
 
         if (move.eq.1 .or. move.eq.(istop+1))
100
 
     $        call dfill(lenpdf, 0.0d0, pdf, 1)
101
 
c     
102
 
c     move the particles and partially update velocities
103
 
c     
104
 
         call domove(3*npart,x,vh,f,side)
105
 
         times(3) = times(3) + timer()
106
 
c     
107
 
c     compute forces in the new positions and accumulate the pdf
108
 
c     virial and potential energy. Have to get the full forces
109
 
c     on each node, hence the dgop (global add).
110
 
c     
111
 
         if (mod(move-1,ineigh).eq.0) then
112
 
            call neigh(npart, x, side, rcoff, ninter, inter, maxint,
113
 
     $           pdf, rdelp, lenpdf)
114
 
            times(8) = times(8) + timer()
115
 
         endif
116
 
         call forces(npart, x, f, vir, epot, side, rcoff, ninter, inter)
117
 
         call dgop(1+MSGDBL, f, 3*npart, '+')
118
 
         times(4) = times(4) + timer()
119
 
c
120
 
c     scale forces, complete update of velocites and compute k.e.
121
 
c
122
 
         call mkekin(npart,f,vh,hsq2,hsq,ekin)
123
 
         ekinavg = ekinavg + ekin
124
 
         times(5) = times(5) + timer()
125
 
c
126
 
c     average the velocity and temperature scale if desired
127
 
c
128
 
         if ((move.le.istop) .and. (mod(move, irep).eq.0)) then
129
 
            call velavg(npart, vh, vaver, count, vel, h)
130
 
            sc = sqrt( tref / (tscale * ekinavg / irep) )
131
 
            call dscal(3*npart, sc, vh, 1)
132
 
            ekinavg = 0.0d0
133
 
         endif
134
 
         times(6) = times(6) + timer()
135
 
c
136
 
c     printout information if desired ... have to do global
137
 
c     sum to get full potential energy and virial
138
 
c
139
 
         if (mod(move, iprint) .eq. 0) then
140
 
            call velavg(npart, vh, vaver, count, vel, h)
141
 
            call dgop(2+MSGDBL, epot, 1, '+')
142
 
            call dgop(2+MSGDBL, vir, 1, '+')
143
 
            if (me.eq.0)
144
 
     $           call prnout(move, ekin, epot, tscale, vir, vel, count,
145
 
     $           npart, den)
146
 
            times(7) = times(7) + timer()
147
 
         endif
148
 
c     
149
 
 200  continue
150
 
c     
151
 
c     print out the pdf at the end of the calculation
152
 
c     have first to get contribution from everyone with global add
153
 
c     
154
 
      call dgop(2+MSGDBL, pdf, lenpdf, '+')
155
 
      if (me.eq.0)
156
 
     $     call prnpdf(lenpdf, pdf, side, delpdf, npart, movemx-istop,
157
 
     $     ineigh)
158
 
      times(7) = times(7) + timer()
159
 
      if (me.eq.0) write(6,431) nnodes(),(times(i),i=1,8)
160
 
 431  format('  nproc    geom  mxwell  domove  forces    ekin  velscl',
161
 
     &     '   print   neigh'
162
 
     &     /1x,i6,8f8.2)
163
 
c
164
 
      if (me.eq.0) call stats
165
 
      call pend
166
 
      call fexit
167
 
      end
168
 
      subroutine mxwell(vh,n3,h,tref)
169
 
      implicit double precision (a-h, o-z)
170
 
      dimension vh(1:n3)
171
 
c
172
 
c sample maxwell distribution at temperature tref
173
 
c
174
 
c alliant 3
175
 
      iseed = 4711
176
 
      ujunk = drand48(iseed)
177
 
      iseed = 0
178
 
      npart = n3/3
179
 
      iof1 = npart
180
 
      iof2 = npart*2
181
 
      tscale = 16.0d0 / (1.0d0 * npart - 1.0d0)
182
 
      do 10 i = 1,n3,2
183
 
c
184
 
c cray 2
185
 
c1         u1 = ranf()
186
 
c          u2 = ranf()
187
 
c alliant 2
188
 
1         u1 = drand48(iseed)
189
 
          u2 = drand48(iseed)
190
 
          v1 = 2.0d0 * u1 - 1.0d0
191
 
          v2 = 2.0d0 * u2 - 1.0d0
192
 
          s = v1*v1 + v2*v2
193
 
          if (s.ge.1.0) goto 1
194
 
c
195
 
          r = sqrt(-2.0d0*dlog(s)/s)
196
 
          vh(i) = v1 * r
197
 
          vh(i+1) = v2 * r
198
 
10    continue
199
 
c
200
 
      ekin = 0.0d0
201
 
      sp = 0.0d0
202
 
      do 20 i = 1,npart
203
 
          sp = sp + vh(i)
204
 
20    continue
205
 
      sp = sp / npart
206
 
      do 21 i = 1,npart
207
 
          vh(i) = vh(i) - sp
208
 
          ekin = ekin + vh(i)*vh(i)
209
 
21    continue
210
 
      sp = 0.0d0
211
 
      do 22 i = iof1 + 1,iof2
212
 
          sp = sp + vh(i)
213
 
22    continue
214
 
      sp = sp / npart
215
 
      do 23 i = iof1 + 1,iof2
216
 
          vh(i) = vh(i) - sp
217
 
          ekin = ekin + vh(i)*vh(i)
218
 
23    continue
219
 
      sp = 0.0d0
220
 
      do 24 i = iof2 + 1,n3
221
 
          sp = sp + vh(i)
222
 
24    continue
223
 
      sp = sp / npart
224
 
      do 25 i = iof2 + 1,n3
225
 
          vh(i) = vh(i) - sp
226
 
          ekin = ekin + vh(i)*vh(i)
227
 
25    continue
228
 
      ts = tscale * ekin
229
 
      sc = h * sqrt(tref/ts)
230
 
      do 30 i = 1,n3
231
 
          vh(i) = vh(i) * sc
232
 
30    continue
233
 
c
234
 
      end
235
 
      subroutine domove(n3,x,vh,f,side)
236
 
      implicit double precision (a-h, o-z)
237
 
      dimension x(n3),vh(n3),f(n3)
238
 
c
239
 
c     move particles
240
 
c
241
 
      do 10 i = 1,n3
242
 
         x(i) = x(i) + vh(i) + f(i)
243
 
c     periodic boundary conditions
244
 
         if (x(i).lt.0.0d0) x(i) = x(i) + side
245
 
         if (x(i).gt.side) x(i) = x(i) - side
246
 
c     partial velocity updates
247
 
         vh(i) = vh(i) + f(i)
248
 
c     initialise forces for next iteration
249
 
         f(i) = 0.0d0
250
 
 10   continue
251
 
c
252
 
      end    
253
 
      subroutine mkekin(npart,f,vh,hsq2,hsq,ekin)
254
 
      implicit double precision (a-h, o-z)
255
 
      dimension f(1:npart,3),vh(1:npart,3)
256
 
c
257
 
c     scale forces, update velocites and compute k.e.
258
 
c
259
 
      sum = 0.0d0
260
 
      do 10 ix = 1,3
261
 
         do 20 i = 1,npart
262
 
            f(i,ix) = f(i,ix) * hsq2
263
 
            vold = vh(i,ix)
264
 
            vh(i,ix) = vh(i,ix) + f(i,ix)
265
 
            sum = sum + vh(i,ix) * vh(i,ix)
266
 
 20      continue
267
 
 10   continue
268
 
      ekin = sum / hsq
269
 
c
270
 
      end
271
 
      subroutine fcc(x, npart, mm, a)
272
 
      implicit double precision (a-h, o-z)
273
 
      dimension x(1:npart, 3)
274
 
c     
275
 
c     generate fcc lattice for atoms inside box
276
 
c     
277
 
      ijk = 0
278
 
      do 10 lg = 0,1
279
 
         do 11 i = 0,mm-1
280
 
            do 12 j = 0,mm-1
281
 
               do 13 k = 0,mm-1
282
 
                  ijk = ijk + 1
283
 
                  x(ijk,1) = i * a + lg * a * 0.5d0
284
 
                  x(ijk,2) = j * a + lg * a * 0.5d0
285
 
                  x(ijk,3) = k * a
286
 
 13            continue
287
 
 12         continue
288
 
 11      continue
289
 
 10   continue
290
 
      do 20 lg = 1,2
291
 
         do 21 i = 0,mm-1
292
 
            do 22 j = 0,mm-1
293
 
               do 23 k = 0,mm-1
294
 
                  ijk = ijk + 1
295
 
                  x(ijk,1) = i * a + (2-lg) * a * 0.5d0
296
 
                  x(ijk,2) = j * a + (lg-1) * a * 0.5d0
297
 
                  x(ijk,3) = k * a + a * 0.5d0
298
 
 23            continue
299
 
 22         continue
300
 
 21      continue
301
 
20    continue
302
 
c
303
 
      end
304
 
      subroutine dfill(n,val,a,ia)
305
 
      implicit double precision (a-h, o-z)
306
 
      dimension a(*)
307
 
c
308
 
c     initialise double precision array to scalar value
309
 
c
310
 
      do 10 i = 1,(n-1)*ia+1,ia
311
 
         a(i) = val
312
 
 10   continue
313
 
c
314
 
      end
315
 
      subroutine prnout(move, ekin, epot, tscale, vir, vel, count,
316
 
     $     npart, den)
317
 
      implicit double precision (a-h, o-z)
318
 
c
319
 
c     printout interesting (?) information at current timestep
320
 
c
321
 
      ek = 24.0d0 * ekin
322
 
      epot = 4.0d0 * epot
323
 
      etot = ek + epot
324
 
      temp = tscale * ekin
325
 
      pres = den * 16.0d0 * (ekin - vir) / npart
326
 
      vel = vel / npart
327
 
      rp = (count / dble(npart)) * 100.0d0
328
 
      write(6,2) move,ek,epot,etot,temp,pres,vel,rp
329
 
 2    format(1x,i6,3f12.4,f10.4,f10.4,f10.4,f6.1)
330
 
c     call flush(6)
331
 
c
332
 
      end
333
 
      subroutine velavg(npart, vh, vaver, count, vel, h)
334
 
      implicit double precision (a-h, o-z)
335
 
      dimension vh(npart, 3)
336
 
c     
337
 
c     compute average velocity
338
 
c     
339
 
      vaverh = vaver*h
340
 
      vel = 0.0d0
341
 
      count = 0.0d0
342
 
      do 10 i = 1,npart
343
 
         sq = sqrt(vh(i,1)**2 + vh(i,2)**2 + vh(i,3)**2)
344
 
         if (sq.gt.vaverh) count = count + 1.0d0
345
 
         vel = vel + sq
346
 
 10   continue
347
 
      vel = vel / h
348
 
c     
349
 
      end
350
 
      subroutine prnpdf(lenpdf, pdf, side, delpdf, npart, nmove, ineigh)
351
 
      implicit double precision (a-h, o-z)
352
 
      dimension pdf(lenpdf)
353
 
c     
354
 
c     final scaling and printout of the pdf
355
 
c     
356
 
      write(6,1)
357
 
 1    format(/' pair distribution function written to file pdf.dat'/)
358
 
      open(1, file='pdf.dat', form='formatted', status='unknown',
359
 
     $     err=999)
360
 
c     
361
 
      coord = 0.0d0
362
 
      volfac = side*side*side / (4.0d0*delpdf*delpdf*delpdf*3.141593d0)
363
 
      facnn = 2.0d0 / dble(npart*nmove/ineigh)
364
 
      facn = 1.0d0 / dble(npart)
365
 
      do 10 i = 1,lenpdf
366
 
         ri = dble(i)
367
 
         grfac = volfac / (ri*ri)
368
 
         func = pdf(i) * facnn
369
 
         coord = coord + func
370
 
         pdf(i) = func * grfac * facn
371
 
         write(1,2) dble(i)*delpdf,pdf(i),coord
372
 
 10   continue
373
 
 2    format(1x,f7.3,f13.6,4x,f9.2)
374
 
      close(1)
375
 
      return
376
 
c
377
 
 999  write(6,*) ' error opening pdf.dat'
378
 
      call parerr(999)
379
 
c     
380
 
      end
381
 
      subroutine neigh(npart, x, side, rcoff, ninter, inter, maxint,
382
 
     $     pdf, rdelp, lenpdf)
383
 
      implicit double precision (a-h, o-z)
384
 
      dimension x(npart, 3), inter(2,maxint), pdf(lenpdf)
385
 
c     
386
 
c     Form my part of the neighbour list and also compute the pair
387
 
c     distribution function.
388
 
c     
389
 
c     npart = no. of particles
390
 
c     x(,)  = coords
391
 
c     side  = side of box
392
 
c     rcoff = cutoff for force
393
 
c     ninter= returns no. of interactions
394
 
c     inter(,) = returns interactions
395
 
c     maxint   = size of inter
396
 
c     
397
 
      me = nodeid()
398
 
      nproc = nnodes()
399
 
c     
400
 
      sideh = 0.5d0*side
401
 
      rcoffs = (rcoff*1.2d0)**2
402
 
      ninter = 0
403
 
c
404
 
c     Get better work distribution by having the same
405
 
c     processor handle particles (i) and (npart-i)
406
 
c     Note that assume that npart is even.
407
 
c
408
 
      do 270 ii = me+1, npart/2, nproc
409
 
         do 275 icase = 1, 2
410
 
            if (icase .eq. 1) then
411
 
               i = ii
412
 
            else
413
 
               i = npart - ii
414
 
            endif
415
 
            xi = x(i,1)
416
 
            yi = x(i,2)
417
 
            zi = x(i,3)
418
 
            do 280 j = i+1,npart
419
 
               ij = ij + 1
420
 
               xx = xi - x(j,1)
421
 
               yy = yi - x(j,2)
422
 
               zz = zi - x(j,3)
423
 
               if (xx .lt. -sideh) xx = xx + side
424
 
               if (xx .gt.  sideh) xx = xx - side
425
 
               if (yy .lt. -sideh) yy = yy + side
426
 
               if (yy .gt.  sideh) yy = yy - side
427
 
               if (zz .lt. -sideh) zz = zz + side
428
 
               if (zz .gt.  sideh) zz = zz - side
429
 
               rd = xx*xx + yy*yy + zz*zz
430
 
               ipdf = min(sqrt(rd),sideh) * rdelp + 1
431
 
               pdf(ipdf) = pdf(ipdf) + 1.0d0
432
 
               if (rd .le. rcoffs) then
433
 
                  ninter = ninter + 1
434
 
                  if (ninter.gt.maxint) then
435
 
                     write(6,*) ' too many interactions ', ninter
436
 
                     call parerr(1)
437
 
                  endif
438
 
                  inter(1,ninter) = i
439
 
                  inter(2,ninter) = j
440
 
               endif
441
 
 280        continue
442
 
 275     continue
443
 
 270  continue
444
 
c     
445
 
c$$$      ij = ninter
446
 
c$$$      call igop(99, ij, 1, '+')
447
 
c$$$      if (me .eq. 0) then
448
 
c$$$      write(6,*) ' No. of interactions per particle = ',
449
 
c$$$     $        ij/npart
450
 
c$$$      endif
451
 
c     
452
 
      end
453
 
      subroutine forces(npart, x, f, vir, epot, side, rcoff, 
454
 
     $     ninter, inter)
455
 
      implicit double precision (a-h, o-z)
456
 
      logical oshift
457
 
      parameter (oshift = .true.)
458
 
      dimension x(npart, 3), f(npart, 3), inter(2, ninter)
459
 
c     
460
 
c     compute forces driven by the neighbour list
461
 
c     
462
 
      vir = 0.0d0
463
 
      epot = 0.0d0
464
 
      sideh = 0.5d0*side
465
 
      rcoffs = rcoff*rcoff
466
 
c
467
 
c     for shifted potential ... set oshift true to enable
468
 
c
469
 
      if (oshift) then
470
 
         rc6  = 1.0d0 / rcoff**6
471
 
         rc12 = 1.0d0 / rcoff**12
472
 
         ecut = rc12 - rc6
473
 
         fcut = (rc12 - 0.5d0*rc6)/rcoff
474
 
         efcut = 12.0d0 * fcut
475
 
      endif
476
 
c
477
 
      do 10 ij = 1, ninter
478
 
         i = inter(1,ij)
479
 
         j = inter(2,ij)
480
 
         xx = x(i,1) - x(j,1)
481
 
         yy = x(i,2) - x(j,2)
482
 
         zz = x(i,3) - x(j,3)
483
 
         if (xx .lt. -sideh) xx = xx + side
484
 
         if (xx .gt.  sideh) xx = xx - side
485
 
         if (yy .lt. -sideh) yy = yy + side
486
 
         if (yy .gt.  sideh) yy = yy - side
487
 
         if (zz .lt. -sideh) zz = zz + side
488
 
         if (zz .gt.  sideh) zz = zz - side
489
 
         rd = xx*xx + yy*yy + zz*zz
490
 
         if (rd .le. rcoffs) then
491
 
            rrd = 1.0d0/rd
492
 
            rrd2 = rrd*rrd
493
 
            rrd3 = rrd2*rrd
494
 
            rrd4 = rrd2*rrd2
495
 
            rrd6 = rrd2*rrd4
496
 
            rrd7 = rrd6*rrd
497
 
            epot = epot + (rrd6 - rrd3)
498
 
            r148 = rrd7 - 0.5d0*rrd4
499
 
            if (oshift) then
500
 
               r = sqrt(rd)
501
 
               epot = epot - ecut + efcut*(r-rcoff)
502
 
               r148 = r148 - fcut / r
503
 
            endif
504
 
            vir = vir - rd*r148
505
 
            forcex = xx * r148
506
 
            f(i,1) = f(i,1) + forcex
507
 
            f(j,1) = f(j,1) - forcex
508
 
            forcey = yy * r148
509
 
            f(i,2) = f(i,2) + forcey
510
 
            f(j,2) = f(j,2) - forcey
511
 
            forcez = zz * r148
512
 
            f(i,3) = f(i,3) + forcez
513
 
            f(j,3) = f(j,3) - forcez
514
 
         endif
515
 
 10   continue
516
 
c     
517
 
      end