1      subroutine dia_gettrj
2c
3c $Id$
4c
5c     read the trajectory
6c     -------------------
7c
8      implicit none
9c
10#include "dia_common.fh"
11#include "dia_params.fh"
12#include "global.fh"
13#include "mafdecls.fh"
14#include "rtdb.fh"
15#include "msgids.fh"
16#include "util.fh"
17c
18      external dia_rdfram
19      logical dia_rdfram
20c
21      integer nread,nframe,itmp(4),ifil,nfram
22      integer meio,meframf,meframl,mefilf,mefill,meskip,melast
23      integer ifram,iffram,ilfram,lq,ifr
24      real*8 t(6)
25      character*255 fname
26c
27c     deallocate any existing global arrays
28c
29      if(nfrtot.gt.0) then
30      if(.not.ga_destroy(ga_trj))
31     + call md_abort('Failed to destroy ga_trj',0)
32      if(lslvnt) then
33      if(.not.ga_destroy(ga_trw))
34     + call md_abort('Failed to destroy ga_trw',0)
35      endif
36      if(.not.ga_destroy(ga_trt))
37     + call md_abort('Failed to destroy ga_trt',0)
38      endif
39c
40      read(card(8:38),'(3i10,l1)') ifrfr,ifrto,ifrsk,lslvnt
41c
42      if(me.eq.0) then
43      write(*,'(/,3(a,i10))') ' Selected frames ',ifrfr,' to ',ifrto,
44     + ' by',ifrsk
45      endif
46      timoff=0.0d0
47      time=0.0d0
48      timr=0.0d0
49c
50c     allocate global array
51c
52c     nfrtot : total number of frames
53c     nfrdim : maximum number of frames per processor
54c     nfrme  : actual number of frames on processor
55c     nfrndx : global index of first element local list
56c     ntrj   : number of solute atoms per frame in memory resident frames
57c
58      nfrtot=int((ifrto-ifrfr+1)/ifrsk)
59      nfrdim=nfrtot/np
60      if(mod(nfrtot,np).ne.0) nfrdim=nfrdim+1
61      nfrme=min(nfrdim,nfrtot-me*nfrdim)
62      if(nfrme.lt.0) nfrme=0
63      nfrndx=me*nfrdim+1
64      ntrj=nsa
65c
66      if(me.eq.0) then
67      write(*,6001) np,nfrtot,nfrdim,nsel
68 6001 format(/,' Number of processors',t45,i5,/,
69     + ' Total number of selected frames',t45,i5,/,
70     + ' Number of frames per processor',t45,i5,//,
71     + ' Number of selected atoms',t45,i5,//)
72      endif
73c
74      if(.not.ga_create(mt_dbl,ntrj*3,nfrtot,'trj',ntrj*3,nfrdim,
75     + ga_trj)) call md_abort('Failed to create ga_trj',0)
76      if(lslvnt) then
77      if(.not.ga_create(mt_dbl,nwm*nwa*3,nfrtot,'trw',nwm*nwa*3,nfrdim,
78     + ga_trw)) call md_abort('Failed to create ga_trw',0)
79      endif
80      if(.not.ga_create(mt_dbl,6,nfrtot,'trt',6,nfrdim,
81     + ga_trt)) call md_abort('Failed to create ga_trt',0)
82c
83c      write(*,'(i5,a,2i5)') me,' allocated ga for frames ',
84c     + me*nfrdim+1,min((me+1)*nfrdim,nfrtot)
85c
86      if(iomode.eq.0) then
87c
88c     read trajectory files sequentially
89c     ----------------------------------
90c
91c     read file header
92c
93      call dia_rdhdr(byte_mb(i_snam))
94      nread=0
95      nwrit=0
96      nframe=0
97 1010 continue
98      if(dia_rdfram(dbl_mb(i_xdat),dbl_mb(i_wdat),t)) then
99      nread=nread+1
100      if(nread.lt.ifrfr) goto 1010
101      if(nframe.eq.(nframe/ifrsk)*ifrsk) then
102      ndata=ndata+1
103      if(me.eq.0) then
104      call ga_put(ga_trj,1,3*nsa,ndata,ndata,dbl_mb(i_xdat),nsa)
105      if(lslvnt) call ga_put(ga_trw,1,3*nwm*nwa,ndata,ndata,
106     + dbl_mb(i_wdat),nwm*nwa)
107      call ga_put(ga_trt,1,6,ndata,ndata,t,6)
108      endif
109      endif
110      nframe=nframe+1
111      if(ifrto.lt.ifrfr.or.nread.lt.ifrto) goto 1010
112      endif
113c
114      if(me.eq.0) then
115      close(unit=lfntrj)
116      write(*,'(a)') ' Closing trj file '
117      endif
118c
119      else
120c
121c     read trajectory files in parallel
122c     ---------------------------------
123c
124c     read file header
125c
126      call dia_rdhdr(byte_mb(i_snam))
127      call dia_getofffr()
128      iofffr=iofffr-ioffhd
129c
130      itmp(1)=nfilf
131      itmp(2)=ioffhd
132      itmp(3)=iofffr
133      itmp(4)=0
134      call ga_brdcst(mag_d01,itmp,4*ma_sizeof(mt_int,1,mt_byte),0)
135      nfilf=itmp(1)
136      ioffhd=itmp(2)
137      iofffr=itmp(3)
138c
139c      write(*,'(i5,a,i10)') me,' number of frames per file is ',nfilf
140c      write(*,'(i5,a,i10)') me,' trajectory header offset  is ',ioffhd
141c      write(*,'(i5,a,i10)') me,' trajectory frame size     is ',iofffr
142c
143c     meio   : processor carrying out i/o for current processor
144c     meframf: first frame to be read from trajectory for this processor
145c     meframl: last frame to be read from trajectory for this processor
146c     mefilf : file that has first frame for this processor
147c     mefill : file that has last frame for this processor
148c     meskip : number of frames to be skipped on first file for this processor
149c     melast : last frame to be read from the last file for this processor
150c
151      meio=(me/iomode)*iomode
152      meframf=ifrfr+me*nfrme*ifrsk
153      meframl=ifrfr+(me+1)*nfrme*ifrsk-1
154      mefilf=((meframf-1)/nfilf)+1
155      mefill=((meframl-1)/nfilf)+1
156      meskip=meframf-(mefilf-1)*nfilf-1
157      melast=meframl-(mefill-1)*nfilf
158c
159cc      write(*,'(i5,a,i5,a,2i5,a,2i5,a,2i5)')
160cc     + me,' p :',meio,':',meframf,meframl,':',mefilf,mefill,':',
161cc     + meskip,melast
162c
163c     since a single processor reads the trajectory for iomode processors
164c
165      if(me.eq.meio) then
166c
167c     meframl: last frame to be read from trajectory by the i/o processor
168c     mefill : file that has last frame for the i/o processor
169c     melast : last frame to be read from the last file for the i/o processor
170c
171      meframl=min(ifrfr+(me+iomode)*nfrme*ifrsk-1,ifrto)
172      mefill=((meframl-1)/nfilf)+1
173      melast=meframl-(mefill-1)*nfilf
174c
175c      write(*,'(i5,a,i5,a,2i5,a,2i5,a,2i5)')
176c     + me,' n :',meio,':',meframf,meframl,':',mefilf,mefill,':',
177c     + meskip,melast
178c
179c     loop ov er the number of frames to be stored by this processor
180c
181c
182      nfram=me*nfrme
183      lq=index(filtrj,'?')
184c
185c     loop over all files for this processor
186c
187      do 101 ifil=mefilf,mefill
188c
189c     construct file name for file ifil
190c
191      if(lq.gt.0) then
192      write(cnum,'(i3.3)') ifil
193      fname=filtrj(1:lq-1)//cnum//filtrj(lq+1:index(filtrj,' ')-1)
194      else
195      fname=filtrj(1:index(filtrj,' ')-1)
196      endif
197c
198c     determine first frame to read from this file
199c
200      iffram=1
201      ilfram=nfilf
202      if(ifil.eq.mefilf) iffram=meskip+1
203      if(ifil.eq.mefill) ilfram=melast
204c
205      open(unit=lfntrj,file=fname(1:index(fname,' ')-1),
206     + status='old',form='formatted')
207      do 102 ifr=iffram,ilfram
208      if(mod(ifr-iffram+1,ifrsk).eq.0) then
209      nfram=nfram+1
210      call fseek(lfntrj,ioffhd+(ifr-1)*iofffr,0)
211c      call fseek(lfntrj,ioffhd+(ifr-1)*iofffr,0,9999)
212      call dia_preadframe(dbl_mb(i_xdat),dbl_mb(i_wdat),t)
213      call ga_put(ga_trj,1,3*nsa,nfram,nfram,dbl_mb(i_xdat),nsa)
214      if(lslvnt) call ga_put(ga_trw,1,3*nwm*nwa,ndata,ndata,
215     + dbl_mb(i_wdat),nwm*nwa)
216      call ga_put(ga_trt,1,6,nfram,nfram,t,6)
217c      write(*,'(i5,6f12.6)') nfram,t
218      endif
219  102 continue
220  103 continue
221c
222      close(unit=lfntrj)
223c
224  101 continue
225      endif
226c
227      endif
228c
229      return
230 9999 continue
231      call md_abort('Error in fseek',0)
232      return
233      end
234      subroutine dia_puttrj()
235c
236c     read the trajectory
237c     -------------------
238c
239      implicit none
240c
241#include "dia_common.fh"
242#include "dia_params.fh"
243#include "global.fh"
244#include "mafdecls.fh"
245#include "rtdb.fh"
246#include "msgids.fh"
247#include "util.fh"
248c
249      return
250      end
251      subroutine dia_wrttrj(ga_t,fname,fmt,nbatch,logw)
252c
253c     read the trajectory
254c     -------------------
255c
256      implicit none
257c
258#include "dia_common.fh"
259#include "dia_params.fh"
260#include "global.fh"
261#include "mafdecls.fh"
262#include "rtdb.fh"
263#include "msgids.fh"
264#include "util.fh"
265c
266      integer ga_t,nbatch
267      character*255 fname
268      real*8 x(msa,3)
269      logical logw
270c
271      integer ifil,ibatch
272      character*3 fmt
273      integer i,ipd,ifr
274c
275      if(nbatch.le.0) return
276      if(me.ne.0) return
277c
278      ibatch=0
279      ifil=0
280      ipd=index(fname,'.')
281      if(ipd.eq.0) ipd=index(fname,' ')
282      lsonly=.true.
283c
284      do 1 ifr=1,nfrtot
285c
286      if(ibatch.eq.0) then
287      ifil=ifil+1
288      write(filcop,'(a,i5.5,a,a)') fname(1:ipd-1),ifil,'.',fmt
289      call dia_wthdr(fmt,byte_mb(i_snam),byte_mb(i_tag),int_mb(i_isel),
290     + logw)
291      endif
292      ibatch=ibatch+1
293c
294c     write frame
295c
296      call ga_get(ga_t,1,3*nsa,ifr,ifr,dbl_mb(i_xdat),nsa)
297      call dia_wtfram(fmt,byte_mb(i_snam),dbl_mb(i_xdat),dbl_mb(i_wdat),
298     + int_mb(i_isel),byte_mb(i_tag),dbl_mb(i_val),int_mb(i_wsel),
299     + int_mb(i_ndxw),ifr)
300c
301      if(ibatch.eq.nbatch.or.ifr.eq.nfrtot) then
302      ibatch=0
303      close(unit=lfncop)
304      endif
305c
306c
307    1 continue
308c
309      return
310      end
311      subroutine dia_setiomode()
312c
313      implicit none
314c
315#include "dia_common.fh"
316#include "dia_params.fh"
317#include "mafdecls.fh"
318#include "global.fh"
319#include "msgids.fh"
320c
321      read(card(8:12),1000) iomode
322 1000 format(i5)
323c
324      call ga_brdcst(mag_d01,iomode,ma_sizeof(mt_int,1,mt_byte),0)
325c
326      return
327      end
328      subroutine dia_copy()
329c
330      implicit none
331c
332#include "dia_common.fh"
333#include "dia_params.fh"
334#include "mafdecls.fh"
335#include "global.fh"
336#include "msgids.fh"
337c
338      character*255 fil,fname
339      character*3 fmt
340      integer i,nbatch
341c
342      read(card(8:14),1000) nbatch
343 1000 format(i7)
344      fil=card(15:80)
345      i=index(fil,'.')
346      if(i.eq.0) then
347      fmt='trj'
348      fname=fil(1:index(fil,' ')-1)//'.trj'
349      i=index(fname,'.')
350      else
351      fname=fil
352      fmt=fname(i+1:i+3)
353      endif
354c
355      call dia_wrttrj(ga_trj,fname(1:i-1),fmt,nbatch,.false.)
356c
357      return
358      end
359