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