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