read2sigma in wham and cluster_wham
[unres.git] / source / wham / src / cxread.F.org
index 80bc1a0..7c7ae50 100644 (file)
@@ -5,6 +5,7 @@
       include 'DIMENSIONS.FREE'
       integer MaxTraj
       parameter (MaxTraj=2050)
+      include 'COMMON.CONTROL'
       include 'COMMON.CHAIN'
       include 'COMMON.INTERACT'
       include 'COMMON.NAMES'
@@ -29,6 +30,7 @@
       integer is(MaxSlice),ie(MaxSlice),nrec_slice
       double precision ts(MaxSlice),te(MaxSlice),time_slice
       integer slice
+      logical conf_check
       call set_slices(is,ie,ts,te,iR,ib,iparm)
 
       do i=1,nQ
@@ -53,53 +55,98 @@ c      print *,"bumbum"
       do while (iret.gt.0) 
 
 #if (defined(AIX) && !defined(JUBL))
+#ifdef DEBUG
+      write (iout,*) "ii",ii," itraj",itraj," it",it
+#endif
       call xdrffloat_(ixdrf, rtime, iret)
-c      print *,"rtime",rtime," iret",iret
       call xdrffloat_(ixdrf, rpotE, iret)
-c      write (iout,*) "rpotE",rpotE," iret",iret
+#ifdef DEBUG
+      write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
+#endif
       call flush(iout)
       call xdrffloat_(ixdrf, ruconst, iret)
       call xdrffloat_(ixdrf, rt_bath, iret)
       call xdrfint_(ixdrf, nss, iret)
+#ifdef DEBUG
+      write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
+#endif
       do j=1,nss
+       if (dyn_ss) then
+        call xdrfint_(ixdrf, idssb(j), iret)
+        call xdrfint_(ixdrf, jdssb(j), iret)
+       idssb(j)=idssb(j)-nres
+       jdssb(j)=jdssb(j)-nres
+       else
         call xdrfint_(ixdrf, ihpb(j), iret)
         call xdrfint_(ixdrf, jhpb(j), iret)
+       endif
       enddo
       call xdrfint_(ixdrf, nprop, iret)
+      if (umbrella(iparm) .or. homol_nset.gt.1 .or. read_iset(iparm) 
+     &  .or. hamil_rep) 
+     &  call xdrfint(ixdrf, iset, iret)
       do i=1,nprop
         call xdrffloat_(ixdrf, rprop(i), iret)
       enddo
 #else
+#ifdef DEBUG
+      write (iout,*) "ii",ii," itraj",itraj," it",it
+#endif
       call xdrffloat(ixdrf, rtime, iret)
       call xdrffloat(ixdrf, rpotE, iret)
-c      write (iout,*) "rpotE",rpotE," iret",iret
+#ifdef DEBUG
+      write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
+#endif
       call flush(iout)
       call xdrffloat(ixdrf, ruconst, iret)
       call xdrffloat(ixdrf, rt_bath, iret)
       call xdrfint(ixdrf, nss, iret)
+#ifdef DEBUG
+      write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
+#endif
       do j=1,nss
+       if (dyn_ss) then
+        call xdrfint(ixdrf, idssb(j), iret)
+        call xdrfint(ixdrf, jdssb(j), iret)
+cc        idssb(j)=idssb(j)-nres
+cc        jdssb(j)=jdssb(j)-nres
+cc        write(iout,*) idssb(j),jdssb(j)
+       else
         call xdrfint(ixdrf, ihpb(j), iret)
         call xdrfint(ixdrf, jhpb(j), iret)
+       endif
       enddo
       call xdrfint(ixdrf, nprop, iret)
 c      write (iout,*) "nprop",nprop
+      if (it.gt.0 .and. nprop.ne.nprop_prev) then
+        write (iout,*) "Warning previous nprop",nprop_prev,
+     &   " current",nprop
+        nprop=nprop_prev
+      else
+        nprop_prev=nprop
+      endif
       call flush(iout)
+      if (umbrella(iparm) .or. homol_nset.gt.1 .or. read_iset(iparm) 
+     &  .or. hamil_rep) 
+     &  call xdrfint(ixdrf, iset, iret)
       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
+      if (iset.eq.0) iset = 1
       call flush(iout)
       it=it+1
       if (itraj.gt.ntraj) ntraj=itraj
       nstep(itraj)=nstep(itraj)+1
+c      rprop(2)=dsqrt(rprop(2))
+c      rprop(3)=dsqrt(rprop(3))
 #ifdef DEBUG
+       write (iout,*) "umbrella ",umbrella
        write (iout,*) rtime,rpotE,rt_bath,nss,
      &     (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
+       write (iout,*) "nprop",nprop," iset",iset," myparm",myparm
        call flush(iout)
 #endif
       prec=10000.0
@@ -127,105 +174,137 @@ c      call flush(iout)
 c      write (iout,*) "islice",islice
 c      call flush(iout)
 
-      if (islice.gt.0 .and. islice.le.nslice) then
+      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
+
+      if (islice.gt.0 .and. islice.le.nslice .and. (.not.separate_parset
+     &    .or. iset.eq.myparm)) 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)
+        if (mod(nstep(itraj),isampl(iparm)).eq.0 .and. 
+     &     conf_check(ll(islice)+1,1)) 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
-            else
-                iib = ib
-            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)
+c               exit
+c               call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
+               ii=ii-1
+               kk(islice)=kk(islice)-1
+               mm(islice)=mm(islice)-1
+               goto 112
+             endif
+          else
+            iib = ib
+          endif
 
-            efree=0.0d0
-            jj(islice)=jj(islice)+1
+          efree=0.0d0
+          jj(islice)=jj(islice)+1
+          if (umbrella(iparm)) then
+            snk(iset,iib,iparm,islice)=snk(iset,iib,iparm,islice)+1
+          else if (hamil_rep) then
+            snk(1,iib,iparm,islice)=snk(1,iib,iparm,islice)+1
+          else
             snk(iR,iib,iparm,islice)=snk(iR,iib,iparm,islice)+1
-            ll(islice)=ll(islice)+1
+          endif
+          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)
+          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
+          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
+          if (umbrella(iparm) .or. homol_nset.gt.1) then
+            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),
+     &        iset,iib,iparm
+          else if (hamil_rep) then
+            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,iset
+          else
             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
+          endif
 #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)
+          write (iout,*) " constr_homology",constr_homology,
+     &      " ll",ll(islice)," iset",iset
+          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
+
+      enddo
       close(ientout)
 #if (defined(AIX) && !defined(JUBL))
       call xdrfclose_(ixdrf, iret)