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