Merge branch 'devel' into AFM
[unres.git] / source / wham / src-NEWSC-NEWCORR / cxread.F.org
diff --git a/source/wham/src-NEWSC-NEWCORR/cxread.F.org b/source/wham/src-NEWSC-NEWCORR/cxread.F.org
new file mode 100644 (file)
index 0000000..80bc1a0
--- /dev/null
@@ -0,0 +1,248 @@
+      subroutine cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,*)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
+      integer MaxTraj
+      parameter (MaxTraj=2050)
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.HEADER'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.PROTFILES'
+      include 'COMMON.OBCINKA'
+      include 'COMMON.FREE'
+      include 'COMMON.VAR'
+      include 'COMMON.GEO'
+      include 'COMMON.PROT'
+      character*64 nazwa,bprotfile_temp
+      real*4 rtime,rpotE,ruconst,rt_bath,rprop(maxQ)
+      double precision time
+      integer iret,itmp,itraj,ntraj
+      real xoord(3,maxres2+2),prec
+      integer nstep(0:MaxTraj-1)
+      integer ilen
+      external ilen
+      integer ii,jj(maxslice),kk(maxslice),ll(maxslice),mm(maxslice)
+      integer is(MaxSlice),ie(MaxSlice),nrec_slice
+      double precision ts(MaxSlice),te(MaxSlice),time_slice
+      integer slice
+      call set_slices(is,ie,ts,te,iR,ib,iparm)
+
+      do i=1,nQ
+        rprop(i)=0.0d0
+      enddo
+      do i=0,MaxTraj-1
+        nstep(i)=0
+      enddo
+      ntraj=0
+      it=0
+      iret=1
+#if (defined(AIX) && !defined(JUBL))
+      call xdrfopen_(ixdrf,nazwa, "r", iret)
+#else
+      call xdrfopen(ixdrf,nazwa, "r", iret)
+#endif
+      if (iret.eq.0) return1
+
+      islice1=1
+      call opentmp(islice1,ientout,bprotfile_temp)
+c      print *,"bumbum"
+      do while (iret.gt.0) 
+
+#if (defined(AIX) && !defined(JUBL))
+      call xdrffloat_(ixdrf, rtime, iret)
+c      print *,"rtime",rtime," iret",iret
+      call xdrffloat_(ixdrf, rpotE, iret)
+c      write (iout,*) "rpotE",rpotE," iret",iret
+      call flush(iout)
+      call xdrffloat_(ixdrf, ruconst, iret)
+      call xdrffloat_(ixdrf, rt_bath, iret)
+      call xdrfint_(ixdrf, nss, iret)
+      do j=1,nss
+        call xdrfint_(ixdrf, ihpb(j), iret)
+        call xdrfint_(ixdrf, jhpb(j), iret)
+      enddo
+      call xdrfint_(ixdrf, nprop, iret)
+      do i=1,nprop
+        call xdrffloat_(ixdrf, rprop(i), iret)
+      enddo
+#else
+      call xdrffloat(ixdrf, rtime, iret)
+      call xdrffloat(ixdrf, rpotE, iret)
+c      write (iout,*) "rpotE",rpotE," iret",iret
+      call flush(iout)
+      call xdrffloat(ixdrf, ruconst, iret)
+      call xdrffloat(ixdrf, rt_bath, iret)
+      call xdrfint(ixdrf, nss, iret)
+      do j=1,nss
+        call xdrfint(ixdrf, ihpb(j), iret)
+        call xdrfint(ixdrf, jhpb(j), iret)
+      enddo
+      call xdrfint(ixdrf, nprop, iret)
+c      write (iout,*) "nprop",nprop
+      call flush(iout)
+      do i=1,nprop
+        call xdrffloat(ixdrf, rprop(i), iret)
+      enddo
+#endif
+      if (iret.eq.0) exit
+      itraj=mod(it,totraj(iR,iparm))
+#ifdef DEBUG
+      write (iout,*) "ii",ii," itraj",itraj
+#endif
+      call flush(iout)
+      it=it+1
+      if (itraj.gt.ntraj) ntraj=itraj
+      nstep(itraj)=nstep(itraj)+1
+#ifdef DEBUG
+       write (iout,*) rtime,rpotE,rt_bath,nss,
+     &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
+       call flush(iout)
+#endif
+      prec=10000.0
+
+      itmp=0
+#if (defined(AIX) && !defined(JUBL))
+      call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
+#else
+      call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
+#endif
+#ifdef DEBUG
+      write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
+#endif
+      if (iret.eq.0) exit
+      if (itmp .ne. nres + nct - nnt + 1) then
+        write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
+        call flush(iout)
+        exit
+      endif
+
+      time=rtime
+c      write (iout,*) "calling slice"
+c      call flush(iout)
+      islice=slice(nstep(itraj),time,is,ie,ts,te)
+c      write (iout,*) "islice",islice
+c      call flush(iout)
+
+      if (islice.gt.0 .and. islice.le.nslice) then
+        ii=ii+1
+        kk(islice)=kk(islice)+1
+        mm(islice)=mm(islice)+1
+        if (mod(nstep(itraj),isampl(iparm)).eq.0) then
+            if (replica(iparm)) then
+               rt_bath=1.0d0/(rt_bath*1.987D-3)
+               do i=1,nT_h(iparm)
+                 if (abs(real(beta_h(i,iparm))-rt_bath).lt.1.0e-4) then
+                   iib = i
+                   goto 22
+                 endif
+               enddo
+  22           continue
+               if (i.gt.nT_h(iparm)) then
+                 write (iout,*) "Error - temperature of conformation",
+     &           ii,1.0d0/(rt_bath*1.987D-3),
+     &           " does not match any of the list"
+                 write (iout,*)
+     &            1.0d0/(rt_bath*1.987D-3),
+     &            (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
+                 call flush(iout)
+                 exit
+                 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+               endif
+            else
+                iib = ib
+            endif
+
+            efree=0.0d0
+            jj(islice)=jj(islice)+1
+            snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
+            ll(islice)=ll(islice)+1
+#ifdef DEBUG
+            write (iout,*) "Writing conformation, record",ll(islice)
+            write (iout,*) "ib",ib," iib",iib
+            write (iout,*) "ntraj",ntraj," itraj",itraj,
+     &        " nstep",nstep(itraj)
+            write (iout,*) "pote",rpotE," time",rtime
+c            if (replica(iparm)) then
+c              write (iout,*) "TEMP",1.0d0/(rt_bath*1.987D-3)
+c              write (iout,*) "TEMP list"
+c              write (iout,*)
+c     &         (1.0d0/(beta_h(i,iparm)*1.987D-3),i=1,nT_h(iparm))
+c            endif
+            write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
+c            write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
+c            write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
+            call flush(iout)
+#endif
+            if (islice.ne.islice1) then
+c              write (iout,*) "islice",islice," islice1",islice1
+              close(ientout) 
+c              write (iout,*) "Closing file ",
+c     &            bprotfile_temp(:ilen(bprotfile_temp))
+              call opentmp(islice,ientout,bprotfile_temp)
+c              write (iout,*) "Opening file ",
+c     &            bprotfile_temp(:ilen(bprotfile_temp))
+              islice1=islice
+            endif
+            write(ientout,rec=ll(islice))
+     &        ((xoord(l,k),l=1,3),k=1,nres),
+     &        ((xoord(l,k),l=1,3),k=nres+1,nres+nct-nnt+1),
+     &        nss,(ihpb(k),jhpb(k),k=1,nss),
+     &        rpotE+0.0d0,efree,rmsdev,(rprop(i)+0.0d0,i=1,nQ),
+     &        iR,iib,iparm
+#ifdef DEBUG
+            do i=1,nres
+              do j=1,3
+                c(j,i)=xoord(j,i)
+              enddo
+            enddo
+            do i=1,nct-nnt+1
+              do j=1,3
+                c(j,i+nres+nnt-1)=xoord(j,i+nres)
+              enddo
+            enddo
+            call int_from_cart1(.false.)
+            write (iout,*) "Writing conformation, record",ll(islice)
+            write (iout,*) "Cartesian coordinates"
+            write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
+            write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
+            write (iout,*) "Internal coordinates"
+            write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
+            write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
+            write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+            write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+            write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+            write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+            write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
+c            write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
+            write (iout,'(16i5)') iscor
+            call flush(iout)
+#endif
+        endif 
+      endif
+
+      enddo
+  112 continue
+      close(ientout)
+#if (defined(AIX) && !defined(JUBL))
+      call xdrfclose_(ixdrf, iret)
+#else
+      call xdrfclose(ixdrf, iret)
+#endif
+      write (iout,'(i10," trajectories found in file.")') ntraj+1
+      write (iout,'(a)') "Numbers of steps in trajectories:"
+      write (iout,'(8i10)') (nstep(i),i=0,ntraj)
+      write (iout,*) ii," conformations read from file",
+     &   nazwa(:ilen(nazwa))
+      do islice=1,nslice
+        write (iout,*) mm(islice)," conformations read so far, slice",
+     &    islice
+        write (iout,*) ll(islice),
+     &  " conformations stored so far, slice",islice
+      enddo
+      call flush(iout)
+      return
+      end