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