← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/tools/ga-5-2/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-2/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