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