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