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*256 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      write (iout,*) "Before boxshift"
 175 c      call flush(iout)
 176 c Box shift
 177       call oligomer
 178 c      write (iout,*) "After oligomer"
 179 c      call flush(iout)
 180       do i=1,nres
 181         do j=1,3
 182           xoord(j,i)=c(j,i)
 183         enddo
 184       enddo
 185       do i=1,nct-nnt+1
 186         do j=1,3
 187           xoord(j,i+nres)=c(j,i+nres+nnt-1)
 188         enddo
 189       enddo
 190 c end change
 191 c      write (iout,*) "Before islice"
 192 c      call flush(iout)
 193       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset
 194      &    .or. iset.eq.myparm)) then
 195         ii=ii+1
 196         kk(islice)=kk(islice)+1
 197         mm(islice)=mm(islice)+1
 198 #ifdef DEBUG
 199         write (iout,*) "islice",islice," ii",ii," kk",kk(islice),
 200      &    " mm",mm(islice)
 201         write (iout,*) "itraj",itraj," nstep",nstep(itraj),
 202      &    " isampl",isampl(iparm)
 203         call flush(iout)
 204 #endif
 205         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and.
 206      &     conf_check(ll(islice)+1,1)) then
 207           if (replica(iparm)) then
 208              if (rt_bath.eq.0.0d0) then
 209                write (iout,*) "ERROR: zero temperature",
 210      &         islice,kk(islice),mm(islice)
 211                call flush(iout)
 212              endif
 213              rt_bath=1.0d0/(rt_bath*1.987D-3)
 214              do i=1,nT_h(iparm)
 215                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 216                  iib = i
 217                  goto 22
 218                endif
 219              enddo
 220   22         continue
 221              if (i.gt.nT_h(iparm)) then
 222                write (iout,*) "Error - temperature of conformation",
 223      &         ii,1.0d0/(rt_bath*1.987D-3),
 224      &         " does not match any of the list"
 225                write (iout,*)
 226      &          1.0d0/(rt_bath*1.987D-3),
 227      &          (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 228                call flush(iout)
 229 c               exit
 230 c               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 231                ii=ii-1
 232                kk(islice)=kk(islice)-1
 233                mm(islice)=mm(islice)-1
 234                goto 112
 235              endif
 236           else
 237             iib = ib
 238           endif
 239
 240           efree=0.0d0
 241           jj(islice)=jj(islice)+1
 242           if (umbrella(iparm)) then
 243             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 244           else if (hamil_rep) then
 245             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
 246           else
 247             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 248           endif
 249           ll(islice)=ll(islice)+1
 250 #ifdef DEBUG
 251           write (iout,*) "Writing conformation, record",ll(islice)
 252           write (iout,*) "ib",ib," iib",iib
 253           write (iout,*) "ntraj",ntraj," itraj",itraj,
 254      &      " nstep",nstep(itraj)
 255           write (iout,*) "pote",rpotE," time",rtime
 256           write (iout,*) "nss",nss
 257           write (iout,*) (ihpb(k),jhpb(k),k=1,nss)
 258 c          if (replica(iparm)) then
 259 c            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 260 c            write (iout,*) "TEMP list"
 261 c            write (iout,*)
 262 c     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 263 c          endif
 264 c          write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 265 c          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 266 c          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 267           call flush(iout)
 268 #endif
 269           if (islice.ne.islice1) then
 270 c            write (iout,*) "islice",islice," islice1",islice1
 271             close(ientout)
 272 c            write (iout,*) "Closing file ",
 273 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 274             call opentmp(islice,ientout,bprotfile_temp)
 275 c            write (iout,*) "Opening file ",
 276 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 277             islice1=islice
 278           endif
 279           if (umbrella(iparm)) then
 280             write(ientout,rec=ll(islice))
 281      &        ((xoord(l,k),l=1,3),k=1,nres),
 282      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 283      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 284      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 285      &        iset,iib,iparm
 286           else if (hamil_rep) then
 287             write(ientout,rec=ll(islice))
 288      &        ((xoord(l,k),l=1,3),k=1,nres),
 289      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 290      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 291      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 292      &        iR,iib,iset
 293           else
 294             write(ientout,rec=ll(islice))
 295      &        ((xoord(l,k),l=1,3),k=1,nres),
 296      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 297      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 298      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 299      &        iR,iib,iparm
 300           endif
 301 #ifdef DEBUG
 302           call int_from_cart1(.false.)
 303           write (iout,*) "Writing conformation, record",ll(islice)
 304           write (iout,*) "Cartesian coordinates"
 305           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 306           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 307           write (iout,*) "Internal coordinates"
 308           call intout
 309           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 310           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 311           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 312           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 313           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 314           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 315           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 316 c          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 317           write (iout,'(16i5)') iscor
 318           call flush(iout)
 319 #endif
 320         endif
 321       endif
 322
 323   112 continue
 324
 325       enddo
 326       close(ientout)
 327 #if (defined(AIX) && !defined(JUBL))
 328       call xdrfclose_(ixdrf, iret)
 329 #else
 330       call xdrfclose(ixdrf, iret)
 331 #endif
 332       write (iout,'(i10," trajectories found in file.")') ntraj+1
 333       write (iout,'(a)') "Numbers of steps in trajectories:"
 334       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 335       write (iout,*) ii," conformations read from file",
 336      &   nazwa(:ilen(nazwa))
 337       do islice=1,nslice
 338         write (iout,*) mm(islice)," conformations read so far, slice",
 339      &    islice
 340         write (iout,*) ll(islice),
 341      &  " conformations stored so far, slice",islice
 342       enddo
 343       call flush(iout)
 344       return
 345       end