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