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