source/wham/src-new/cxread.F.orig

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