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