source/wham/src-M/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            if (dyn_ss) then
  67             call xdrfint(ixdrf, idssb(j)+nres, iret)
  68             call xdrfint(ixdrf, jdssb(j)+nres, iret)
  69            else
  70             call xdrfint_(ixdrf, ihpb(j), iret)
  71             call xdrfint_(ixdrf, jhpb(j), iret)
  72            endif
  73       enddo
  74       call xdrfint_(ixdrf, nprop, iret)
  75       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
  76      &  call xdrfint(ixdrf, iset, iret)
  77       do i=1,nprop
  78         call xdrffloat_(ixdrf, rprop(i), iret)
  79       enddo
  80 #else
  81       call xdrffloat(ixdrf, rtime, iret)
  82       call xdrffloat(ixdrf, rpotE, iret)
  83 c      write (iout,*) "rpotE",rpotE," iret",iret
  84       call flush(iout)
  85       call xdrffloat(ixdrf, ruconst, iret)
  86       call xdrffloat(ixdrf, rt_bath, iret)
  87       call xdrfint(ixdrf, nss, iret)
  88       do j=1,nss
  89            if (dyn_ss) then
  90             call xdrfint(ixdrf, idssb(j)+nres, iret)
  91             call xdrfint(ixdrf, jdssb(j)+nres, iret)
  92            else
  93             call xdrfint(ixdrf, ihpb(j), iret)
  94             call xdrfint(ixdrf, jhpb(j), iret)
  95            endif
  96       enddo
  97       call xdrfint(ixdrf, nprop, iret)
  98 c      write (iout,*) "nprop",nprop
  99       if (it.gt.0 .and. nprop.ne.nprop_prev) then
 100         write (iout,*) "Warning previous nprop",nprop_prev,
 101      &   " current",nprop
 102         nprop=nprop_prev
 103       else
 104         nprop_prev=nprop
 105       endif
 106       call flush(iout)
 107       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
 108      &  call xdrfint(ixdrf, iset, iret)
 109       do i=1,nprop
 110         call xdrffloat(ixdrf, rprop(i), iret)
 111       enddo
 112 #endif
 113       if (iret.eq.0) exit
 114       itraj=mod(it,totraj(iR,iparm))
 115 #ifdef DEBUG
 116       write (iout,*) "ii",ii," itraj",itraj," it",it
 117 #endif
 118       if (iset.eq.0) iset = 1
 119       call flush(iout)
 120       it=it+1
 121       if (itraj.gt.ntraj) ntraj=itraj
 122       nstep(itraj)=nstep(itraj)+1
 123 c      rprop(2)=dsqrt(rprop(2))
 124 c      rprop(3)=dsqrt(rprop(3))
 125 #ifdef DEBUG
 126        write (iout,*) "umbrella ",umbrella
 127        write (iout,*) rtime,rpotE,rt_bath,nss,
 128      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
 129        write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
 130        call flush(iout)
 131 #endif
 132       prec=10000.0
 133
 134       itmp=0
 135 #if (defined(AIX) && !defined(JUBL))
 136       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 137 #else
 138       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 139 #endif
 140 #ifdef DEBUG
 141       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
 142 #endif
 143       if (iret.eq.0) exit
 144       if (itmp .ne. nres + nct - nnt + 1) then
 145         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 146         call flush(iout)
 147         exit
 148       endif
 149
 150       time=rtime
 151 c      write (iout,*) "calling slice"
 152 c      call flush(iout)
 153       islice=slice(nstep(itraj),time,is,ie,ts,te)
 154 c      write (iout,*) "islice",islice
 155 c      call flush(iout)
 156
 157       do i=1,nres
 158         do j=1,3
 159           c(j,i)=xoord(j,i)
 160         enddo
 161       enddo
 162       do i=1,nct-nnt+1
 163         do j=1,3
 164           c(j,i+nres+nnt-1)=xoord(j,i+nres)
 165         enddo
 166       enddo
 167
 168       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset
 169      &    .or. iset.eq.myparm)) then
 170         ii=ii+1
 171         kk(islice)=kk(islice)+1
 172         mm(islice)=mm(islice)+1
 173         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and.
 174      &     conf_check(ll(islice)+1,1)) then
 175           if (replica(iparm)) then
 176              rt_bath=1.0d0/(rt_bath*1.987D-3)
 177              do i=1,nT_h(iparm)
 178                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 179                  iib = i
 180                  goto 22
 181                endif
 182              enddo
 183   22         continue
 184              if (i.gt.nT_h(iparm)) then
 185                write (iout,*) "Error - temperature of conformation",
 186      &         ii,1.0d0/(rt_bath*1.987D-3),
 187      &         " does not match any of the list"
 188                write (iout,*)
 189      &          1.0d0/(rt_bath*1.987D-3),
 190      &          (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 191                call flush(iout)
 192 c               exit
 193 c               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 194                ii=ii-1
 195                kk(islice)=kk(islice)-1
 196                mm(islice)=mm(islice)-1
 197                goto 112
 198              endif
 199           else
 200             iib = ib
 201           endif
 202
 203           efree=0.0d0
 204           jj(islice)=jj(islice)+1
 205           if (umbrella(iparm)) then
 206             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 207           else if (hamil_rep) then
 208             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
 209           else
 210             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 211           endif
 212           ll(islice)=ll(islice)+1
 213 #ifdef DEBUG
 214           write (iout,*) "Writing conformation, record",ll(islice)
 215           write (iout,*) "ib",ib," iib",iib
 216           write (iout,*) "ntraj",ntraj," itraj",itraj,
 217      &      " nstep",nstep(itraj)
 218           write (iout,*) "pote",rpotE," time",rtime
 219 c          if (replica(iparm)) then
 220 c            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 221 c            write (iout,*) "TEMP list"
 222 c            write (iout,*)
 223 c     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 224 c          endif
 225           write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 226 c          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 227 c          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 228           call flush(iout)
 229 #endif
 230           if (islice.ne.islice1) then
 231 c            write (iout,*) "islice",islice," islice1",islice1
 232             close(ientout)
 233 c            write (iout,*) "Closing file ",
 234 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 235             call opentmp(islice,ientout,bprotfile_temp)
 236 c            write (iout,*) "Opening file ",
 237 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 238             islice1=islice
 239           endif
 240           if (umbrella(iparm)) then
 241             write(ientout,rec=ll(islice))
 242      &        ((xoord(l,k),l=1,3),k=1,nres),
 243      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 244      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 245      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 246      &        iset,iib,iparm
 247           else if (hamil_rep) then
 248             write(ientout,rec=ll(islice))
 249      &        ((xoord(l,k),l=1,3),k=1,nres),
 250      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 251      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 252      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 253      &        iR,iib,iset
 254           else
 255             write(ientout,rec=ll(islice))
 256      &        ((xoord(l,k),l=1,3),k=1,nres),
 257      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 258      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 259      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 260      &        iR,iib,iparm
 261           endif
 262 #ifdef DEBUG
 263           call int_from_cart1(.false.)
 264           write (iout,*) "Writing conformation, record",ll(islice)
 265           write (iout,*) "Cartesian coordinates"
 266           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 267           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 268           write (iout,*) "Internal coordinates"
 269           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 270           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 271           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 272           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 273           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 274           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 275           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 276 c          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 277           write (iout,'(16i5)') iscor
 278           call flush(iout)
 279 #endif
 280         endif
 281       endif
 282
 283   112 continue
 284
 285       enddo
 286       close(ientout)
 287 #if (defined(AIX) && !defined(JUBL))
 288       call xdrfclose_(ixdrf, iret)
 289 #else
 290       call xdrfclose(ixdrf, iret)
 291 #endif
 292       write (iout,'(i10," trajectories found in file.")') ntraj+1
 293       write (iout,'(a)') "Numbers of steps in trajectories:"
 294       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 295       write (iout,*) ii," conformations read from file",
 296      &   nazwa(:ilen(nazwa))
 297       do islice=1,nslice
 298         write (iout,*) mm(islice)," conformations read so far, slice",
 299      &    islice
 300         write (iout,*) ll(islice),
 301      &  " conformations stored so far, slice",islice
 302       enddo
 303       call flush(iout)
 304       return
 305       end