1 function mc_main(rtdb,grad,thr) 2 implicit none 3 integer rtdb 4 5#include "inp.fh" 6#include "mafdecls.fh" 7#include "rtdb.fh" 8#include "stdio.fh" 9#include "errquit.fh" 10#include "util.fh" 11#include "global.fh" 12#include "geom.fh" 13#include "util_sgroup.fh" 14#include "const_data.fh" 15#include "eaf.fh" 16#include "subgr.fh" 17 logical mc_main 18 logical md_driver,mc_driver,mc_init 19 integer mcsteps 20 logical grad, task_energy 21 external grad, task_energy 22 23 logical status 24 integer natom 25 integer c_start 26 integer vel_init,acc_init 27 logical forward,backward ! tells you which part of the traj is being called 28 logical mc_data_set_ifirc, mc_data_get_ifirc 29 integer mc_data_get_forside, mc_data_get_forside0 30 integer mc_data_get_backside, mc_data_get_backside0 31 integer i,j,k,naccept 32 integer c_array, vel, acc, prp, prp_tmp,sidef,sideb 33 integer mc_data_get_i_c_array, md_data_get_i_v, md_data_get_i_a 34 integer md_data_get_i_c, mc_data_get_i_in_vel,mc_data_get_i_in_acc 35 integer mc_data_get_natom, mc_data_get_i_prp 36 integer mc_data_get_i_prp_tmp 37 integer tag,charge,md_data_get_i_t,md_data_get_i_q 38 character*32 thr 39 logical ircflag,mc_data_set_forward,mc_data_set_backward 40 integer i_sc, mc_data_get_i_sc, nonreact,reactive,reverse 41 integer i_s, mc_data_get_i_s,mc_data_get_mcsteps 42 double precision avg,rc,rc1,temp,mc_data_get_temp 43 logical md_data_set_temp 44 logical mc_data_set_side_prev, mc_data_set_nxing 45 logical mc_data_set_trajnum 46 double precision mdtemp,mdtemp10, md_data_get_temp 47 integer md_data_get_nsteps, nsteps 48 logical dbug 49c<<<<<<< .mine 50 logical restart_dump 51 52c======= 53 character*(nw_max_path_len) fnameout 54 integer tunita 55c 56c>>>>>>> .r18494 57 integer g_a,ncpu 58C PROCESS AND SUBGROUP VARS 59 integer myid,idbig,inodesbig 60 integer idmedium, idzero, g_pr, numgroups 61 integer l_groupnums, k_groupnums 62C I/O vars 63 character*4 prcfil 64 character*4 prcfil2 65 character*256 fprefix, fprefix2 66 character*256 mc_dir 67 integer iw 68 mc_main=.false. 69 iw=6 70 call ga_sync() 71C INFO ON PROCESSES in the big group 72 myid = ga_nodeid() 73 idbig = ga_pgroup_get_default() 74 inodesbig = ga_nnodes() 75 76 if (myid.eq.0) then 77 open(unit=ir,file='restart.pun',status='replace') 78 endif 79 if (.not. rtdb_get(rtdb,'subgroups_number', MT_INT, 1, numgroups)) 80 & numgroups = ga_nnodes() 81 82C ***** Can't have more subgroups then procs ***** 83 if (numgroups.gt.ga_nnodes()) then 84 numgroups = ga_nnodes() 85 write (iw,*) 'Requested number of subgroups too big, 86 & reseting to:', numgroups 87 endif 88C Allocate memory for array of group numbersi, only if more then 1 subgroup requested 89 if (numgroups.gt.1) then 90 if (.not. ma_push_get(MT_INT, numgroups, 'mc: groupnums', 91 & l_groupnums, k_groupnums)) 92 & call errquit('mc_allocate_arrays: error groupnums',0,MA_ERR) 93 94 95 if (ga_nodeid().eq.0) then 96 write(iw,*) 'SPLITTING INTO SUBGROUPS' 97 endif 98 call ga_sync() 99 call util_flush(iw) 100 close(iw) 101C SPLITTING INTO SUBGROUPS 102 103C Close the output file and reopen it as series of ".out" files 104 call mc_setup_ga 105 & (rtdb, myid, idbig, inodesbig, 106 & idmedium, idzero, g_pr, int_mb(k_groupnums),numgroups) ! Now in Subgroups 107 108C Initialize group related information 109C mostly file opening for each group,etc. 110 endif 111 call mc_setup_group(numgroups) !deals with file i/o 112c & (rtdb, myid, int_mb(k_groupnums),prcfil, prcfil2, 113c & fprefix, numgroups) 114 115 write(iw,*) 'DONE SETTING I/O files for subgroups' 116C-- the stdout is now redirected to other .out output files for each trajectory 117C from here we can start mc_loop 118c goto 100 119c status = task_energy(rtdb) 120c goto 100 121 122 write(iw,*) 'testing energy calculation\n' 123 status = task_energy(rtdb) 124 125 write(iw,*) '\n\t\t METROPOLIS MONTE CARLO RECROSSING \n' 126 127 write(iw,*) 'Subgroup is: ',util_sgroup_mygroup() 128 dbug=.false. 129 forward=.true. 130 backward=.false. 131 132C********* MC driver ****** 133 134C******* RESTART point 1****** 135 if(.not.mc_driver(rtdb,naccept)) call errquit('mc_driver:failed', 136 & UNKNOWN_ERR) 137 write(iw,*) 'Number of initial points generated: ', naccept 138 139 if(util_sgroup_mygroup().eq.1.and.ga_nodeid().eq.0) then 140 status=restart_dump() 141 endif 142C******** End of MC driver ******* 143c go to 100 144 ircflag=.true. 145 status=mc_data_set_ifirc(.true.) 146 mdtemp10 =10.0 147 mdtemp = 298.15 148 149 natom=mc_data_get_natom() 150 151 call md_data_set_natom(natom) 152 call md_data_allocate() 153 call md_set(rtdb) ! sets md_temp, nsteps and mc_stepsize 154 mdtemp=md_data_get_temp() 155 tag=md_data_get_i_t() 156 charge=md_data_get_i_q() 157 c_start=md_data_get_i_c() 158 vel=md_data_get_i_v() 159 acc=md_data_get_i_a() 160c 100 continue 161 write(iw,*) '\nStarting the IRC trajectory\n' 162 163c ---- do the IRC trajectory first, at low temp 164C ---- this requires quenching of velocities so that Ek is close to 0 165C ---- quenching not implemented yet 166 status=md_data_set_temp(mdtemp10) 167C --- set vel and acc 168 169 status=mc_init(rtdb,grad,ircflag) 170 171 c_array=mc_data_get_i_c_array() 172 vel_init=mc_data_get_i_in_vel() 173 acc_init=mc_data_get_i_in_acc() 174 i_sc=mc_data_get_i_sc() 175 i_s=mc_data_get_i_s() 176 prp=mc_data_get_i_prp() 177 prp_tmp = mc_data_get_i_prp_tmp() 178 mcsteps=mc_data_get_mcsteps() 179 180 do k=1,3*natom 181 dbl_mb(c_start+k-1)=dbl_mb(c_array+k-1) 182 dbl_mb(vel+k-1)=dbl_mb(vel_init+k-1) 183 dbl_mb(acc+k-1)=dbl_mb(acc_init+k-1) 184 enddo 185 186 if(dbug ) then 187 write(iw,*) 'In MC main coordinates' 188 do i=1,naccept+1 189 write(iw,*) 'point ',i 190 do k=1,3*natom 191 write(iw,*) dbl_mb(c_array+(i-1)*3*natom+k-1) 192 enddo 193 write(iw,*) 'vel ',i 194 do k=1,3*natom 195 write(iw,*) dbl_mb(vel_init+(i-1)*3*natom+k-1) 196 enddo 197 write(iw,*) 'acc ',i 198 do k=1,3*natom 199 write(iw,*) dbl_mb(acc_init+(i-1)*3*natom+k-1) 200 enddo 201 202 enddo 203 endif 204C --- run the IRC traj forward 205 206 status=mc_data_set_forward(.true.) 207 status=mc_data_set_backward(.false.) 208 status=mc_data_set_side_prev(0) 209 status=mc_data_set_nxing(0) 210 status=mc_data_set_trajnum(0) 211 212 if(.not.md_driver(rtdb,grad,thr)) 213 & call errquit('Error in MDDriver',0,0) 214 status=mc_data_set_side_prev(0) 215 status=mc_data_set_nxing(0) 216C --- run the IRC traj backward 217 call md_set(rtdb) 218 do k=1,3*natom 219 dbl_mb(c_start+k-1)=dbl_mb(c_array+k-1) 220 dbl_mb(vel+k-1)=-dbl_mb(vel_init+k-1) 221 dbl_mb(acc+k-1)=-dbl_mb(acc_init+k-1) 222 enddo 223 224 status=mc_data_set_forward(.false.) 225 status=mc_data_set_backward(.true.) 226 227 if(.not.md_driver(rtdb,grad,thr)) 228 & call errquit('Error in MDDriver',0,0) 229 230 call ga_sync() 231 232C ****** RESTART point 2 ******* 233c ---- dump the data for restart 234 if(util_sgroup_mygroup().eq.1.and.ga_nodeid().eq.0) then 235 status=restart_dump() 236 endif 237 238 write(iw,*) 'Done with IRC' 239 ircflag=.false. 240 status=mc_data_set_ifirc(.false.) 241C --- generate initial velocities for all trajectories 242 status=mc_init(rtdb,grad,ircflag) 243 244 do i=1,naccept+1 245 write(iw,*) 'traj, prp: ',i,dbl_mb(prp+i-1) 246 dbl_mb(prp_tmp+i-1) = dbl_mb(prp+i-1) 247 enddo 248 249 write(iw,*) 'STARTING MD TRAJECTORIES ' 250 251 status=mc_data_set_forward(.true.) 252 status=mc_data_set_backward(.false.) 253 254C --- START LOOP OVER MD TRAJ - PARALLEL 255 call ga_sync() 256 ncpu=numgroups 257 258C ****** RESTART point 3 ******* 259 do i=1,naccept+1 260C ---- decide where to run the traj 261 if(util_sgroup_mygroup().eq.mod(i,ncpu)+1) then 262 263C ---- set the appropriate coordinate 264 call md_set(rtdb) 265 do k=1,3*natom 266 dbl_mb(c_start+k-1)=dbl_mb(c_array+(i-1)*3*natom+k-1) 267 dbl_mb(vel+k-1)=dbl_mb(vel_init+(i-1)*3*natom+k-1) 268 dbl_mb(acc+k-1)=dbl_mb(acc_init+(i-1)*3*natom+k-1) 269 enddo 270 if (dbug ) then 271 write(iw,*) 'Coordinates:' 272 do k=1,3*natom 273 write(iw,*) dbl_mb(c_start+k-1) 274 enddo 275 write(iw,*) 'Velocieties:' 276 do k=1,3*natom 277 write(iw,*) dbl_mb(vel+k-1) 278 enddo 279 write(iw,*) 'Acceleration:' 280 do k=1,3*natom 281 write(iw,*) dbl_mb(acc+k-1) 282 enddo 283 endif 284 do j=1,3 285 if (forward.and. .not.backward) then 286 write(iw,*) 'Forward: traj',i 287 288 status=mc_data_set_side_prev(0) 289 status=mc_data_set_nxing(0) 290 status=mc_data_set_trajnum(i) 291 292 if(.not.md_driver(rtdb,grad,thr)) 293 & call errquit('Error in MDDriver',0,0) 294 forward=.false. 295 backward=.true. 296 status=mc_data_set_forward(forward) 297 status=mc_data_set_backward(backward) 298 status=mc_data_set_side_prev(0) 299 status=mc_data_set_nxing(0) 300 301 elseif(backward.and. .not.forward) then 302 write(iw,*) 'Backward: traj',i 303C --- reset the geometry to c_start 304C --- get the initial conditions generated in the forward run 305 call md_set(rtdb) 306 do k=1,3*natom 307 dbl_mb(c_start+k-1)=dbl_mb(c_array+(i-1)*3*natom+k-1) 308 dbl_mb(vel+k-1)=-dbl_mb(vel_init+(i-1)*3*natom+k-1) 309 dbl_mb(acc+k-1)=-dbl_mb(acc_init+(i-1)*3*natom+k-1) 310 enddo 311 status=mc_data_set_side_prev(0) 312 status=mc_data_set_nxing(0) 313 314 if(.not.md_driver(rtdb,grad,thr)) 315 & call errquit('Error in MDDriver',0,0) 316 forward=.false. 317 backward=.false. 318 status=mc_data_set_side_prev(0) 319 status=mc_data_set_nxing(0) 320 321 write(iw,9000) 'traj ifirc forside0 backside0 forside 322 & backside', 323 & i, 324 & mc_data_get_ifirc(), 325 & mc_data_get_forside0(),mc_data_get_backside0(), 326 & mc_data_get_forside(),mc_data_get_ backside() 327 328 status=mc_data_set_forward(forward) 329 status=mc_data_set_backward(backward) 330 else 331 forward=.true. 332 backward=.false. 333 status=mc_data_set_forward(forward) 334 status=mc_data_set_backward(backward) 335 336C --- frajectory finished, moved to the other trajectory 337C --- here you would also check if the trajectory was reactive or not 338 endif 339 enddo 340c else 341c write(*,*) 'SKIP TRAJECTORY ',i,' node ',ga_nodeid() 342 endif 343c write(*,*) 'LOOP ',i,' node ',ga_nodeid() 344c call ga_sync() 345c idmedium = ga_pgroup_get_default() 346c write(iw,*) 'big,med,node',idbig,idmedium,ga_nodeid() 347c flush(iw) 348c call ga_pgroup_set_default(idbig) 349* CCCCCC igop needs to sync all procs, so it makes the subgroups wait for eachother 350* CCCCCc may need to have it as a global array 351c call ga_igop(MT_F_INT,int_mb(i_sc),(mcsteps+1)*2,'max') 352c call ga_pgroup_set_default(idmedium) 353c if(util_sgroup_mygroup().eq.1.and.ga_nodeid().eq.0) then 354c status=restart_dump() 355c endif 356 enddo 357c write(*,*) ' Out of loop. Node id :',ga_nodeid() 358C --- loop over trajectories finished 359 call ga_sync() 360 call ga_pgroup_set_default(idbig) 361 call ga_sync() 362 nsteps=md_data_get_nsteps() 363 call ga_igop(MT_F_INT,int_mb(i_sc),(mcsteps+1)*2*nsteps,'max') 364 call ga_sync() 365 status=restart_dump() 366 write(iw,*) ' Recrossing Results:' 367 368 369 370 if (ga_nodeid().eq.0) then 371c loop over each mdstep to get kappa for each step 372 rc1=0.0 373 do k=0,nsteps-1 374 nonreact=0 375 reactive=0 376 reverse=0 377 do i=1,naccept+1 378c int_mb(i_sc+(trajnum-1)*2*nsteps+mdstep-1) = side 379 sidef = int_mb(i_sc+2*(i-1)*nsteps+k) 380 sideb = int_mb(i_sc+2*(i-1)*nsteps+nsteps+k) 381 if(dbug) 382 & write(iw,9001) 'step traj sidef sideb ',k,i,sidef,sideb 383 if (sidef.eq.sideb) then 384 nonreact=nonreact+1 385 dbl_mb(prp_tmp+i-1)=0.0 386c if (dbug) 387 write(iw,*) 'tr number for non-reactive traj: ', i 388 else if (sidef.eq.mc_data_get_forside0()) then 389 dbl_mb(prp_tmp+i-1)=0.5*dbl_mb(prp_tmp+i-1) 390c reactive=reactive+1 391 else 392c reverse=reverse+1 393 dbl_mb(prp_tmp+i-1)=-0.5*dbl_mb(prp_tmp+i-1) 394 endif 395 if (dbl_mb(prp_tmp+i-1).gt.0.0) then 396 reactive=reactive+1 397 else if (dbl_mb(prp_tmp+i-1).lt.0.0) then 398 reverse=reverse+1 399 endif 400 enddo 401 402 j=0 403 avg=0 404 405 406 do i=0,mcsteps 407 if (int_mb(i_s+i).eq.1) then 408 j=j+1 409 endif 410c if(dbug) then 411 write(iw,*) 'Mc step: ',i,'adding: ',j,' val: ', 412 & dbl_mb(prp_tmp+j-1) 413c endif 414 avg=avg+dbl_mb(prp_tmp+j-1) 415 enddo 416 417 avg=avg/(mcsteps+1) 418 temp = mc_data_get_temp() 419 rc=avg 420c rc =avg/sqrt((temp*boltz)/(twopi)) 421 422 if(k.eq.0) rc1=rc 423 424 write(iw,9001) 'non-reactive, reactive,reverse traj: ', 425 & nonreact,reactive,reverse 426 write(iw,*) ' Average prp value: ',avg 427 write(iw,*) 'Analytical avg: ',sqrt((temp*boltz)/(twopi)) 428 write(iw,9002) 'MD Step, Recrossing coef: ',k+1,rc/rc1 429 write(iw,9002) 'MD Step, Recrossing coef anl: ',k+1, 430 $ rc/sqrt((temp*boltz)/(twopi)) 431C --- reset the original prp values 432 do i=1,naccept+1 433 dbl_mb(prp_tmp+i-1) = dbl_mb(prp+i-1) 434 enddo 435 enddo ! end of loop over nsteps 436 endif 437c 100 continue 438 call ga_pgroup_set_default(idbig) 439 call ga_sync() 440 call mc_data_free_all() 441 call md_data_free_all() 442c --- close files that were open 443 close(iw) 444 close(12) 445 close(ir) 446 call util_sgend(rtdb) 447 mc_main = .true. 448 9000 format(A,I6,L3,4(I6)) 449 9001 format(A,3(I6)) 450 9002 format(A,I6,F10.6) 451 return 452 end 453 454C*********************************************************************** 455 subroutine mc_setup_ga 456 & (rtdb, myid, idbig, inodesbig, 457 & idmedium, idzero, g_pr, groupnums, numgroups) 458C*********************************************************************** 459 Implicit none 460C Include Files 461#include "rtdb.fh" 462#include "errquit.fh" 463#include "mafdecls.fh" 464#include "tcgmsg.fh" 465#include "global.fh" 466#include "msgids.fh" 467#include "msgtypesf.h" 468#include "util.fh" 469#include "util_sgroup.fh" 470#include "subgr.fh" 471#include "pstat.fh" 472 473C Variable Declarations 474 integer rtdb ! input 475C integer nob, nspc 476 integer myid, idbig, inodesbig, numgroups ! input 477C Depreciated 478 logical status ! internal use 479C Depreciated 480 integer idmedium 481 integer ld(2)! internal use 482 integer g_proc, procnums, groupnums, bigproc !p_proc is internal 483C procnums is internal 484 integer i, j, idzero ! i and j are internal indexes 485 integer ndim, dims(2), chunk(2), g_pr !n(dim)s and chunck are internal 486 external util_sggo 487C Additional variables 488 integer groups_want, array_cpu(1), method 489 integer dir 490C End Additions 491 492C Dimensions 493 dimension procnums(numgroups) 494 dimension groupnums(numgroups) 495 dimension bigproc(inodesbig,1) 496 497C timers common 498C Create big GAs 499 if (.not.ga_create(MT_DBL, inodesbig, 1, "proc list", 500 & 1, -1, g_proc)) 501 & call errquit('dntmc_setup_ga:g_proc create error', 0, 502 & GA_ERR) 503 call ga_fill(g_proc, -1) 504 505C Create Subgroups 506 if (.not.rtdb_get(rtdb, 'subgroups_number', mt_int, 1, 507 & groups_want))groups_want=ga_nnodes() 508C write(*,*) 'groups_want =',groups_want 509C groups_want is the number of subgroups, so initialy we want each CPU to be subgroup 510C that is all done in the input with set "subgroups_number" 511 512C Setting Method = 1 513c 1 -- use groups_want to generate equal sized groups (array_cpu ignored) 514C Simplest. 515c 2 -- turn each SMP box into a group (array_cpu and groups_want ignored) 516C This uses GA to tell it about the cluster. 517c 3 -- use array_cpu(groups_want) to define number of nodes per group 518c 4 -- use array_cpu(groups_want+nnodes) to define which nodes per group 519C This is just option 3, but you get to lay the groups out exactly. 520 method = 1 521 array_cpu(1) = 0 522 dir = 1 ! Write group rtdb's in scratch directories 523 524 call util_sggo(rtdb,groups_want,method,array_cpu,dir) 525 526 if (util_print('debug',print_debug)) then 527 write(0,*)'Now in Subgroups' 528 call flush(6) 529 endif 530 531C I this is the handle to get into subgroups subgroup 532 idmedium = ga_pgroup_get_default() ! Now in Subgroups 533C This is the handele to a group of 0th cpus 534C for example: 6 cpus in 3 groups, will have zerogroup(0,2,4) 535 idzero = util_sgroup_zero_group() 536C write(6,*) 'idbig,idmedium, idzero:',idbig,idmedium,idzero 537 538C Get Back to Big Group 539 call ga_pgroup_set_default(idbig) ! Now in Big Group 540C Create Processor Zero and Group Processor Lists 541 if (ga_pgroup_nodeid(idmedium).eq.0) then 542 ld(1) = 1 !must be physical dimension of local array 543 ld(2) = 1 544 call ga_put(g_proc,myid+1,myid+1,1,1,util_sgroup_mygroup(),ld) 545 call ga_sync() 546 ld(1) = inodesbig !must be physical dimension of local array 547 ld(2) = 1 548 call ga_get(g_proc,1,inodesbig,1,1,bigproc,ld) 549 j = 0 550 do i = 1, inodesbig 551C write(*,*) 'big proc (i,1): ',bigproc(i,1) 552 if (bigproc(i,1).ne.-1) then 553 j = j + 1 554 groupnums(j) = bigproc(i,1) 555 procnums(j) = i-1 556 endif 557 enddo 558 if (j.ne.util_sgroup_numgroups()) 559 & call errquit('mc_setup_ga:zero node creation problem', 560 & j, UNKNOWN_ERR) 561 else 562 call ga_sync() 563 endif 564c write(6,*)'Now broad:',idmedium,msg_dntmc3+MSGINT, procnums 565c write(6,*)'Now broad:',mitob(numgroups) 566C Broadcase Results, array of groups( starting from 0) 567 call ga_pgroup_brdcst(idmedium,msg_dntmc3+MSGINT, procnums, 568 & mitob(numgroups), 0) 569C Broadcase Results, array of groups( starting from 1) 570c write(6,*)'Now broadcasting,again' 571 call ga_pgroup_brdcst(idmedium,msg_dntmc7+MSGINT, groupnums, 572 & mitob(numgroups), 0) 573c write(6,*)'Now in destroying list' 574C Destroy GA "proc list" 575 if (.not. ga_destroy(g_proc)) call errquit('mc_setup_ga: 576 &ga_destroy(g_proc) failed', 0, GA_ERR) 577 578C Create Zeros GAs 579C this may not be needed, or will be needed for setting hared data for a subgroup 580 581c write(6,*)'Now starting subgroups ' 582C Start Subgroups 583 call ga_pgroup_set_default(idmedium) ! Now in Subgroups 584c write(*,*) 'In isubgroups nnodes ',ga_nnodes() 585 call ga_sync() 586C Test write 587c write(*,*) 'Depth: ',depth 588c write(6,*) '@ proc ',myid,' gr ',idmedium,'gr proc',ga_nodeid() 589c write(6,*) '@ proc ',myid,' gr2 ',my_ga_grp(depth), 590c & ' @ with zero group ',util_sgroup_zero_group() 591c if (myid .eq. 0) then 592c write(6,*)'@ group and proc arrays:numgroups',numgroups 593c do i=1,numgroups 594c write(6,*)'@ ',groupnums(i),procnums(i) 595c EndDo 596c endif 597 598 return 599 end 600 601 602C*********************************************************************** 603 subroutine mc_setup_group(numgroups) 604C gets the appropriate path for the output files and opens them 605C*********************************************************************** 606 Implicit none 607 608C Include Statements 609#include "errquit.fh" 610#include "rtdb.fh" 611#include "inp.fh" 612#include "subgr.fh" 613#include "util.fh" 614#include "util_sgroup.fh" 615#include "global.fh" 616C Variable declarations 617c integer rtdb ! input 618c integer myid ! input 619c integer numgroups ! input 620c integer groupnums(numgroups) ! input 621c character*4 prcfil, prcfil2 622c integer i 623C Indect output of file units 624c character*256 fprefix, fprefixcat ! internal use only !fprefix output 625c character*256 fprefix2, mc_dir 626 integer numgroups 627 character*256 xyzname, outname 628 character*256 pxyzname, poutname 629 630C Begin Main Program 631c call mc_write_prcfil(groupnums, prcfil, numgroups) 632 633c call mc_build_prcfil(prcfil2, myid) 634c write(fprefix2,'(256(a))') (' ', i=1,256) 635c if(.not. rtdb_cget(rtdb,'file_prefix',1,fprefix2)) 636c & call errquit('dntmc_setup_group:rtdb get file_prefix failed' 637c & ,0,RTDB_ERR) 638c write(mc_dir,'(256(a))') (' ', i=1,256) 639c if (.not. rtdb_cget(rtdb, 'mc:directory',1,mc_dir)) 640c & mc_dir(1:2)='./' 641c write(fprefix,'(256(a))') (' ', i=1, 256) 642c write(fprefix,'(3(a))') 643c & mc_dir(1:inp_strlen(mc_dir)), 644c & '/', 645c & fprefix2(1:inp_strlen(fprefix2)) 646C Open file Units 647C general output ! Only group Zeros write 648 if (ga_nodeid().eq.0 ) then 649 call util_file_name('xyz', .false., .false., xyzname) 650 call util_file_name('out', .false., .false., outname) 651 call util_pname(xyzname,pxyzname) 652 call util_pname(outname,poutname) 653 write(6,*) 'opening files', pxyzname, poutname 654c write(fprefixcat,'(256(a))') (' ', i=1, 256) 655c write(fprefixcat,'(3(a))') 656c & fprefix(1:inp_strlen(fprefix)), 657c & '.xyz.',prcfil(1:4) 658c write(6,*) 'Filename', xyzname 659c -- open a file for trajectory for each subgroup 660 if(numgroups.gt.1) then 661 OPEN(UNIT = 12,FILE = pxyzname,STATUS = 'REPLACE') 662 663C ---- added this so that the output file is closed and output is open for each 664C ---- subgroup 665 OPEN(UNIT = 6,FILE = poutname,STATUS = 'REPLACE') 666 else 667 OPEN(UNIT = 12,FILE = xyzname,STATUS = 'REPLACE') 668 OPEN(UNIT = 6,FILE = outname,STATUS = 'REPLACE') 669 endif 670 else 671C all non-zeros write to dev/null 672 OPEN(6,FILE='/dev/null',STATUS='UNKNOWN') 673 OPEN(12,FILE='/dev/null',STATUS='UNKNOWN') 674 endif 675 676 677 return 678 end 679 680C*********************************************************************** 681 subroutine mc_write_prcfil(groupnums,prcfil,numgroups) 682C finds the appropriate groupid, to get the right output file 683C not using idmedium, need to remove it 684C*********************************************************************** 685 Implicit none 686 687C Include Statements 688#include "subgr.fh" 689#include "errquit.fh" 690#include "util_sgroup.fh" 691 692C Variable Declarations 693 integer i 694 integer numgroups 695 integer groupnums(numgroups) 696 integer groupid 697 character*4 prcfil ! only output 698 699C Main Program 700 groupid = -1 701 do i=1, util_sgroup_numgroups() 702 if (groupnums(i) .eq. util_sgroup_mygroup()) then 703 groupid = i 704 endif 705 enddo 706 707 if (groupid .eq. -1) 708 & call errquit('dntmc_write_prcfil:failed to allocate idgroup',0, 709 & GA_ERR) 710 711 call mc_build_prcfil(prcfil, groupid) 712 713C End Main Program 714 return 715 end 716 717*********************************************************************** 718 subroutine mc_build_prcfil(prcfil, i) 719C formats the name of the procesor file: ex. 0001,0010,0100 720C so you can have up to 9999 cpus and files 721C*********************************************************************** 722 Implicit none 723 724C Variable Declarations 725 integer i 726 character*4 prcfil ! only output 727 728C Main Program 729 write(prcfil(1:4), '(i4)') i 730 if (i .le. 9) then 731 prcfil(1:3) = '000' 732 endif 733 if (i .le. 99) then 734 prcfil(1:2) = '00' 735 endif 736 if (i .le. 999) then 737 prcfil(1:1) = '0' 738 endif 739 740C End Main Program 741 return 742 end 743 744 function restart_dump() 745 Implicit none 746#include "inp.fh" 747#include "stdio.fh" 748#include "util.fh" 749#include "util_sgroup.fh" 750#include "const_data.fh" 751#include "subgr.fh" 752#include "errquit.fh" 753#include "global.fh" 754#include "mafdecls.fh" 755 logical restart_dump 756C: save the array of initial velocities and initial structures and initial prp values 757C: save the array with 1 and 0 for iniital values 758C: before the run, save values to -1, so at restart we know which ones to get. 759C: restart file has handle ir=11 (set as parameter in const_data.fh 760 integer prp, vel, acc, coor, sidecnt, mcsteps, naccept, natom 761 integer nsteps 762 integer mc_data_get_i_prp !property (naccept+2) 763 integer mc_data_get_i_sc !side count (mcsteps+1)*2*nsteps 764 integer mc_data_get_i_in_vel ! initial velocities (naccept+2)*3*natom 765 integer mc_data_get_i_in_acc ! initial accelerations (naccept+2)*3*natom 766 integer mc_data_get_i_c_array ! initial coordinates (mcsteps+1)*3*natom 767 integer mc_data_get_natom ! number of atoms 768 integer mc_data_get_mcsteps ! mc steps mcsteps+1 769 integer mc_data_get_naccept ! number of accepted steps 770 integer myid, mysubgr 771 integer md_data_get_nsteps ! number of steps in MD trajectory 772 integer i 773C******* double check the size of the array for vel,acc and prp**** 774 restart_dump=.false. 775 natom=mc_data_get_natom() 776 mcsteps=mc_data_get_mcsteps() 777 naccept=mc_data_get_naccept() 778 prp=mc_data_get_i_prp() 779 vel=mc_data_get_i_in_vel() 780 acc=mc_data_get_i_in_acc() 781 coor=mc_data_get_i_c_array() 782 sidecnt=mc_data_get_i_sc() 783 nsteps = md_data_get_nsteps() 784C: Check again if it is subgroup 1, node 0 before you write 785 myid = ga_nodeid() 786 mysubgr = util_sgroup_mygroup() 787c write(ir,*) 'writing to restart file' 788 if (myid.eq.0.and.mysubgr.eq.1) then 789 rewind(ir) 790 write(ir,*) 'natom,mcsteps,naccept,nsteps' 791 write(ir,*) natom,mcsteps,naccept,nsteps 792C figure out a better print 793 write(ir,*) 'structures' 794 call print_array_dbl(dbl_mb(coor),(mcsteps+1)*3*natom,ir) 795 write(ir,*) 'sidecount' 796 call print_array_int(int_mb(sidecnt), (mcsteps+1)*2*nsteps,ir) 797 write(ir,*) 'prp' 798 call print_array_dbl(dbl_mb(prp),(naccept+1),ir) 799 write(ir,*) 'vel' 800 call print_array_dbl(dbl_mb(vel),(naccept+1)*3*natom,ir) 801 write(ir,*) 'acc' 802 call print_array_dbl(dbl_mb(acc),(naccept+1)*3*natom,ir) 803 call util_flush(ir) 804 endif 805 restart_dump=.true. 806 return 807 end 808C ******** End of writing restart file ********* 809C ******** Start reading restart file ********* 810 function restart_read() 811 Implicit none 812#include "inp.fh" 813#include "stdio.fh" 814#include "util.fh" 815#include "util_sgroup.fh" 816#include "const_data.fh" 817#include "subgr.fh" 818#include "errquit.fh" 819#include "global.fh" 820#include "mafdecls.fh" 821 logical restart_read 822C: save the array of initial velocities and initial structures and initial prp values 823C: save the array with 1 and 0 for iniital values 824C: before the run, save values to -1, so at restart we know which ones to get. 825C: restart file has handle ir=11 (set as parameter in const_data.fh 826 integer prp, vel, acc, coor, sidecnt, mcsteps, naccept, natom 827 integer nsteps 828 integer mc_data_get_i_prp !property (naccept+2) 829 integer mc_data_get_i_sc !side count (mcsteps+1)*2*nsteps 830 integer mc_data_get_i_in_vel ! initial velocities (naccept+2)*3*natom 831 integer mc_data_get_i_in_acc ! initial accelerations (naccept+2)*3*natom 832 integer mc_data_get_i_c_array ! initial coordinates (mcsteps+1)*3*natom 833 logical mc_data_set_naccept ! number of accepted steps 834 logical md_data_set_nsteps 835 integer myid, mysubgr 836 integer i 837 logical status 838C******* double check the size of the array for vel,acc and prp**** 839 restart_read=.false. 840 prp=mc_data_get_i_prp() 841 vel=mc_data_get_i_in_vel() 842 acc=mc_data_get_i_in_acc() 843 coor=mc_data_get_i_c_array() 844 sidecnt=mc_data_get_i_sc() 845C: Check again if it is subgroup 1, node 0 before you write 846c write(ir,*) 'writing to restart file' 847 rewind(ir) 848 read(ir,*) natom,mcsteps,naccept,nsteps 849C figure out a better print 850 call read_array_dbl(dbl_mb(coor),(mcsteps+1)*3*natom,ir) 851 call read_array_int(int_mb(sidecnt), (mcsteps+1)*2*nsteps,ir) 852 call read_array_dbl(dbl_mb(prp),(naccept+1),ir) 853 call read_array_dbl(dbl_mb(vel),(naccept+1)*3*natom,ir) 854 call read_array_dbl(dbl_mb(acc),(naccept+1)*3*natom,ir) 855 status=mc_data_set_naccept(naccept) 856 restart_read=.true. 857 return 858 end 859C ******** End reading restart file ********** 860 861 862 863 subroutine print_array_int(array,asize,fileout) 864 implicit none 865 integer asize,fileout 866 integer array(asize) 867 write(fileout,*) array 868 return 869 end 870 subroutine print_array_dbl(array,asize,fileout) 871 implicit none 872 integer asize,fileout 873 double precision array(asize) 874 write(fileout,*) array 875 return 876 end 877 subroutine read_array_int(array,asize,filein) 878 implicit none 879 integer asize,filein 880 integer array(asize) 881 read(filein,*) array 882 return 883 end 884 subroutine read_array_dbl(array,asize,filein) 885 implicit none 886 integer asize,filein 887 double precision array(asize) 888 read(filein,*) array 889 return 890 end 891 892 893c $Id$ 894