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