source/wham/src-M-SAXS.safe/cxread.F.safe

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