source/wham/src/cxread.F.org

   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       call set_slices(is,ie,ts,te,iR,ib,iparm)
  33
  34       do i=1,nQ
  35         rprop(i)=0.0d0
  36       enddo
  37       do i=0,MaxTraj-1
  38         nstep(i)=0
  39       enddo
  40       ntraj=0
  41       it=0
  42       iret=1
  43 #if (defined(AIX) && !defined(JUBL))
  44       call xdrfopen_(ixdrf,nazwa, "r", iret)
  45 #else
  46       call xdrfopen(ixdrf,nazwa, "r", iret)
  47 #endif
  48       if (iret.eq.0) return1
  49
  50       islice1=1
  51       call opentmp(islice1,ientout,bprotfile_temp)
  52 c      print *,"bumbum"
  53       do while (iret.gt.0)
  54
  55 #if (defined(AIX) && !defined(JUBL))
  56       call xdrffloat_(ixdrf, rtime, iret)
  57 c      print *,"rtime",rtime," iret",iret
  58       call xdrffloat_(ixdrf, rpotE, iret)
  59 c      write (iout,*) "rpotE",rpotE," iret",iret
  60       call flush(iout)
  61       call xdrffloat_(ixdrf, ruconst, iret)
  62       call xdrffloat_(ixdrf, rt_bath, iret)
  63       call xdrfint_(ixdrf, nss, iret)
  64       do j=1,nss
  65         call xdrfint_(ixdrf, ihpb(j), iret)
  66         call xdrfint_(ixdrf, jhpb(j), iret)
  67       enddo
  68       call xdrfint_(ixdrf, nprop, iret)
  69       do i=1,nprop
  70         call xdrffloat_(ixdrf, rprop(i), iret)
  71       enddo
  72 #else
  73       call xdrffloat(ixdrf, rtime, iret)
  74       call xdrffloat(ixdrf, rpotE, iret)
  75 c      write (iout,*) "rpotE",rpotE," iret",iret
  76       call flush(iout)
  77       call xdrffloat(ixdrf, ruconst, iret)
  78       call xdrffloat(ixdrf, rt_bath, iret)
  79       call xdrfint(ixdrf, nss, iret)
  80       do j=1,nss
  81         call xdrfint(ixdrf, ihpb(j), iret)
  82         call xdrfint(ixdrf, jhpb(j), iret)
  83       enddo
  84       call xdrfint(ixdrf, nprop, iret)
  85 c      write (iout,*) "nprop",nprop
  86       call flush(iout)
  87       do i=1,nprop
  88         call xdrffloat(ixdrf, rprop(i), iret)
  89       enddo
  90 #endif
  91       if (iret.eq.0) exit
  92       itraj=mod(it,totraj(iR,iparm))
  93 #ifdef DEBUG
  94       write (iout,*) "ii",ii," itraj",itraj
  95 #endif
  96       call flush(iout)
  97       it=it+1
  98       if (itraj.gt.ntraj) ntraj=itraj
  99       nstep(itraj)=nstep(itraj)+1
 100 #ifdef DEBUG
 101        write (iout,*) rtime,rpotE,rt_bath,nss,
 102      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
 103        call flush(iout)
 104 #endif
 105       prec=10000.0
 106
 107       itmp=0
 108 #if (defined(AIX) && !defined(JUBL))
 109       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 110 #else
 111       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 112 #endif
 113 #ifdef DEBUG
 114       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
 115 #endif
 116       if (iret.eq.0) exit
 117       if (itmp .ne. nres + nct - nnt + 1) then
 118         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 119         call flush(iout)
 120         exit
 121       endif
 122
 123       time=rtime
 124 c      write (iout,*) "calling slice"
 125 c      call flush(iout)
 126       islice=slice(nstep(itraj),time,is,ie,ts,te)
 127 c      write (iout,*) "islice",islice
 128 c      call flush(iout)
 129
 130       if (islice.gt.0 .and. islice.le.nslice) then
 131         ii=ii+1
 132         kk(islice)=kk(islice)+1
 133         mm(islice)=mm(islice)+1
 134         if (mod(nstep(itraj),isampl(iparm)).eq.0) then
 135             if (replica(iparm)) then
 136                rt_bath=1.0d0/(rt_bath*1.987D-3)
 137                do i=1,nT_h(iparm)
 138                  if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 139                    iib = i
 140                    goto 22
 141                  endif
 142                enddo
 143   22           continue
 144                if (i.gt.nT_h(iparm)) then
 145                  write (iout,*) "Error - temperature of conformation",
 146      &           ii,1.0d0/(rt_bath*1.987D-3),
 147      &           " does not match any of the list"
 148                  write (iout,*)
 149      &            1.0d0/(rt_bath*1.987D-3),
 150      &            (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 151                  call flush(iout)
 152                  exit
 153                  call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 154                endif
 155             else
 156                 iib = ib
 157             endif
 158
 159             efree=0.0d0
 160             jj(islice)=jj(islice)+1
 161             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 162             ll(islice)=ll(islice)+1
 163 #ifdef DEBUG
 164             write (iout,*) "Writing conformation, record",ll(islice)
 165             write (iout,*) "ib",ib," iib",iib
 166             write (iout,*) "ntraj",ntraj," itraj",itraj,
 167      &        " nstep",nstep(itraj)
 168             write (iout,*) "pote",rpotE," time",rtime
 169 c            if (replica(iparm)) then
 170 c              write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 171 c              write (iout,*) "TEMP list"
 172 c              write (iout,*)
 173 c     &         (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 174 c            endif
 175             write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 176 c            write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 177 c            write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 178             call flush(iout)
 179 #endif
 180             if (islice.ne.islice1) then
 181 c              write (iout,*) "islice",islice," islice1",islice1
 182               close(ientout)
 183 c              write (iout,*) "Closing file ",
 184 c     &            bprotfile_temp(:ilen(bprotfile_temp))
 185               call opentmp(islice,ientout,bprotfile_temp)
 186 c              write (iout,*) "Opening file ",
 187 c     &            bprotfile_temp(:ilen(bprotfile_temp))
 188               islice1=islice
 189             endif
 190             write(ientout,rec=ll(islice))
 191      &        ((xoord(l,k),l=1,3),k=1,nres),
 192      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 193      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 194      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 195      &        iR,iib,iparm
 196 #ifdef DEBUG
 197             do i=1,nres
 198               do j=1,3
 199                 c(j,i)=xoord(j,i)
 200               enddo
 201             enddo
 202             do i=1,nct-nnt+1
 203               do j=1,3
 204                 c(j,i+nres+nnt-1)=xoord(j,i+nres)
 205               enddo
 206             enddo
 207             call int_from_cart1(.false.)
 208             write (iout,*) "Writing conformation, record",ll(islice)
 209             write (iout,*) "Cartesian coordinates"
 210             write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 211             write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 212             write (iout,*) "Internal coordinates"
 213             write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 214             write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 215             write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 216             write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 217             write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 218             write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 219             write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 220 c            write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 221             write (iout,'(16i5)') iscor
 222             call flush(iout)
 223 #endif
 224         endif
 225       endif
 226
 227       enddo
 228   112 continue
 229       close(ientout)
 230 #if (defined(AIX) && !defined(JUBL))
 231       call xdrfclose_(ixdrf, iret)
 232 #else
 233       call xdrfclose(ixdrf, iret)
 234 #endif
 235       write (iout,'(i10," trajectories found in file.")') ntraj+1
 236       write (iout,'(a)') "Numbers of steps in trajectories:"
 237       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 238       write (iout,*) ii," conformations read from file",
 239      &   nazwa(:ilen(nazwa))
 240       do islice=1,nslice
 241         write (iout,*) mm(islice)," conformations read so far, slice",
 242      &    islice
 243         write (iout,*) ll(islice),
 244      &  " conformations stored so far, slice",islice
 245       enddo
 246       call flush(iout)
 247       return
 248       end