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