source/wham/src/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 #ifdef DEBUG
  58       write (iout,*) "ii",ii," itraj",itraj," it",it
  59 #endif
  60       call xdrffloat_(ixdrf, rtime, iret)
  61       call xdrffloat_(ixdrf, rpotE, iret)
  62 #ifdef DEBUG
  63       write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
  64 #endif
  65       call flush(iout)
  66       call xdrffloat_(ixdrf, ruconst, iret)
  67       call xdrffloat_(ixdrf, rt_bath, iret)
  68       call xdrfint_(ixdrf, nss, iret)
  69 #ifdef DEBUG
  70       write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
  71 #endif
  72       do j=1,nss
  73         call xdrfint_(ixdrf, ihpb(j), iret)
  74         call xdrfint_(ixdrf, jhpb(j), iret)
  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 #ifdef DEBUG
  84       write (iout,*) "ii",ii," itraj",itraj," it",it
  85 #endif
  86       call xdrffloat(ixdrf, rtime, iret)
  87       call xdrffloat(ixdrf, rpotE, iret)
  88 #ifdef DEBUG
  89       write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
  90 #endif
  91       call flush(iout)
  92       call xdrffloat(ixdrf, ruconst, iret)
  93       call xdrffloat(ixdrf, rt_bath, iret)
  94       call xdrfint(ixdrf, nss, iret)
  95 #ifdef DEBUG
  96       write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
  97 #endif
  98       do j=1,nss
  99         call xdrfint(ixdrf, ihpb(j), iret)
 100         call xdrfint(ixdrf, jhpb(j), iret)
 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       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