source/wham/src-M-Pawel/cxread.F

   1       subroutine cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,*)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'DIMENSIONS.FREE'
   6       integer MaxTraj
   7       parameter (MaxTraj=2050)
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.NAMES'
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.HEADER'
  13       include 'COMMON.SBRIDGE'
  14       include 'COMMON.PROTFILES'
  15       include 'COMMON.OBCINKA'
  16       include 'COMMON.FREE'
  17       include 'COMMON.VAR'
  18       include 'COMMON.GEO'
  19       include 'COMMON.PROT'
  20       character*64 nazwa,bprotfile_temp
  21       real*4 rtime,rpotE,ruconst,rt_bath,rprop(maxQ)
  22       double precision time
  23       integer iret,itmp,itraj,ntraj
  24       real xoord(3,maxres2+2),prec
  25       integer nstep(0:MaxTraj-1)
  26       integer ilen
  27       external ilen
  28       integer ii,jj(maxslice),kk(maxslice),ll(maxslice),mm(maxslice)
  29       integer is(MaxSlice),ie(MaxSlice),nrec_slice
  30       double precision ts(MaxSlice),te(MaxSlice),time_slice
  31       integer slice
  32       logical conf_check
  33       call set_slices(is,ie,ts,te,iR,ib,iparm)
  34
  35       do i=1,nQ
  36         rprop(i)=0.0d0
  37       enddo
  38       do i=0,MaxTraj-1
  39         nstep(i)=0
  40       enddo
  41       ntraj=0
  42       it=0
  43       iret=1
  44 #if (defined(AIX) && !defined(JUBL))
  45       call xdrfopen_(ixdrf,nazwa, "r", iret)
  46 #else
  47       call xdrfopen(ixdrf,nazwa, "r", iret)
  48 #endif
  49       if (iret.eq.0) return1
  50
  51       islice1=1
  52       call opentmp(islice1,ientout,bprotfile_temp)
  53 c      print *,"bumbum"
  54       do while (iret.gt.0)
  55
  56 #if (defined(AIX) && !defined(JUBL))
  57       call xdrffloat_(ixdrf, rtime, iret)
  58 c      print *,"rtime",rtime," iret",iret
  59       call xdrffloat_(ixdrf, rpotE, iret)
  60 c      write (iout,*) "rpotE",rpotE," iret",iret
  61       call flush(iout)
  62       call xdrffloat_(ixdrf, ruconst, iret)
  63       call xdrffloat_(ixdrf, rt_bath, iret)
  64       call xdrfint_(ixdrf, nss, iret)
  65       do j=1,nss
  66         call xdrfint_(ixdrf, ihpb(j), iret)
  67         call xdrfint_(ixdrf, jhpb(j), iret)
  68       enddo
  69       call xdrfint_(ixdrf, nprop, iret)
  70       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
  71      &  call xdrfint(ixdrf, iset, iret)
  72       do i=1,nprop
  73         call xdrffloat_(ixdrf, rprop(i), iret)
  74       enddo
  75 #else
  76       call xdrffloat(ixdrf, rtime, iret)
  77       call xdrffloat(ixdrf, rpotE, iret)
  78 c      write (iout,*) "rpotE",rpotE," iret",iret
  79       call flush(iout)
  80       call xdrffloat(ixdrf, ruconst, iret)
  81       call xdrffloat(ixdrf, rt_bath, iret)
  82       call xdrfint(ixdrf, nss, iret)
  83       do j=1,nss
  84         call xdrfint(ixdrf, ihpb(j), iret)
  85         call xdrfint(ixdrf, jhpb(j), iret)
  86       enddo
  87       call xdrfint(ixdrf, nprop, iret)
  88 c      write (iout,*) "nprop",nprop
  89       if (it.gt.0 .and. nprop.ne.nprop_prev) then
  90         write (iout,*) "Warning previous nprop",nprop_prev,
  91      &   " current",nprop
  92         nprop=nprop_prev
  93       else
  94         nprop_prev=nprop
  95       endif
  96       call flush(iout)
  97       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
  98      &  call xdrfint(ixdrf, iset, iret)
  99       do i=1,nprop
 100         call xdrffloat(ixdrf, rprop(i), iret)
 101       enddo
 102 #endif
 103       if (iret.eq.0) exit
 104       itraj=mod(it,totraj(iR,iparm))
 105 #ifdef DEBUG
 106       write (iout,*) "ii",ii," itraj",itraj," it",it
 107 #endif
 108       if (iset.eq.0) iset = 1
 109       call flush(iout)
 110       it=it+1
 111       if (itraj.gt.ntraj) ntraj=itraj
 112       nstep(itraj)=nstep(itraj)+1
 113 c      rprop(2)=dsqrt(rprop(2))
 114 c      rprop(3)=dsqrt(rprop(3))
 115 #ifdef DEBUG
 116        write (iout,*) "umbrella ",umbrella
 117        write (iout,*) rtime,rpotE,rt_bath,nss,
 118      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
 119        write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
 120        call flush(iout)
 121 #endif
 122       prec=10000.0
 123
 124       itmp=0
 125 #if (defined(AIX) && !defined(JUBL))
 126       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 127 #else
 128       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 129 #endif
 130 #ifdef DEBUG
 131       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
 132 #endif
 133       if (iret.eq.0) exit
 134       if (itmp .ne. nres + nct - nnt + 1) then
 135         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 136         call flush(iout)
 137         exit
 138       endif
 139
 140       time=rtime
 141 c      write (iout,*) "calling slice"
 142 c      call flush(iout)
 143       islice=slice(nstep(itraj),time,is,ie,ts,te)
 144 c      write (iout,*) "islice",islice
 145 c      call flush(iout)
 146
 147       do i=1,nres
 148         do j=1,3
 149           c(j,i)=xoord(j,i)
 150         enddo
 151       enddo
 152       do i=1,nct-nnt+1
 153         do j=1,3
 154           c(j,i+nres+nnt-1)=xoord(j,i+nres)
 155         enddo
 156       enddo
 157
 158       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset
 159      &    .or. iset.eq.myparm)) then
 160         ii=ii+1
 161         kk(islice)=kk(islice)+1
 162         mm(islice)=mm(islice)+1
 163         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and.
 164      &     conf_check(ll(islice)+1,1)) then
 165           if (replica(iparm)) then
 166              rt_bath=1.0d0/(rt_bath*1.987D-3)
 167              do i=1,nT_h(iparm)
 168                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 169                  iib = i
 170                  goto 22
 171                endif
 172              enddo
 173   22         continue
 174              if (i.gt.nT_h(iparm)) then
 175                write (iout,*) "Error - temperature of conformation",
 176      &         ii,1.0d0/(rt_bath*1.987D-3),
 177      &         " does not match any of the list"
 178                write (iout,*)
 179      &          1.0d0/(rt_bath*1.987D-3),
 180      &          (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 181                call flush(iout)
 182 c               exit
 183 c               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 184                ii=ii-1
 185                kk(islice)=kk(islice)-1
 186                mm(islice)=mm(islice)-1
 187                goto 112
 188              endif
 189           else
 190             iib = ib
 191           endif
 192
 193           efree=0.0d0
 194           jj(islice)=jj(islice)+1
 195           if (umbrella(iparm)) then
 196             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 197           else if (hamil_rep) then
 198             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
 199           else
 200             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 201           endif
 202           ll(islice)=ll(islice)+1
 203 #ifdef DEBUG
 204           write (iout,*) "Writing conformation, record",ll(islice)
 205           write (iout,*) "ib",ib," iib",iib
 206           write (iout,*) "ntraj",ntraj," itraj",itraj,
 207      &      " nstep",nstep(itraj)
 208           write (iout,*) "pote",rpotE," time",rtime
 209 c          if (replica(iparm)) then
 210 c            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 211 c            write (iout,*) "TEMP list"
 212 c            write (iout,*)
 213 c     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 214 c          endif
 215           write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 216 c          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 217 c          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 218           call flush(iout)
 219 #endif
 220           if (islice.ne.islice1) then
 221 c            write (iout,*) "islice",islice," islice1",islice1
 222             close(ientout)
 223 c            write (iout,*) "Closing file ",
 224 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 225             call opentmp(islice,ientout,bprotfile_temp)
 226 c            write (iout,*) "Opening file ",
 227 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 228             islice1=islice
 229           endif
 230           if (umbrella(iparm)) then
 231             write(ientout,rec=ll(islice))
 232      &        ((xoord(l,k),l=1,3),k=1,nres),
 233      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 234      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 235      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 236      &        iset,iib,iparm
 237           else if (hamil_rep) then
 238             write(ientout,rec=ll(islice))
 239      &        ((xoord(l,k),l=1,3),k=1,nres),
 240      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 241      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 242      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 243      &        iR,iib,iset
 244           else
 245             write(ientout,rec=ll(islice))
 246      &        ((xoord(l,k),l=1,3),k=1,nres),
 247      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 248      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 249      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 250      &        iR,iib,iparm
 251           endif
 252 #ifdef DEBUG
 253           call int_from_cart1(.false.)
 254           write (iout,*) "Writing conformation, record",ll(islice)
 255           write (iout,*) "Cartesian coordinates"
 256           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 257           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 258           write (iout,*) "Internal coordinates"
 259           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 260           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 261           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 262           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 263           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 264           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 265           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 266 c          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 267           write (iout,'(16i5)') iscor
 268           call flush(iout)
 269 #endif
 270         endif
 271       endif
 272
 273   112 continue
 274
 275       enddo
 276       close(ientout)
 277 #if (defined(AIX) && !defined(JUBL))
 278       call xdrfclose_(ixdrf, iret)
 279 #else
 280       call xdrfclose(ixdrf, iret)
 281 #endif
 282       write (iout,'(i10," trajectories found in file.")') ntraj+1
 283       write (iout,'(a)') "Numbers of steps in trajectories:"
 284       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 285       write (iout,*) ii," conformations read from file",
 286      &   nazwa(:ilen(nazwa))
 287       do islice=1,nslice
 288         write (iout,*) mm(islice)," conformations read so far, slice",
 289      &    islice
 290         write (iout,*) ll(islice),
 291      &  " conformations stored so far, slice",islice
 292       enddo
 293       call flush(iout)
 294       return
 295       end