source/wham/src-NEWSC-NEWCORR/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 #ifdef MPI
   7       include "mpif.h"
   8       include "COMMON.MPI"
   9 #endif
  10       integer MaxTraj
  11       parameter (MaxTraj=2050)
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.INTERACT'
  14       include 'COMMON.NAMES'
  15       include 'COMMON.IOUNITS'
  16       include 'COMMON.HEADER'
  17       include 'COMMON.SBRIDGE'
  18       include 'COMMON.PROTFILES'
  19       include 'COMMON.OBCINKA'
  20       include 'COMMON.FREE'
  21       include 'COMMON.VAR'
  22       include 'COMMON.GEO'
  23       include 'COMMON.PROT'
  24       character*64 nazwa,bprotfile_temp
  25       real*4 rtime,rpotE,ruconst,rt_bath,rprop(maxQ)
  26       double precision time
  27       integer iret,itmp,itraj,ntraj
  28       real xoord(3,maxres2+2),prec
  29       integer nstep(0:MaxTraj-1)
  30       integer ilen
  31       external ilen
  32       integer ii,jj(maxslice),kk(maxslice),ll(maxslice),mm(maxslice)
  33       integer is(MaxSlice),ie(MaxSlice),nrec_slice
  34       double precision ts(MaxSlice),te(MaxSlice),time_slice
  35       integer slice
  36       logical conf_check
  37       character*4 lt_bath
  38       character*256 pdbfilename
  39       character*50 tytul
  40       call set_slices(is,ie,ts,te,iR,ib,iparm)
  41
  42       do i=1,nQ
  43         rprop(i)=0.0d0
  44       enddo
  45       do i=0,MaxTraj-1
  46         nstep(i)=0
  47       enddo
  48       ntraj=0
  49       it=0
  50       iret=1
  51 #if (defined(AIX) && !defined(JUBL))
  52       call xdrfopen_(ixdrf,nazwa, "r", iret)
  53 #else
  54       call xdrfopen(ixdrf,nazwa, "r", iret)
  55 #endif
  56       if (iret.eq.0) return1
  57
  58       islice1=1
  59       call opentmp(islice1,ientout,bprotfile_temp)
  60 c      print *,"bumbum"
  61       do while (iret.gt.0)
  62
  63 #if (defined(AIX) && !defined(JUBL))
  64 #ifdef DEBUG
  65       write (iout,*) "ii",ii," itraj",itraj," it",it
  66 #endif
  67       call xdrffloat_(ixdrf, rtime, iret)
  68       call xdrffloat_(ixdrf, rpotE, iret)
  69 #ifdef DEBUG
  70       write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
  71 #endif
  72       call flush(iout)
  73       call xdrffloat_(ixdrf, ruconst, iret)
  74       call xdrffloat_(ixdrf, rt_bath, iret)
  75       call xdrfint_(ixdrf, nss, iret)
  76 #ifdef DEBUG
  77       write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
  78 #endif
  79       do j=1,nss
  80         call xdrfint_(ixdrf, ihpb(j), iret)
  81         call xdrfint_(ixdrf, jhpb(j), iret)
  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         call xdrfint(ixdrf, ihpb(j), iret)
 107         call xdrfint(ixdrf, jhpb(j), iret)
 108       enddo
 109       call xdrfint(ixdrf, nprop, iret)
 110 c      write (iout,*) "nprop",nprop
 111       if (it.gt.0 .and. nprop.ne.nprop_prev) then
 112         write (iout,*) "Warning previous nprop",nprop_prev,
 113      &   " current",nprop
 114         nprop=nprop_prev
 115       else
 116         nprop_prev=nprop
 117       endif
 118       call flush(iout)
 119       if (umbrella(iparm) .or. read_iset(iparm) .or. hamil_rep)
 120      &  call xdrfint(ixdrf, iset, iret)
 121       do i=1,nprop
 122         call xdrffloat(ixdrf, rprop(i), iret)
 123       enddo
 124 #endif
 125       if (iret.eq.0) exit
 126       itraj=mod(it,totraj(iR,iparm))
 127       if (iset.eq.0) iset = 1
 128       call flush(iout)
 129       it=it+1
 130       if (itraj.gt.ntraj) ntraj=itraj
 131       nstep(itraj)=nstep(itraj)+1
 132 c      rprop(2)=dsqrt(rprop(2))
 133 c      rprop(3)=dsqrt(rprop(3))
 134 #ifdef DEBUG
 135        write (iout,*) "umbrella ",umbrella
 136        write (iout,*) rtime,rpotE,rt_bath,nss,
 137      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
 138        write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
 139        call flush(iout)
 140 #endif
 141       prec=10000.0
 142
 143       itmp=0
 144 #if (defined(AIX) && !defined(JUBL))
 145       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 146 #else
 147       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
 148 #endif
 149 #ifdef DEBUG
 150       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
 151 #endif
 152       if (iret.eq.0) exit
 153       if (itmp .ne. nres + nct - nnt + 1) then
 154         write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
 155         call flush(iout)
 156         exit
 157       endif
 158
 159       time=rtime
 160 c      write (iout,*) "calling slice"
 161 c      call flush(iout)
 162       islice=slice(nstep(itraj),time,is,ie,ts,te)
 163 c      write (iout,*) "islice",islice
 164 c      call flush(iout)
 165
 166       do i=1,nres
 167         do j=1,3
 168           c(j,i)=xoord(j,i)
 169         enddo
 170       enddo
 171       do i=1,nct-nnt+1
 172         do j=1,3
 173           c(j,i+nres+nnt-1)=xoord(j,i+nres)
 174         enddo
 175       enddo
 176
 177       if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset
 178      &    .or. iset.eq.myparm)) then
 179         ii=ii+1
 180         kk(islice)=kk(islice)+1
 181         mm(islice)=mm(islice)+1
 182         if (mod(nstep(itraj),isampl(iparm)).eq.0 .and.
 183      &     conf_check(ll(islice)+1,1)) then
 184           if (replica(iparm)) then
 185              rt_bath=1.0d0/(rt_bath*1.987D-3)
 186              do i=1,nT_h(iparm)
 187                if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
 188                  iib = i
 189                  goto 22
 190                endif
 191              enddo
 192   22         continue
 193              if (i.gt.nT_h(iparm)) then
 194                write (iout,*) "Error - temperature of conformation",
 195      &         ii,1.0d0/(rt_bath*1.987D-3),
 196      &         " does not match any of the list"
 197                write (iout,*)
 198      &          1.0d0/(rt_bath*1.987D-3),
 199      &          (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 200                call flush(iout)
 201 c               exit
 202 c               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
 203                ii=ii-1
 204                kk(islice)=kk(islice)-1
 205                mm(islice)=mm(islice)-1
 206                goto 112
 207              endif
 208           else
 209             iib = ib
 210           endif
 211
 212           efree=0.0d0
 213           jj(islice)=jj(islice)+1
 214           if (umbrella(iparm)) then
 215             snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
 216           else if (hamil_rep) then
 217             snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
 218           else
 219             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
 220           endif
 221           ll(islice)=ll(islice)+1
 222 #ifdef DEBUG
 223           write (iout,*) "Writing conformation, record",ll(islice)
 224           write (iout,*) "ib",ib," iib",iib
 225           write (iout,*) "ntraj",ntraj," itraj",itraj,
 226      &      " nstep",nstep(itraj)
 227           write (iout,*) "pote",rpotE," time",rtime
 228 c          if (replica(iparm)) then
 229 c            write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
 230 c            write (iout,*) "TEMP list"
 231 c            write (iout,*)
 232 c     &       (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
 233 c          endif
 234           write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
 235 c          write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
 236 c          write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
 237           call flush(iout)
 238 #endif
 239           if (islice.ne.islice1) then
 240 c            write (iout,*) "islice",islice," islice1",islice1
 241             close(ientout)
 242 c            write (iout,*) "Closing file ",
 243 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 244             call opentmp(islice,ientout,bprotfile_temp)
 245 c            write (iout,*) "Opening file ",
 246 c     &          bprotfile_temp(:ilen(bprotfile_temp))
 247             islice1=islice
 248           endif
 249           if (umbrella(iparm)) then
 250             write(ientout,rec=ll(islice))
 251      &        ((xoord(l,k),l=1,3),k=1,nres),
 252      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 253      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 254      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 255      &        iset,iib,iparm
 256           else if (hamil_rep) then
 257             write(ientout,rec=ll(islice))
 258      &        ((xoord(l,k),l=1,3),k=1,nres),
 259      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 260      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 261      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 262      &        iR,iib,iset
 263           else
 264             write(ientout,rec=ll(islice))
 265      &        ((xoord(l,k),l=1,3),k=1,nres),
 266      &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
 267      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 268      &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
 269      &        iR,iib,iparm
 270           endif
 271 #ifdef PDBOUT
 272 #ifdef MPI
 273           if (me.eq.Master) then
 274 #endif
 275           write (iout,*) "PDBOUT"
 276           write (iout,*) "temperature",1.0d0/(rt_bath*1.987D-3)
 277           call flush(iout)
 278           write (lt_bath,'(f4.0)') 1.0d0/(rt_bath*1.987D-3)
 279           write (iout,*) "lt_bath ",lt_bath
 280           pdbfilename=prefix(:ilen(prefix))//"_"//lt_bath//"pdb"
 281           write (iout,*) "pdb ",pdbfilename
 282           call flush(iout)
 283           open(ipdb,file=pdbfilename,position="append")
 284 c          write (tytul,'("Conformation",i10," T=",f5.1)')
 285 c     &      kk(islice),rt_bath
 286           call pdbout(kk(islice),1.0d0/(rt_bath*1.987D-3),
 287      &      efree+0.0d0,rpotE+0.0d0,efree+0.0d0,rmsdev+0.0d0)
 288           close(ipdb)
 289 #ifdef MPI
 290           endif
 291 #endif
 292 #endif
 293 #ifdef DEBUG
 294           call int_from_cart1(.false.)
 295           write (iout,*) "Writing conformation, record",ll(islice)
 296           write (iout,*) "Cartesian coordinates"
 297           write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
 298           write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
 299           write (iout,*) "Internal coordinates"
 300           write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
 301           write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
 302           write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
 303           write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
 304           write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
 305           write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
 306           write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
 307 c          write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
 308           write (iout,'(16i5)') iscor
 309           call flush(iout)
 310 #endif
 311         endif
 312       endif
 313
 314   112 continue
 315
 316       enddo
 317       close(ientout)
 318 #if (defined(AIX) && !defined(JUBL))
 319       call xdrfclose_(ixdrf, iret)
 320 #else
 321       call xdrfclose(ixdrf, iret)
 322 #endif
 323       write (iout,'(i10," trajectories found in file.")') ntraj+1
 324       write (iout,'(a)') "Numbers of steps in trajectories:"
 325       write (iout,'(8i10)') (nstep(i),i=0,ntraj)
 326       write (iout,*) ii," conformations read from file",
 327      &   nazwa(:ilen(nazwa))
 328       do islice=1,nslice
 329         write (iout,*) mm(islice)," conformations read so far, slice",
 330      &    islice
 331         write (iout,*) ll(islice),
 332      &  " conformations stored so far, slice",islice
 333       enddo
 334       call flush(iout)
 335       return
 336       end