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