source/wham/src-M/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.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       call xdrffloat_(ixdrf, rtime, iret)
  59 c      print *,"rtime",rtime," iret",iret
  60       call xdrffloat_(ixdrf, rpotE, iret)
  61 c      write (iout,*) "rpotE",rpotE," iret",iret
  62       call flush(iout)
  63       call xdrffloat_(ixdrf, ruconst, iret)
  64       call xdrffloat_(ixdrf, rt_bath, iret)
  65       call xdrfint_(ixdrf, nss, iret)
  66       do j=1,nss
  67            if (dyn_ss) then
  68             call xdrfint(ixdrf, idssb(j), iret)
  69             call xdrfint(ixdrf, jdssb(j), iret)
  70         idssb(j)=idssb(j)-nres
  71         jdssb(j)=jdssb(j)-nres
  72            else
  73             call xdrfint_(ixdrf, ihpb(j), iret)
  74             call xdrfint_(ixdrf, jhpb(j), iret)
  75            endif
  76       enddo
  77       call xdrfint_(ixdrf, nprop, iret)
  78       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
  79      &  call xdrfint(ixdrf, iset, iret)
  80       do i=1,nprop
  81         call xdrffloat_(ixdrf, rprop(i), iret)
  82       enddo
  83 #else
  84       call xdrffloat(ixdrf, rtime, iret)
  85       call xdrffloat(ixdrf, rpotE, iret)
  86 c      write (iout,*) "rpotE",rpotE," iret",iret
  87       call flush(iout)
  88       call xdrffloat(ixdrf, ruconst, iret)
  89       call xdrffloat(ixdrf, rt_bath, iret)
  90       call xdrfint(ixdrf, nss, iret)
  91       do j=1,nss
  92            if (dyn_ss) then
  93             call xdrfint(ixdrf, idssb(j), iret)
  94             call xdrfint(ixdrf, jdssb(j), iret)
  95            else
  96             call xdrfint(ixdrf, ihpb(j), iret)
  97             call xdrfint(ixdrf, jhpb(j), iret)
  98            endif
  99       enddo
 100       call xdrfint(ixdrf, nprop, iret)
 101 c      write (iout,*) "nprop",nprop
 102       if (it.gt.0 .and. nprop.ne.nprop_prev) then
 103         write (iout,*) "Warning previous nprop",nprop_prev,
 104      &   " current",nprop
 105         nprop=nprop_prev
 106       else
 107         nprop_prev=nprop
 108       endif
 109       call flush(iout)
 110       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
 111      &  call xdrfint(ixdrf, iset, iret)
 112       do i=1,nprop
 113         call xdrffloat(ixdrf, rprop(i), iret)
 114       enddo
 115 #endif
 116       if (iret.eq.0) exit
 117       itraj=mod(it,totraj(iR,iparm))
 118 #ifdef DEBUG
 119       write (iout,*) "ii",ii," itraj",itraj," it",it
 120 #endif
 121       if (iset.eq.0) iset = 1
 122       call flush(iout)
 123       it=it+1
 124       if (itraj.gt.ntraj) ntraj=itraj
 125       nstep(itraj)=nstep(itraj)+1
 126 c      rprop(2)=dsqrt(rprop(2))
 127 c      rprop(3)=dsqrt(rprop(3))
 128 #ifdef DEBUG
 129        write (iout,*) "umbrella ",umbrella
 130        write (iout,*) rtime,rpotE,rt_bath,nss,
 131      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
 132        write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
 133        call flush(iout)
 134 #endif
 135       prec=10000.0
 136
 137       itmp=0
 138 #if (defined(AIX) && !defined(JUBL))
 139       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 140 #else
 141       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 142 #endif
 143 #ifdef DEBUG
 144       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
 145 #endif
 146       if (iret.eq.0) exit
 147       if (itmp .ne. nres + nct - nnt + 1) then
 148         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 149         call flush(iout)
 150         exit
 151       endif
 152
 153       time=rtime
 154 c      write (iout,*) "calling slice"
 155 c      call flush(iout)
 156       islice=slice(nstep(itraj),time,is,ie,ts,te)
 157 c      write (iout,*) "islice",islice
 158 c      call flush(iout)
 159
 160       do i=1,nres
 161         do j=1,3
 162           c(j,i)=xoord(j,i)
 163         enddo
 164       enddo
 165       do i=1,nct-nnt+1
 166         do j=1,3
 167           c(j,i+nres+nnt-1)=xoord(j,i+nres)
 168         enddo
 169       enddo
 170
 171 c Box shift
 172       call oligomer
 173       do i=1,nres
 174         do j=1,3
 175           xoord(j,i)=c(j,i)
 176         enddo
 177       enddo
 178       do i=1,nct-nnt+1
 179         do j=1,3
 180           xoord(j,i+nres)=c(j,i+nres+nnt-1)
 181         enddo
 182       enddo
 183 c end change
 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