1#if HAVE_CONFIG_H
2#   include "config.fh"
3#endif
4      subroutine cluster_com
5#include "common.fh"
6c
7      integer i
8      double precision masstot, dx, dy, dz
9      double precision r,rdot,com(10)
10c
11c   This subroutine calculates the center of mass (COM) of the cluster of
12c   particles of type 2, excluding the lone particle.
13c
14      if (nocluster) return
15      cl_cmx = 0.0d00
16      cl_cmy = 0.0d00
17      cl_cmz = 0.0d00
18      cl_acmx = 0.0d00
19      cl_acmy = 0.0d00
20      cl_acmz = 0.0d00
21      cl_vcmx = 0.0d00
22      cl_vcmy = 0.0d00
23      cl_vcmz = 0.0d00
24      masstot = 0.0d00
25c
26      do i = 1, antot
27c      if (at(i).eq.2.and.aidx(i).ne.cl_lone_particle) then
28        if (at(i).eq.2) then
29          cl_cmx = cl_cmx + mass(i)*ra(i,1,6)
30          cl_cmy = cl_cmy + mass(i)*ra(i,2,6)
31          cl_cmz = cl_cmz + mass(i)*ra(i,3,6)
32          cl_acmx = cl_acmx + ra(i,1,4)
33          cl_acmy = cl_acmy + ra(i,2,4)
34          cl_acmz = cl_acmz + ra(i,3,4)
35          cl_vcmx = cl_vcmx + mass(i)*ra(i,1,2)
36          cl_vcmy = cl_vcmy + mass(i)*ra(i,2,2)
37          cl_vcmz = cl_vcmz + mass(i)*ra(i,3,2)
38          masstot = masstot + mass(i)
39        endif
40      end do
41      com(1) = cl_cmx
42      com(2) = cl_cmy
43      com(3) = cl_cmz
44      com(4) = cl_vcmx
45      com(5) = cl_vcmy
46      com(6) = cl_vcmz
47      com(7) = cl_acmx
48      com(8) = cl_acmy
49      com(9) = cl_acmz
50      com(10) = masstot
51      call ga_dgop(3,com,10,'+')
52      cl_cmx  = com(1)
53      cl_cmy  = com(2)
54      cl_cmz  = com(3)
55      cl_acmx  = com(7)
56      cl_acmy  = com(8)
57      cl_acmz  = com(9)
58      cl_vcmx  = com(4)
59      cl_vcmy  = com(5)
60      cl_vcmz  = com(6)
61      masstot = com(10)
62      if (masstot.gt.0.0d00) then
63        cl_cmx = cl_cmx / masstot
64        cl_cmy = cl_cmy / masstot
65        cl_cmz = cl_cmz / masstot
66        cl_acmx = cl_acmx / masstot
67        cl_acmy = cl_acmy / masstot
68        cl_acmz = cl_acmz / masstot
69        cl_vcmx = cl_vcmx / masstot
70        cl_vcmy = cl_vcmy / masstot
71        cl_vcmz = cl_vcmz / masstot
72      endif
73      cl_mass = masstot
74c
75      return
76      end
77c
78      subroutine cluster_therm
79#include "common.fh"
80c
81      integer i
82      double precision dx, dy, dz
83      double precision r,rdot
84c
85c  remove component of velocity that lies along the vector from COM
86c  to the particle, if particle is outside cluster radius. This gradually
87c  forces particles to form a drop.
88c
89      if (nocluster) return
90      if (istep.ge.equil_1) return
91      call cluster_com
92      do i = 1, antot
93c        if (at(i).eq.2.and.aidx(i).ne.cl_lone_particle) then
94        if (at(i).eq.2) then
95          dx = ra(i,1,1) - cl_cmx
96          dy = ra(i,2,1) - cl_cmy
97          dz = ra(i,3,1) - cl_cmz
98          r = sqrt(dx**2 + dy**2 + dz**2)
99          dx = dx/r
100          dy = dy/r
101          dz = dz/r
102          if (r.gt.r_cluster) then
103            rdot = ra(i,1,2)*dx+ra(i,2,2)*dy+ra(i,3,2)*dz
104            if (rdot.gt.0.0d00) then
105              ra(i,1,2) = ra(i,1,2) - 2.0d00*rdot*dx
106              ra(i,2,2) = ra(i,2,2) - 2.0d00*rdot*dy
107              ra(i,3,2) = ra(i,3,2) - 2.0d00*rdot*dz
108            endif
109          endif
110        endif
111      end do
112c
113      return
114      end
115c
116      subroutine cluster_center
117#include "common.fh"
118c
119c  move center of mass to origin (only perform this if it is immediately
120c  followed by a call to update subroutine)
121c
122      integer i,me
123      logical debug
124      if (nocluster) return
125      if (istep.eq.0.or.(mod(istep,ilist).eq.0.and.
126     +    t_rmndr.eq.0.0d00)) then
127        me = ga_nodeid()
128        if (istep.ge.6932366) then
129          debug = .false.
130        else
131          debug = .false.
132        endif
133        if (debug) then
134          write(6,*) me,' mod(istep,ilist) = ',mod(istep,ilist)
135          write(6,*) me,' ilist = ',ilist
136        endif
137        call cluster_com
138        if (debug) write(6,*) me,'cl_cmx (a) = ',cl_cmx,istep
139        if (debug) write(6,*) me,'cl_cmy (a) = ',cl_cmy,istep
140        if (debug) write(6,*) me,'cl_cmz (a) = ',cl_cmz,istep
141        do i = 1, antot
142          ra(i,1,6) = ra(i,1,6) - cl_cmx
143          ra(i,2,6) = ra(i,2,6) - cl_cmy
144          ra(i,3,6) = ra(i,3,6) - cl_cmz
145        end do
146        call fixper
147        call cluster_com
148        if (debug) write(6,*) me,'cl_cmx (b) = ',cl_cmx,istep
149        if (debug) write(6,*) me,'cl_cmy (b) = ',cl_cmy,istep
150        if (debug) write(6,*) me,'cl_cmz (b) = ',cl_cmz,istep
151      endif
152      return
153      end
154c
155      subroutine cluster_old_at
156#include "common.fh"
157c
158c   store original coordinates of cluster atoms
159c
160      integer i, icnt, jcnt
161      if (nocluster) return
162      icnt = 0
163      jcnt = 0
164      do i = 1, antot
165        if (at(i).eq.2) then
166          icnt = icnt + 1
167          cl_at(icnt) = i
168          cl_old(icnt,1,1) = ra(i,1,1)
169          cl_old(icnt,2,1) = ra(i,2,1)
170          cl_old(icnt,3,1) = ra(i,3,1)
171          cl_old(icnt,1,2) = ra(i,1,6)
172          cl_old(icnt,2,2) = ra(i,2,6)
173          cl_old(icnt,3,2) = ra(i,3,6)
174          cl_old(icnt,1,3) = ra(i,1,2)
175          cl_old(icnt,2,3) = ra(i,2,2)
176          cl_old(icnt,3,3) = ra(i,3,2)
177        endif
178        if (at(i).eq.1) then
179          jcnt = jcnt + 1
180          sl_at(jcnt) = i
181          sl_old(jcnt,1,1) = ra(i,1,1)
182          sl_old(jcnt,2,1) = ra(i,2,1)
183          sl_old(jcnt,3,1) = ra(i,3,1)
184          sl_old(jcnt,1,2) = ra(i,1,6)
185          sl_old(jcnt,2,2) = ra(i,2,6)
186          sl_old(jcnt,3,2) = ra(i,3,6)
187          sl_old(jcnt,1,3) = ra(i,1,2)
188          sl_old(jcnt,2,3) = ra(i,2,2)
189          sl_old(jcnt,3,3) = ra(i,3,2)
190        endif
191      end do
192c
193      cl_cm_old(1) = cl_cmx
194      cl_cm_old(2) = cl_cmy
195      cl_cm_old(3) = cl_cmz
196      cl_vcm_old(1) = cl_vcmx
197      cl_vcm_old(2) = cl_vcmy
198      cl_vcm_old(3) = cl_vcmz
199c
200      do i = 1, 3
201        cl_alen1_old(i) = alen1(i)
202        cl_alen2_old(i) = alen2(i)
203      end do
204      cl_vol1_old = vol1
205      cl_vol2_old = vol2
206      cl_box_old(1) = xbox
207      cl_box_old(2) = ybox
208      cl_box_old(3) = zbox
209      cl_scal1_old = scal1
210      cl_scal2_old = scal2
211c
212      cl_tot = icnt
213      sl_tot = jcnt
214      return
215      end
216c
217      subroutine cluster_check_cllsn
218#include "common.fh"
219c
220c   This subroutine checks to see if any of the cluster particles have moved
221c   beyond the cutoff radius from the center of mass of the cluster. If they
222c   have, then the initial guess for the particle coordinates and velocities
223c   is modified to reflect the collision with the restraining sphere.
224c
225      integer icnt, iat, i, j, ii, jj, iloc, imin, scndat, reset
226      double precision r, r2, rx, ry, rz, dtmax, tsav
227      double precision rrx, rry, rrz, vvx, vvy, vvz
228      double precision vn, vx, vy, vz, fx, fy, fz, mu, v2, f2
229      double precision vrdot, frdot, vfdot, a, b, c
230      double precision t1, t2, tcut, comm(4), tmax, htausq
231      double precision rnx, rny, rnz, vpx, vpy, vpz, vllx, vlly, vllz
232      double precision t3, t4, ax, ay, az, tchk, tmin, t_est
233      double precision r1x, r1y, r1z, r2x, r2y, r2z, rmax, rskin, rn
234      double precision cluster_find_tau
235      integer failsafe, ibuf(MD_MAXPROC), me, nproc, twohit
236      logical debug
237      integer iter
238c
239      iter = 0
240      l_cllsn = .false.
241      if (nocluster) then
242        t_rmndr = 0.0d00
243        t_done = tau
244        return
245      endif
246      me = ga_nodeid()
247      nproc = ga_nnodes()
248      if (istep.ge.68365721.and.istep.le.68365723) then
249        debug = .false.
250      else
251        debug = .false.
252      endif
253      rskin = 0.02d00
254      tmax = t_rmndr
255      if (t_rmndr.eq.tau) cllsn_isav = 0
256  100 dtmax = 2.0d00*t_rmndr
257      twohit = 0
258      failsafe = 0
259      iat = 0
260      iloc = 0
261      tsav = 0.0d00
262      rmax = r_cluster
263      if (debug) write(6,*) me,' atot = ',atot,istep
264      if (debug) write(6,*) me,' antot = ',antot,istep
265      if (debug) write(6,*) me, ' r_cluster = ',r_cluster,istep
266      if (debug) write(6,*) me,' t_rmndr = ',tmax,istep
267      if (debug) write(6,*) me,' cl_cmx(1) = ',cl_cmx,istep
268      if (debug) write(6,*) me,' cl_cmy(1) = ',cl_cmy,istep
269      if (debug) write(6,*) me,' cl_cmz(1) = ',cl_cmz,istep
270      do i = 1, cl_tot
271        ii = cl_at(i)
272        rx = ra(ii,1,6) - cl_cmx
273        ry = ra(ii,2,6) - cl_cmy
274        rz = ra(ii,3,6) - cl_cmz
275        r2 = rx**2 + ry**2 + rz**2
276        r = sqrt(r2)
277c        if (debug) write(6,*) me, 'ra(ii,1,6) ',ra(ii,1,6),istep
278c        if (debug) write(6,*) me, 'ra(ii,2,6) ',ra(ii,2,6),istep
279c        if (debug) write(6,*) me, 'ra(ii,3,6) ',ra(ii,3,6),istep
280        if (r.gt.r_cluster) then
281          if (r.gt.rmax) rmax = r
282          if (debug) write(6,*) me,' tmax = ',tmax,istep
283          if (debug) write(6,*) me,' dtmax = ',dtmax,istep
284          if (debug) write(6,*) me,' ii = ',ii,istep
285          if (debug) write(6,*) me,' antot = ',antot,istep
286          if (debug) write(6,*) me,' r-r_cluster ',r-r_cluster,istep
287          if (debug) write(6,*) me,' r_cluster = ',r_cluster,istep
288          if (debug) write(6,*) me,' sqrt(r2) = ',sqrt(r2),istep
289          if (debug) write(6,*) me,' ra(1) = ',ra(ii,1,6),istep
290          if (debug) write(6,*) me,' ra(2) = ',ra(ii,2,6),istep
291          if (debug) write(6,*) me,' ra(3) = ',ra(ii,3,6),istep
292          if (debug) write(6,*) me,' cl_cmx = ',cl_cmx,istep
293          if (debug) write(6,*) me,' cl_cmy = ',cl_cmy,istep
294          if (debug) write(6,*) me,' cl_cmz = ',cl_cmz,istep
295          if (debug) write(6,*) me,' r_cluster_old = ',r_cluster_old,
296     +      istep
297          if (debug) write(6,*) me,' cl_old(1) = ',cl_old(i,1,2),istep
298          if (debug) write(6,*) me,' cl_old(2) = ',cl_old(i,2,2),istep
299          if (debug) write(6,*) me,' cl_old(3) = ',cl_old(i,3,2),istep
300          if (debug) write(6,*) me,' cl_cm_old(1) = ',cl_cm_old(1),istep
301          if (debug) write(6,*) me,' cl_cm_old(2) = ',cl_cm_old(2),istep
302          if (debug) write(6,*) me,' cl_cm_old(3) = ',cl_cm_old(3),istep
303          if (debug) write(6,*) me,' r_old = ', sqrt(
304     +      + (cl_old(i,1,2)-cl_cm_old(1))**2
305     +      + (cl_old(i,2,2)-cl_cm_old(2))**2
306     +      + (cl_old(i,3,2)-cl_cm_old(3))**2)
307          rn = sqrt(rx**2+ry**2+rz**2)
308          vn = ra(ii,1,2)*rx+ra(ii,2,2)*ry+ra(ii,3,2)*rz
309          vn = vn/rn
310          t_est = vn*(rn-r_cluster)
311          tcut = cluster_find_tau(ii,i,cllsn_isav,tmax,.false.)
312          if (debug) write(6,*) me,' tcut = ',tcut,istep
313          if (debug) write(6,*) me,' t_est = ',t_est,istep
314          if (debug) write(6,155) me,1,tcut,istep
315  155     format(i3,' tcut at ',i1,': ',f12.4,' step: ',i8)
316c
317c   Check to see if collision is earlier than previously found collisions
318c   and that it is also greater than zero.
319c
320          if (tcut.lt.dtmax.and.tcut.gt.0.0d00) then
321c
322c   Try to protect against numerical roundoff by checking that different
323c   particle collides with sphere or that collision time is significantly
324c   greater than zero if it is the same particle
325c
326            if (.not.(ii.eq.cllsn_isav.and.tcut.lt.1.0d-03)) then
327              dtmax = tcut
328              tsav = tcut
329              iat = ii
330              iloc = i
331            else
332              write(6,101) ii, tcut, istep
333  101         format('Rejected collision of atom ',i8,' at time ',
334     +               f16.8,' at step ',i8)
335              failsafe = 1
336            endif
337          else if (tcut.lt.0.0d00) then
338            write(6,*) ga_nodeid(),' vn ',vn,istep
339            write(6,*) ga_nodeid(),' t_est ',t_est,istep
340            write(6,*) ga_nodeid(),' r-r_cluster ',r-r_cluster,istep
341            write(6,*) ga_nodeid(),' Returned negative value at 1',istep
342            failsafe = 1
343          endif
344#if 0
345c
346c   Check for trajectories that may have gone out and then back into
347c   confining sphere
348c
349        else if (r_cluster-r.lt.rskin.and.r_cluster-r.gt.0.0d00) then
350          vn = rx*ra(ii,1,2)+ry*ra(ii,2,2)+rz*ra(ii,3,2)
351          if (vn.le.0.0d00) then
352            tcut = cluster_find_tau(ii,i,cllsn_isav,tmax,.true.)
353            if (tcut.lt.dtmax.and.tcut.gt.0.0d00) then
354              write(6,*) 'Found looping trajectory',istep
355              write(6,*) 'dtmax = ',dtmax,istep
356              write(6,*) 'tcut = ',tcut,istep
357              write(6,*) 'ii = ',ii,istep
358              write(6,*) 'cllsn_isav = ',cllsn_isav,istep
359              write(6,*) 'r = ',r,istep
360              write(6,*) 'r_cluster = ',r_cluster,istep
361              if (.not.(ii.eq.cllsn_isav.and.tcut.lt.1.0d-10)) then
362                write(6,*) 'Resetting trajectory',istep
363                dtmax = tcut
364                tsav = tcut
365                iat = ii
366                iloc = i
367              else
368                write(6,*)'Not resetting trajectory (loop is too short)'
369              endif
370            endif
371          endif
372#endif
373        endif
374      end do
375c
376c   Check to see if a collision occured on another processor
377c
378      ibuf(1) = iat
379      ibuf(2) = failsafe
380      if (debug) write(6,*) me,' iat = ',iat,istep
381      if (debug) write(6,*) me,' failsafe = ',failsafe,istep
382      if (debug) write(6,*) me,' ibuf(1) = ',ibuf(1),istep
383      if (debug) write(6,*) me,' ibuf(2) = ',ibuf(2),istep
384      call ga_igop(5,ibuf,2,'+')
385      if (debug) write(6,*) me,' ibuf(1) = ',ibuf(1),istep
386      if (debug) write(6,*) me,' ibuf(2) = ',ibuf(2),istep
387c
388c   No collisions detected. Step is complete so just return.
389c
390      if (ibuf(1).eq.0.and.ibuf(2).eq.0) then
391        t_rmndr = 0.0d00
392        t_done = tau
393        if (debug) write(6,*) me,' Returning with no hits ',istep
394        return
395      endif
396c
397c   If no solution found for atom outside confining sphere, then bail
398c   completely and increase diameter of confining sphere so that all
399c   atoms are enclosed.
400c
401      if (ibuf(2).gt.0) then
402        call ga_dgop(4,rmax,1,'max')
403        if (debug) write(6,*) me,' r_cluster (a) = ',r_cluster
404        r_cluster = r_cluster + (rmax - r_cluster)*1.1d00
405        if (debug) write(6,*) me,' r_cluster (b) = ',r_cluster
406        if (ga_nodeid().eq.0) then
407          write(6,*) ga_nodeid(),' Returned negative value at 2',istep
408        endif
409        failcount = failcount + 1
410        t_rmndr = 0.0d00
411        t_done = tau
412        return
413      endif
414c
415c   At this point, a collision has been detected on at least on
416c   processor. Check to find minimum time if collisions occur on more
417c   than one processor.
418c
419      if (iat.eq.0) then
420        tsav = tau
421      endif
422      a = tsav
423      call ga_dgop(6,a,1,'min')
424c
425c   Find mass of particle that has first collision with confining sphere
426c
427      if (a.eq.tsav) then
428        b = mass(iat)
429      else
430        b = 0.0d00
431      endif
432      call ga_dgop(7,b,1,'max')
433 1000 mmass = cl_mass - b
434      tcut = a
435      if (debug) write(6,155) me,2,tcut,istep
436      if (tsav.ne.a) then
437        iat = 0
438        iloc = 0
439        cllsn_isav = 0
440        if (debug) write(6,*) me,' cllsn_isav(1): ',cllsn_isav
441      else
442        cllsn_isav = iat
443        if (debug) write(6,*) me,' cllsn_isav(2): ',cllsn_isav
444      endif
445c
446c  Handle degenerate case of only two atoms in cluster, which will
447c  generally have two simultaneous collisions
448c
449      if (ctot.eq.2) then
450        do i = 1, nproc
451          if (i-1.eq.me) then
452            ibuf(i) = iat
453          else
454            ibuf(i) = 0
455          endif
456        end do
457        call ga_igop(3,ibuf,nproc,'+')
458        j = 0
459        iat = 0
460        cllsn_isav = 0
461        do i = 1, nproc
462          if (ibuf(i).gt.0.and.i-1.eq.me.and.j.eq.0) then
463            iat = ibuf(i)
464            cllsn_isav = ibuf(i)
465            j = 1
466          else if (ibuf(i).gt.0) then
467            j = 1
468          endif
469        end do
470c
471c  Make sure that correct mass is being used for collision
472c
473        if (iat.ne.0) then
474          b = mass(iat)
475        else
476          b = 0.0d00
477        endif
478        call ga_dgop(6,b,1,'max')
479        mmass = cl_mass - b
480      endif
481c
482c   Now recalculate coordinates of all particles in the cluster
483c
484      call cluster_reset(tcut)
485      call cluster_com !CHECK
486c
487c   Check to see if there were any trajectories that pass through confining
488c   sphere twice (i.e. goes out and then comes back in). Don't bother with
489c   check if there are only 2 particles in system.
490c
491      if (ctot.eq.2) then
492        go to 2000
493      endif
494      reset = 0
495      tchk = 0.0d00
496      scndat = 0
497      tmin = tcut
498      tsav = tcut
499      failsafe = 0
500      do i = 1, cl_tot
501        ii = cl_at(i)
502        rx = ra(ii,1,6) - cl_cmx
503        ry = ra(ii,2,6) - cl_cmy
504        rz = ra(ii,3,6) - cl_cmz
505        r = sqrt(rx**2+ry**2+rz**2)
506        if (r.ge.r_cluster.and.ii.ne.cllsn_isav) then
507          tchk = cluster_find_tau(ii,i,cllsn_isav,tcut,.false.)
508          if (debug) write(6,155) me,7,tchk,istep
509          if (debug) write(6,*) me, 'ii: ',ii
510          if (debug) write(6,*) me, 'cllsn_isav: ',cllsn_isav
511          reset = 1
512          if (tchk.lt.tmin.and.tchk.gt.0.0d00) then
513            if (debug) write(6,*) me,' Shorter collision found'
514            tmin = tchk
515            imin = ii
516            if (debug) write(6,155) me,4,tchk,istep
517            if (debug) write(6,155) me,5,tmin,istep
518          else if (tchk.lt.0.0d00) then
519            write(6,*) ga_nodeid(),' Returned negative value at 3',istep
520            failsafe = 1
521          else
522c
523c  This case should theoretically be impossible, but it may occur because of
524c  roundoff
525c
526            write(6,*) me,' tchk greater than or equal to tmin',istep
527            write(6,*) me,' tchk = ',tchk,istep
528            write(6,*) me,' tmin = ',tmin,istep
529            failsafe = 1
530          endif
531        endif
532      end do
533      ibuf(1) = reset
534      ibuf(2) = failsafe
535      call ga_igop(2,ibuf,2,'+')
536c
537c   No solution for tau found, so increase confining sphere and then bail
538c
539      if (ibuf(2).gt.0) then
540        if (debug) write(6,*) me,' rmax (c) = ',rmax
541        call ga_dgop(4,rmax,1,'max')
542        if (debug) write(6,*) me,' rmax (d) = ',rmax
543        if (debug) write(6,*) me,' r_cluster (c) = ',r_cluster
544        r_cluster = r_cluster + (rmax - r_cluster)*1.1d00
545        if (debug) write(6,*) me,' r_cluster (d) = ',r_cluster
546        if (ga_nodeid().eq.0) then
547          write(6,*) 'Bogus value of collision time at 2'
548        endif
549        call cluster_reset(t_rmndr)
550        t_rmndr = 0.0d00
551        t_done = tau
552        return
553      endif
554      reset = ibuf(1)
555c
556c   Shorter collision found so recalculate positions
557c
558      if (reset.gt.0) then
559        if (tmin.le.tsav) then
560          if (debug) write(6,*) 'tmin = ',tmin,istep
561          if (debug) write(6,*) 'tsav = ',tsav,istep
562          a = tmin
563        else
564          a = tau
565        endif
566        call ga_dgop(2,a,1,'min')
567c
568        if (a.le.tsav) then
569          if (a.ne.tmin) then
570            iat = 0
571            iloc = 0
572            cllsn_isav = 0
573            if (debug) write(6,*) me,'cllsn_isav(3): ',cllsn_isav
574          else
575            cllsn_isav = imin
576            if (debug) write(6,*) me,'cllsn_isav(4): ',cllsn_isav
577          endif
578          tcut = a
579          tsav = a   !BJP is this correct?
580          if (debug) write(6,155) me,3,tcut,istep
581          if (iat.gt.0) then
582            b = mass(iat)
583          else
584            b = 0.0d00
585          endif
586          call ga_dgop(7,b,1,'+')
587          call cluster_reset(0.0d00)
588          iter = iter + 1
589          if (iter.gt.100) debug = .true.
590          if (iter.gt.1000) call ga_error('Too many iterations',0)
591          go to 1000
592        endif
593      endif
594c
595 2000 l_cllsn = .true.
596      cllsn_idx = cllsn_isav
597c
598      t_rmndr = tmax - tcut
599      t_done = tau - t_rmndr
600      if (debug) write(6,*) me, 't_rmndr = ',t_rmndr
601      if (debug) write(6,*) me, 't_done = ',t_done
602c
603      return
604      end
605c
606      double precision function cluster_find_tau(iat,iloc,isav,tmax,
607     +                                           forward)
608#include "common.fh"
609c
610c   This function calculates the time at which a particle outside the confining
611c   sphere would have collided with the sphere.
612c
613      integer iat, iloc, isav
614      double precision r, r2, rx, ry, rz, dr
615      double precision rrx, rry, rrz, vvx, vvy, vvz
616      double precision vx, vy, vz, fx, fy, fz, v2, f2
617      double precision clcut, frdot, vrdot, vfdot, a, b, c
618      double precision t, t1, t2, tcut, tmax, thi, tlo, tmid
619      double precision ttmin,ttmax
620      integer i,j,me,iter
621      logical forward, debug
622c
623c  A particle lies outside the cutoff. Find particle that crosses
624c  cutoff first.
625c
626      if (istep.ge.68365721.and.istep.le.68365723) then
627        debug = .false.
628      else
629        debug = .false.
630      endif
631      me = ga_nodeid()
632      mmass = cl_mass - mass(iat)
633c
634c  Calculate center of mass and velocity of center of mass of remaining
635c  particles in cluster (these are the rr and vv vectors)
636c
637      rrx = (cl_mass*cl_cm_old(1) - mass(iat)*cl_old(iloc,1,2))/mmass
638      rry = (cl_mass*cl_cm_old(2) - mass(iat)*cl_old(iloc,2,2))/mmass
639      rrz = (cl_mass*cl_cm_old(3) - mass(iat)*cl_old(iloc,3,2))/mmass
640      vvx = (cl_mass*cl_vcm_old(1) - mass(iat)*cl_old(iloc,1,3))/mmass
641      vvy = (cl_mass*cl_vcm_old(2) - mass(iat)*cl_old(iloc,2,3))/mmass
642      vvz = (cl_mass*cl_vcm_old(3) - mass(iat)*cl_old(iloc,3,3))/mmass
643c
644c  The r and v vectors are the relative coordinates and velocities of
645c  particle iat and the center of mass vectors. The vector f is the
646c  force along this relative coordinate.
647c
648      rx = cl_old(iloc,1,2) - cl_cm_old(1)
649      ry = cl_old(iloc,2,2) - cl_cm_old(2)
650      rz = cl_old(iloc,3,2) - cl_cm_old(3)
651      vx = cl_old(iloc,1,3) - cl_vcm_old(1)
652      vy = cl_old(iloc,2,3) - cl_vcm_old(2)
653      vz = cl_old(iloc,3,3) - cl_vcm_old(3)
654      fx = ra(iat,1,3) - cl_acmx
655      fy = ra(iat,2,3) - cl_acmy
656      fz = ra(iat,3,3) - cl_acmz
657c
658      r = sqrt(rx**2 + ry**2 + rz**2)
659c
660c  Note that we are using acceleration instead of force,
661c  so factors of reduced mass disappear
662c
663      vrdot = vx*rx + vy*ry + vz*rz
664      frdot = fx*rx + fy*ry + fz*rz
665      vfdot = vx*fx + vy*fy + vz*fz
666      clcut = r_cluster
667      v2 = vx**2 + vy**2 + vz**2
668      a = v2+frdot
669      b = 2.0d00*vrdot
670c
671c      c = rx**2 + ry**2 + rz**2 - clcut**2
672c
673c  Calculate c to minimize roundoff error from subtracting large numbers
674c
675      dr = r - clcut
676      c = 2.0d00*r*dr+dr**2
677      f2 = fx**2 + fy**2 + fz**2
678c
679c    Scan forwards from 0 to find a final value of t2. If no final
680c    value found, then return with value of -1.
681c
682              if (debug) write(6,*) me, 'iat: ',iat
683              if (debug) write(6,*) me, 'a: ',a
684              if (debug) write(6,*) me, 'b: ',b
685              if (debug) write(6,*) me, 'c: ',c
686              if (debug) write(6,*) me, 'f2: ',f2
687              if (debug) write(6,*) me, 'vfdot: ',vfdot
688      if (forward) then
689        t = tmax
690        t1 = 0.0d00
691        ttmin = t1
692        ttmax = t
693c        tlo = c + t1*(b + t1*(a+t1*(vfdot+t1*0.25d00*f2)))
694        tlo = (-c-t1**2*(a+t1*(vfdot+t1*0.25d00*f2)))/b - t1
695        do i = 1, 1000
696          if (t.eq.tmax) then
697            t2= tmax*dble(i)/1000.0d00
698c            thi = c + t2*(b + t2*(a+t2*(vfdot+t2*0.25d00*f2)))
699            thi = (-c-t2**2*(a+t2*(vfdot+t2*0.25d00*f2)))/b - t2
700            if ((thi.ge.0.0d00.and.tlo.le.0.0d00).or.
701     +          (thi.le.0.0d00.and.tlo.ge.0.0d00)) then
702              t = t2
703              write(6,*) ga_nodeid(),' tmin = ',ttmin,istep
704              write(6,*) ga_nodeid(),' tmax = ',ttmax,istep
705              write(6,*) ga_nodeid(),' tlo = ',tlo,istep
706              write(6,*) ga_nodeid(),' thi = ',thi,istep
707              write(6,*) ga_nodeid(),' t = ',t,istep
708              go to 500
709            endif
710          endif
711        end do
712        cluster_find_tau = -1.0d00
713        return
714c
715c    scan backwards from tmax to find an initial value of t1
716c
717      else
718        t = 0.0d00
719        t2 = tmax
720c        thi = c + t2*(b + t2*(a+t2*(vfdot+t2*0.25d00*f2)))
721        thi = (-c-t2**2*(a+t2*(vfdot+t2*0.25d00*f2)))/b - t2
722        do i = 1, 1000
723          if (t.eq.0.0d00) then
724            t1= tmax - tmax*dble(i)/1000.0d00
725c            tlo = c + t1*(b + t1*(a+t1*(vfdot+t1*0.25d00*f2)))
726            tlo = (-c-t1**2*(a+t1*(vfdot+t1*0.25d00*f2)))/b - t1
727            if ((thi.ge.0.0d00.and.tlo.le.0.0d00).or.
728     +          (thi.le.0.0d00.and.tlo.ge.0.0d00)) then
729              t = t1
730              if (debug) write(6,*) me, 'tlo: ',tlo
731              if (debug) write(6,*) me, 'thi: ',thi
732              if (debug) write(6,*) me, 't: ',t
733              go to 500
734            endif
735          endif
736        end do
737      endif
738c
739c  Use bisections to find accurate solution
740c
741  500 if (forward) then
742        t1 = 0.0d00
743        t2 = t
744      else
745        t1 = t
746        t2 = tmax
747        if (debug) write(6,*) me, 't1: ',t1
748        if (debug) write(6,*) me, 't2: ',t2
749      endif
750      iter = 0
751c  600 tlo = c + t1*(b + t1*(a+t1*(vfdot+t1*0.25d00*f2)))
752c      thi = c + t2*(b + t2*(a+t2*(vfdot+t2*0.25d00*f2)))
753  600 tlo = (-c-t1**2*(a+t1*(vfdot+t1*0.25d00*f2)))/b - t1
754      thi = (-c-t2**2*(a+t2*(vfdot+t2*0.25d00*f2)))/b - t2
755      if ((thi.ge.0.0d00.and.tlo.le.0.0d00).or.
756     +    (thi.le.0.0d00.and.tlo.ge.0.0d00)) then
757        t = 0.5d00*(t1+t2)
758c        tmid = c + t*(b + t*(a+t*(vfdot+t*0.25d00*f2)))
759        tmid = (-c-t**2*(a+t*(vfdot+t*0.25d00*f2)))/b - t
760        if ((tmid.ge.0.0d00.and.thi.le.0.0d00).or.
761     +      (tmid.le.0.0d00.and.thi.ge.0.0d00)) then
762            t1 = t
763        else if ((tmid.ge.0.0d00.and.tlo.le.0.0d00).or.
764     +           (tmid.le.0.0d00.and.tlo.ge.0.0d00)) then
765            t2 = t
766        endif
767        tcut = 0.5d00*(t1+t2)
768        iter = iter + 1
769        if (abs(t1-t2)/tau.gt.1.0d-14.and.iter.lt.1000) go to 600
770c
771c  Protect against numerical roundoff errors
772c
773        if (iter.ge.1000) then
774          write(6,*) ga_nodeid(),' Maxed out on iterations'
775        endif
776        if (forward.and.(tcut.lt.1.0d-10.or.
777     +      abs(tcut-tmax).lt.1.0d-10)) then
778          cluster_find_tau = -1.0d00
779          return
780        endif
781      else
782        open(unit=2,file='tau.dat',status='unknown')
783        do j = 0, 1000
784          t = tmax*dble(j)/1000.0d00
785          t1 = (-c-t**2*(a+t*(vfdot+t*0.25d00*f2)))/b - t
786          tmid = c + t*(b + t*(a+t*(vfdot+t*0.25d00*f2)))
787          write(2,300) t,t1, tmid
788  300     format(3('    ',f16.8))
789        end do
790        close(2)
791        if (iat.ne.isav) then
792          if (abs(thi).lt.0.000001d00) then
793            tcut = tmax
794          else if (abs(tlo).lt.0.000001d00) then
795            if (vrdot.gt.0.0d00) then
796              tcut = 0.0d00
797            else
798              tcut = tmax
799            endif
800          else
801            if (forward) write(6,*) 'Searching forward'
802            write(6,*) 'Collision time does not converge at step ',istep
803            tcut = -1.0d00
804          endif
805        else
806          tcut = 4.0d00*tau
807        endif
808      endif
809c
810      if (forward) write(6,*) ga_nodeid(),' tcut = ',tcut
811      cluster_find_tau = tcut
812      return
813      end
814c
815      double precision function cluster_check_radius()
816#include "common.fh"
817      double precision rx, ry, rz, r2, cl2
818      integer i
819c
820c   check to make sure that all cluster particles are within cluster
821c   radius
822c
823      cl2 = 0.0d00
824      do i = 1, antot
825        if (at(i).eq.2) then
826          rx = ra(i,1,6) - cl_cmx
827          ry = ra(i,2,6) - cl_cmy
828          rz = ra(i,3,6) - cl_cmz
829          r2 = rx**2 + ry**2 + rz**2
830          if (r2.gt.cl2) cl2 = r2
831        endif
832      end do
833      call ga_dgop(4,cl2,1,'max')
834c
835      cluster_check_radius = sqrt(cl2)
836      return
837      end
838c
839      subroutine cluster_mc
840#include "common.fh"
841c
842c  Perform Monte Carlo adjustment of volume on confining volume
843c
844      double precision r_old, delta_r, r2, cl2, rx, ry, rz, pi
845      double precision vol_new, vol_old, x, ran1, ratio
846      integer i, icnt, iat, me
847      logical force_move,debug
848c
849      if (nocluster) return
850      if (istep.ge.6932366) then
851        debug = .false.
852      else
853        debug = .false.
854      endif
855      call cluster_com
856      me = ga_nodeid()
857      r_cluster_old = r_cluster
858      r_old = r_cluster
859      force_move = .false.
860      delta_r = 0.2
861      if (debug) write(6,*)me,' (1) r_cluster ',r_cluster,istep
862      if (me.eq.0) then
863        r_cluster = r_cluster + delta_r*(ran1(0)-0.5d00)
864      else
865        r_cluster = 0.0d00
866      endif
867      if (debug) write(6,*)me,' (2) r_cluster ',r_cluster,istep
868      call ga_dgop(1,r_cluster,1,'+')
869      if (debug) write(6,*)me,' (3) r_cluster ',r_cluster,istep
870c
871c  If new value of r_cluster falls outside of allowed interval then
872c  it is rejected, unless old value was aready outside interval.
873c
874      if (r_cluster.ge.cl_upper.or.r_cluster.le.cl_lower) then
875c
876c  If new values move cluster radius closer to cutoffs from the outside,
877c  then accept them, reject them otherwise. Before accepting move, make
878c  sure that there are no illegal overlaps.
879c
880        force_move = .true.
881        if (r_cluster.gt.r_old.and.r_cluster.gt.cl_upper) then
882          r_cluster = r_old
883        else if (r_cluster.lt.r_old.and.r_cluster.lt.cl_lower) then
884          r_cluster = r_old
885        endif
886        if (r_cluster.eq.r_old) return
887      endif
888      if (debug) write(6,*)me,' (4) r_cluster ',r_cluster,istep
889c
890c  Check to see if any particles are outside new value of r_cluster.
891c  If so, then new value is rejected.
892c
893      if (r_cluster.lt.r_old) then
894        icnt = 0
895        cl2 = r_cluster**2
896        do i = 1, antot
897          if (at(i).eq.2) then
898            rx = ra(i,1,6) - cl_cmx
899            ry = ra(i,2,6) - cl_cmy
900            rz = ra(i,3,6) - cl_cmz
901            r2 = rx**2 + ry**2 + rz**2
902            if (r2.gt.cl2) then
903              icnt = icnt + 1
904            endif
905          endif
906        end do
907        call ga_igop(9,icnt,1,'+')
908        if (icnt.gt.0) then
909          r_cluster = r_old
910          return
911        endif
912      endif
913      if (debug) write(6,*)me,' (5) r_cluster ',r_cluster,istep
914      if (force_move) return
915c
916c  Accept with Monte Carlo probability.
917c
918      pi = 4.0d00*atan(1.0d00)
919      vol_new = 4.0d00*pi*r_cluster**3/3.0d00
920      vol_old = 4.0d00*pi*r_old**3/3.0d00
921      ratio = (r_cluster/r_old)**2
922c
923      x = ratio*exp(-cl_prssr*(vol_new-vol_old)/mc_tmprtr)
924      if (me.eq.0) then
925        if (x.ge.1.0d00) then
926          icnt = 1
927        else
928          if (ran1(0).lt.x) then
929            icnt = 1
930          else
931            icnt = 0
932          endif
933        endif
934      else
935        icnt = 0
936      endif
937      call ga_igop(1,icnt,1,'+')
938      if (icnt.eq.0) then
939        r_cluster = r_old
940      endif
941      if (debug) write(6,*)me,' (6) r_cluster ',r_cluster,istep
942      return
943      end
944c
945      subroutine cluster_binr
946#include "common.fh"
947c
948c  Bin the current value of r_cluster
949c
950      double precision dr
951      integer ir, isample
952      if (nocluster) return
953c
954c  Find bin
955c
956      dr = r_cluster - cl_lower
957      ir = int(dr/mc_step) + 1
958      if (ir.le.0) ir = 1
959      if (ir.gt.mcbins) ir = mcbins
960      mc_cnt = mc_cnt + 1
961c
962c  Find slice for statistics
963c
964      if (istep.gt.mc_start) then
965        isample = int(10.0d00*(dble(istep-mc_start)-0.5d00)
966     +          / dble(nstep-mc_start)) + 1
967        if (isample.gt.10) isample = 10
968      else
969        isample = 1
970      endif
971      r_cnt(isample) = r_cnt(isample) + 1
972      r_distr(ir,isample) = r_distr(ir,isample) + 1
973      return
974      end
975c
976      subroutine cluster_print_binr
977#include "common.fh"
978c
979c  Print out the distribution of r_cluster
980c
981      double precision rdist(MAXBINS), norm, r, normi
982      double precision rdist2(MAXBINS), sum
983      double precision sigr(MAXBINS),sigr2(MAXBINS)
984      double precision sig2(MAXBINS),sig2_2(MAXBINS)
985      double precision rdisti(MAXBINS,10),rdisti2(MAXBINS,10)
986      integer itot, i, j, isample
987      character*32 filename
988c
989      if (nocluster) return
990      if (task_id.lt.10) then
991        write(filename,100) task_id
992      else if (task_id.ge.10.and.task_id.lt.100) then
993        write(filename,101) task_id
994      else if (task_id.ge.100.and.task_id.lt.1000) then
995        write(filename,102) task_id
996      else if (task_id.ge.1000.and.task_id.lt.10000) then
997        write(filename,103) task_id
998      endif
999  100 format('cl.distr.',i1)
1000  101 format('cl.distr.',i2)
1001  102 format('cl.distr.',i3)
1002  103 format('cl.distr.',i4)
1003c
1004c  Evaluate normalized distribution and uncertainty for P(r)
1005c
1006      if (mc_cnt.lt.0) return
1007      norm = dble(mc_cnt)
1008      do i = 1, mcbins
1009        rdist(i) = 0.0d00
1010        rdist2(i) = 0.0d00
1011        sig2(i) = 0.0d00
1012        sig2_2(i) = 0.0d00
1013        sigr(i) = 0.0d00
1014        sigr2(i) = 0.0d00
1015        do j = 1, 10
1016          rdisti(i,j) = 0.0d00
1017          rdisti2(i,j) = 0.0d00
1018        end do
1019      end do
1020      do j = 1, 10
1021        normi = dble(r_cnt(j))
1022        do i = 1, mcbins
1023          rdist(i) = rdist(i) + dble(r_distr(i,j))/norm
1024          rdisti(i,j) = dble(r_distr(i,j))/normi
1025        end do
1026      end do
1027      do j = 1, 10
1028        do i = 1, mcbins
1029          sig2(i) = sig2(i) + rdisti(i,j)**2
1030        end do
1031      end do
1032      do i = 1, mcbins
1033        sigr(i) = (sig2(i)-10.0d00*rdist(i)**2)/9.0d00
1034      end do
1035c
1036c  Uncertainty at the 95% confidence level
1037c
1038      do i = 1, mcbins
1039        sigr(i) = 2.23d00*sqrt(sigr(i)/10.d00)
1040      end do
1041      sum = 0.0d00
1042      do i = 1, mcbins
1043        sum = sum + rdist(i)
1044      end do
1045      do i = 1, mcbins
1046        rdist(i) = rdist(i)/(sum*mc_step)
1047        sigr(i) = sigr(i)/(sum*mc_step)
1048      end do
1049c
1050c  Evaluate normalized distribution and uncertainty for P(V)
1051c
1052      do j = 1, 10
1053        normi = dble(r_cnt(j))
1054        do i = 1, mcbins
1055          r = cl_lower + (dble(i)-0.5d00)*mc_step
1056          rdist2(i) = rdist2(i) + dble(r_distr(i,j))/(r**2*norm)
1057          rdisti2(i,j) = dble(r_distr(i,j))/(r**2*normi)
1058        end do
1059      end do
1060      do j = 1, 10
1061        do i = 1, mcbins
1062          sig2_2(i) = sig2_2(i) + rdisti2(i,j)**2
1063        end do
1064      end do
1065      do i = 1, mcbins
1066        sigr2(i) = (sig2_2(i)-10.0d00*rdist2(i)**2)/9.0d00
1067      end do
1068c
1069c  Uncertainty at the 95% confidence level
1070c
1071      do i = 1, mcbins
1072        sigr2(i) = 2.23d00*sqrt(sigr2(i)/10.d00)
1073      end do
1074      sum = 0.0d00
1075      do i = 1, mcbins
1076        sum = sum + rdist2(i)
1077      end do
1078      do i = 1, mcbins
1079        rdist2(i) = rdist2(i)/(sum*mc_step)
1080        sigr2(i) = sigr2(i)/(sum*mc_step)
1081      end do
1082c
1083      if (ga_nodeid().eq.0) then
1084        open(unit=2,file=filename,status='unknown')
1085        do i = 1, mcbins
1086          r = cl_lower + (dble(i)-0.5d00)*mc_step
1087          write(2,200) r, rdist(i), rdist2(i),sigr(i),sigr2(i)
1088  200     format(f12.4,'    ',f16.8,'    ',f16.8,'       ',
1089     +           f16.8,'       'f16.8)
1090        end do
1091        close(2)
1092      endif
1093c
1094      return
1095      end
1096c
1097      subroutine cluster_clear_binr
1098#include "common.fh"
1099c
1100c  clear the distribution of r_cluster
1101c
1102      integer itot, i, j
1103c
1104      if (nocluster) return
1105      do j = 1, 10
1106        do i = 1, MAXBINS
1107          r_distr(i,j) = 0
1108        end do
1109        r_cnt(j) = 0
1110      end do
1111      mc_cnt = 0
1112      mc_start = istep
1113c
1114      return
1115      end
1116c
1117      subroutine cluster_reset_binr(iflg)
1118#include "common.fh"
1119c
1120c   recalculate bounds for confining sphere
1121c
1122      integer iflg,ilow,ihi,imax,i
1123      if (nocluster) return
1124      if (iflg.eq.1) then
1125c
1126c   find minimum and maximum values of distribution and reset
1127c   windows to just bound these values
1128c
1129        if (ga_nodeid().eq.0.and.l_rad) then
1130          open(unit=3,file="win1.dat",status='unknown')
1131          do i = 1, mcbins
1132            write(3,*) i, r_distr(i,1)
1133          end do
1134        endif
1135        ilow = 0
1136        ihi = 0
1137        do i = 1, mcbins
1138          if (r_distr(i,1).gt.0.and.ilow.eq.0) then
1139            ilow = i-1
1140            if (ilow.lt.1) ilow = 1
1141          endif
1142          if (ilow.gt.0.and.ihi.eq.0.and.r_distr(i,1).eq.0) then
1143            ihi = i
1144          endif
1145        end do
1146        if (ihi.eq.0) ihi = mcbins
1147        if (ihi.lt.mcbins) then
1148          cl_upper = cl_lower + mc_step*dble(ihi)
1149        endif
1150        if (ihi.gt.1) then
1151          cl_lower = cl_lower + mc_step*dble(ilow-1)
1152        endif
1153      else if (iflg.eq.2) then
1154c
1155c   reset bounds so that maximum is set at radius value that
1156c   is the maximum value of distribution
1157c
1158        ilow = 1
1159        ihi = 0
1160        imax = 0
1161        do i = 1, mcbins
1162c          if (r_distr(i,1).gt.0.and.ilow.eq.0) then
1163c            ilow = i-1
1164c            if (ilow.lt.1) ilow = 1
1165c          endif
1166          if (r_distr(i,1).gt.imax) then
1167            ihi = i
1168            imax = r_distr(i,1)
1169          endif
1170        end do
1171        if (ihi.eq.0) ihi = mcbins
1172        if (ihi.lt.mcbins) then
1173          cl_upper = cl_lower + mc_step*dble(ihi)
1174        endif
1175        if (ilow.gt.1) then
1176          cl_lower = cl_lower + mc_step*dble(ilow-1)
1177        endif
1178      endif
1179      cl_upper = 3.5d00
1180      cl_lower = 0.5d00
1181      cl_upper = 4.5d00
1182      cl_lower = 0.5d00
1183      mc_step = (cl_upper-cl_lower)/dble(mcbins)
1184      call cluster_clear_binr
1185      return
1186      end
1187c
1188      subroutine cluster_do_cllsn
1189#include "common.fh"
1190      integer iat, i, j, ii, jj, iloc, isav, scndat
1191      double precision r, r2, rc2, rx, ry, rz
1192      double precision rrx, rry, rrz, vvx, vvy, vvz
1193      double precision vx, vy, vz, fx, fy, fz, mu, v2, f2
1194      double precision comm(4), tmax, tcut, vrdot
1195      double precision rnx, rny, rnz, vpx, vpy, vpz, vllx, vlly, vllz
1196c
1197c  Re-evaluate all particle velocities. First, recalculate center of
1198c  mass and center of mass velocity of cluster at new coordinates
1199c  and velocities.
1200c
1201      if (nocluster) return
1202      cllsn_cnt = cllsn_cnt + 1
1203      iat = cllsn_idx
1204      call cluster_com
1205      if (iat.gt.0) then
1206        rrx = (cl_mass*cl_cmx - mass(iat)*ra(iat,1,6))/mmass
1207        rry = (cl_mass*cl_cmy - mass(iat)*ra(iat,2,6))/mmass
1208        rrz = (cl_mass*cl_cmz - mass(iat)*ra(iat,3,6))/mmass
1209        vvx = (cl_mass*cl_vcmx - mass(iat)*ra(iat,1,2))/mmass
1210        vvy = (cl_mass*cl_vcmy - mass(iat)*ra(iat,2,2))/mmass
1211        vvz = (cl_mass*cl_vcmz - mass(iat)*ra(iat,3,2))/mmass
1212c
1213c  The r and v vectors are the relative coordinates and velocities of
1214c  particle iat and the center of mass vectors. The vector f is the
1215c  force along this relative coordinate.
1216c
1217        rx = ra(iat,1,6) - rrx
1218        ry = ra(iat,2,6) - rry
1219        rz = ra(iat,3,6) - rrz
1220        vx = ra(iat,1,2) - vvx
1221        vy = ra(iat,2,2) - vvy
1222        vz = ra(iat,3,2) - vvz
1223c
1224c  Decompose velocity v into components that are parallel and perpendicular
1225c  to r. Save this information for computing post-collision velocities.
1226c
1227        r = sqrt(rx**2 + ry**2 + rz**2)
1228        rnx = rx/r
1229        rny = ry/r
1230        rnz = rz/r
1231        vrdot = vx*rnx + vy*rny + vz*rnz
1232        vllx = vrdot * rnx
1233        vlly = vrdot * rny
1234        vllz = vrdot * rnz
1235        vpx = vx - vllx
1236        vpy = vy - vlly
1237        vpz = vz - vllz
1238        comm(1) = (vpx-vllx-vx)/cl_mass
1239        comm(2) = (vpy-vlly-vy)/cl_mass
1240        comm(3) = (vpz-vllz-vz)/cl_mass
1241        comm(4) = mmass
1242      else
1243        comm(1) = 0.0d00
1244        comm(2) = 0.0d00
1245        comm(3) = 0.0d00
1246        comm(4) = 0.0d00
1247      endif
1248      call ga_dgop(8,comm,4,'+')
1249      vvx = comm(1)
1250      vvy = comm(2)
1251      vvz = comm(3)
1252      mmass = comm(4)
1253c
1254      do i = 1, cl_tot
1255        ii = cl_at(i)
1256        if (ii.eq.iat) then
1257          ra(ii,1,2) = ra(ii,1,2) + mmass*vvx
1258          ra(ii,2,2) = ra(ii,2,2) + mmass*vvy
1259          ra(ii,3,2) = ra(ii,3,2) + mmass*vvz
1260        else
1261          ra(ii,1,2) = ra(ii,1,2) - mass(ii)*vvx
1262          ra(ii,2,2) = ra(ii,2,2) - mass(ii)*vvy
1263          ra(ii,3,2) = ra(ii,3,2) - mass(ii)*vvz
1264        endif
1265      end do
1266      call cluster_old_at
1267c
1268      return
1269      end
1270c
1271      subroutine cluster_reset(dt)
1272#include "common.fh"
1273      double precision dt, htausq, scale
1274      integer i, ii
1275c
1276c  Reset coordinates and velocities using old values from beginning of
1277c  step.
1278c
1279      htausq = 0.5d00*dt**2
1280      do i = 1, cl_tot
1281        ii = cl_at(i)
1282        ra(ii,1,1) = cl_old(i,1,1)+dt*cl_old(i,1,3)+htausq*ra(ii,1,3)
1283        ra(ii,2,1) = cl_old(i,2,1)+dt*cl_old(i,2,3)+htausq*ra(ii,2,3)
1284        ra(ii,3,1) = cl_old(i,3,1)+dt*cl_old(i,3,3)+htausq*ra(ii,3,3)
1285        ra(ii,1,6) = cl_old(i,1,2)+dt*cl_old(i,1,3)+htausq*ra(ii,1,3)
1286        ra(ii,2,6) = cl_old(i,2,2)+dt*cl_old(i,2,3)+htausq*ra(ii,2,3)
1287        ra(ii,3,6) = cl_old(i,3,2)+dt*cl_old(i,3,3)+htausq*ra(ii,3,3)
1288        ra(ii,1,2) = cl_old(i,1,3)+dt*ra(ii,1,3)
1289        ra(ii,2,2) = cl_old(i,2,3)+dt*ra(ii,2,3)
1290        ra(ii,3,2) = cl_old(i,3,3)+dt*ra(ii,3,3)
1291      end do
1292#if 1
1293      do i = 1, sl_tot
1294        ii = sl_at(i)
1295        ra(ii,1,1) = sl_old(i,1,1)+dt*sl_old(i,1,3)+htausq*ra(ii,1,3)
1296        ra(ii,2,1) = sl_old(i,2,1)+dt*sl_old(i,2,3)+htausq*ra(ii,2,3)
1297        ra(ii,3,1) = sl_old(i,3,1)+dt*sl_old(i,3,3)+htausq*ra(ii,3,3)
1298        ra(ii,1,6) = sl_old(i,1,2)+dt*sl_old(i,1,3)+htausq*ra(ii,1,3)
1299        ra(ii,2,6) = sl_old(i,2,2)+dt*sl_old(i,2,3)+htausq*ra(ii,2,3)
1300        ra(ii,3,6) = sl_old(i,3,2)+dt*sl_old(i,3,3)+htausq*ra(ii,3,3)
1301        ra(ii,1,2) = sl_old(i,1,3)+dt*ra(ii,1,3)
1302        ra(ii,2,2) = sl_old(i,2,3)+dt*ra(ii,2,3)
1303        ra(ii,3,2) = sl_old(i,3,3)+dt*ra(ii,3,3)
1304      end do
1305#endif
1306#if 1
1307c
1308c   resete extended Hamiltonian parameters
1309c
1310      if ((prsflg.or.ptflg).and.ipmode.eq.0) then
1311        vol1 = cl_vol1_old + dt * cl_vol2_old + htausq * vol3
1312        vol2 = cl_vol2_old + dt * vol3
1313        scale = vol1 / (cl_box_old(1) * cl_box_old(2) * cl_box_old(3))
1314        scale = exp(log(scale)/3.0d00)
1315        xbox = scale * cl_box_old(1)
1316        ybox = scale * cl_box_old(2)
1317        zbox = scale * cl_box_old(3)
1318        xbox2 = 0.5d00 * xbox
1319        ybox2 = 0.5d00 * ybox
1320        zbox2 = 0.5d00 * zbox
1321      endif
1322c
1323      if ((prsflg.or.ptflg).and.ipmode.eq.1) then
1324        do 500 i = 1, 3
1325          alen1(i) = cl_alen1_old(i) + dt * cl_alen2_old(i)
1326     +             + htausq * alen3(i)
1327          alen2(i) = cl_alen2_old(i) + dt * alen3(i)
1328  500   continue
1329        xbox = alen1(1)
1330        ybox = alen1(2)
1331        zbox = alen1(3)
1332        xbox2 = 0.5d00 * xbox
1333        ybox2 = 0.5d00 * ybox
1334        zbox2 = 0.5d00 * zbox
1335      endif
1336c
1337      if ((prsflg.or.ptflg).and.ipmode.eq.2) then
1338        do 600 i = 1, 2
1339          alen1(i) = cl_alen1_old(i) + dt * cl_alen2_old(i)
1340     +             + htausq * alen3(i)
1341          alen2(i) = cl_alen2_old(i) + dt * alen3(i)
1342  600   continue
1343        xbox = alen1(1)
1344        ybox = alen1(2)
1345        xbox2 = 0.5d00 * xbox
1346        ybox2 = 0.5d00 * ybox
1347      endif
1348c
1349c   predictor step for time scale
1350c
1351      if (tmpflg.or.ptflg) then
1352        scal1 = cl_scal1_old + dt * cl_scal2_old + htausq * scal3
1353        scal2 = cl_scal2_old + dt * scal3
1354      endif
1355c
1356c   center of mass parameters
1357c
1358      cl_cmx = cl_cm_old(1)
1359      cl_cmy = cl_cm_old(2)
1360      cl_cmz = cl_cm_old(3)
1361      cl_vcmx = cl_vcm_old(1)
1362      cl_vcmy = cl_vcm_old(2)
1363      cl_vcmz = cl_vcm_old(3)
1364#endif
1365      return
1366      end
1367