Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / wham / io_database.F90
index c072b2a..94a2957 100644 (file)
           print *,"Calling cxread"
           call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,&
              *1111)
+           print *,"after call cxread"
 write(iout,*)"after call cxread"
           close(ientout)
           write (iout,*) "exit cxread"
@@ -441,7 +442,7 @@ write(iout,*) "end of read database"
 !      rtime=0.0d0
 !      rpotE=0.0d0
 !      rt_bath=0.0d0
-
+      rmsdev=0.0d0
       call set_slices(is,ie,ts,te,iR,ib,iparm)
       nprop_prev=0
       do i=1,nQ
@@ -486,12 +487,14 @@ write(iout,*) "end of read database"
       enddo
 #else
       call xdrffloat(ixdrf, rtime, iret)
+      print *,"rtime",rtime," iret",iret !d
       call xdrffloat(ixdrf, rpotE, iret)
 !      write (iout,*) "rpotE",rpotE," iret",iret !d
-      call flush(iout)
+!      call flush(iout)
       call xdrffloat(ixdrf, ruconst, iret)
       call xdrffloat(ixdrf, rt_bath, iret)
       call xdrfint(ixdrf, nss, iret)
+      print *,"nss",nss
       do j=1,nss
         call xdrfint(ixdrf, ihpb(j), iret)
         call xdrfint(ixdrf, jhpb(j), iret)
@@ -512,6 +515,7 @@ write(iout,*) "end of read database"
         call xdrffloat(ixdrf, rprop(i), iret)
       enddo
 #endif
+      print *,"iret",iret
       if (iret.eq.0) exit
       itraj=mod(it,totraj(iR,iparm))
 !#define DEBUG
@@ -537,7 +541,9 @@ write(iout,*) "end of read database"
 #if (defined(AIX) && !defined(JUBL))
       call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
 #else
+      print *,"before xdrf3dcoord"
       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
+      print *,"after xdrf3dcoord", iret
 #endif
 #ifdef DEBUG
       write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,2*nres+2)
@@ -704,6 +710,7 @@ write(iout,*) "end of read database"
         " conformations stored so far, slice",islice
       enddo
       call flush(iout)
+      print *,"before cxread return"
 !#undef DEBUG
       return
       end subroutine cxread
@@ -910,7 +917,8 @@ write(iout,*) "end of read database"
       use geometry_data
       use control_data, only:indpdb
       use w_compar_data
-      use conform_compar, only:conf_compar
+      use conform_compar, only:conf_compar,rmsnat,qwolynes
+      use energy_data, only:icont,ncont,nnt,nct,maxcont!,&
 !      implicit none
 !      include "DIMENSIONS"
 !      include "DIMENSIONS.ZSCOPT"
@@ -948,13 +956,15 @@ write(iout,*) "end of read database"
       integer :: i,itj,ii,iii,j,k,l
       integer :: ixdrf,iret
       integer :: iscor,islice
-      real(kind=8) :: rmsdev,efree,eini
+      real(kind=8) :: rmsdev,efree,eini,qnat2
       real(kind=4) :: csingle(3,nres*2)
       real(kind=8) :: energ
+       
 !      integer ilen,iroof
 !      external ilen,iroof
       integer :: ir,ib,iparm
       integer :: isecstr(nres)
+      logical :: test
       write (licz2,'(bz,i2.2)') islice
       call opentmp(islice,ientout,bprotfile_temp)
       write (iout,*) "bprotfile_temp ",bprotfile_temp
@@ -973,6 +983,7 @@ write(iout,*) "end of read database"
           write (licz,'(bz,i3.3)') myparm
           bxname = prefix(:ilen(prefix))//liczba//"_par"//licz//".bx"
         endif
+        print *,bxname
         open (ientin,file=bxname,status="unknown",&
           form="unformatted",access="direct",recl=lenrec1)
       endif
@@ -1000,6 +1011,7 @@ write(iout,*) "end of read database"
         endif
 !el      endif 
 #endif
+      print *,indpdb
       if (indpdb.gt.0) then
         if (nslice.eq.1) then
 #ifdef MPI
@@ -1011,7 +1023,7 @@ write(iout,*) "end of read database"
            statname=prefix(:ilen(prefix))//'_par'//licz//'_'// &
             pot(:ilen(pot))//liczba//'.stat'
          endif
-
+         print *,statname
 #else
           statname=prefix(:ilen(prefix))//'_'//pot(:ilen(pot))//'.stat'
 #endif
@@ -1030,14 +1042,17 @@ write(iout,*) "end of read database"
             //"_slice_"//licz2//'.stat'
 #endif
         endif
+        print *,istat,statname
         open(istat,file=statname,status="unknown")
       endif
-
+      print *,"Tu dochodze"
+      print *,scount(me)
 #ifdef MPI
       do i=1,scount(me)
 #else
       do i=1,ntot(islice)
 #endif
+        print *,"before ientout read"
         read(ientout,rec=i,err=101) &
           ((csingle(l,k),l=1,3),k=1,nres),&
           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
@@ -1054,14 +1069,24 @@ write(iout,*) "end of read database"
 !        write (iout,*) "Calling conf_compar",i
 !        call flush(iout)
          anatemp= 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+        print *,"before conf_compar"
         if (indpdb.gt.0) then
-          call conf_compar(i,.false.,.true.)
+        print *,"just before conf_compar",i
+        print *,icont,ncont,nnt,nct,"maxcont",maxcont
+        test=.false.
+!          call conf_compar(i,.false.,.true.)
+!          call conf_compar(i)
+!           call rmsnat(i)
+           rms_nat=rmsnat(i)
+           qnat2=qwolynes(0,0) 
+         print *,"just after conf_compar"
 !        else
 !            call elecont(.false.,ncont,icont,nnt,nct)
 !            call secondary2(.false.,.false.,ncont,icont,isecstr)
         endif
 !        write (iout,*) "Exit conf_compar",i
 !        call flush(iout)
+         print *,"before ientin"
         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i) &
           ((csingle(l,k),l=1,3),k=1,nres),&
           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
@@ -1078,6 +1103,7 @@ write(iout,*) "end of read database"
       close(istat)
       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 #ifdef MPI
+      print *,"before MPI_barrier"
       call MPI_Barrier(WHAM_COMM,IERROR)
       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile &
          .and. ensembles.eq.0) return
@@ -1159,7 +1185,9 @@ write(iout,*) "end of read database"
               nss,(ihpb(k),jhpb(k),k=1,nss),&
               eini,efree,rmsdev,iscor
           endif
+!           print *,"before cxwrite"
           if(cxfile)call cxwrite(ixdrf,csingle,eini,efree,rmsdev,iscor)
+!          print *,"after cxwrite"
 #ifdef DEBUG
           do k=1,2*nres
             do l=1,3
@@ -1257,24 +1285,34 @@ write(iout,*) "end of read database"
 !      write (iout,*) "xdrf3dfcoord"
 !      call flush(iout)
       call xdrfint_(ixdrf, nss, iret)
+            write (iout,*) "iret",iret
+            write (iout,*) "nss",nss,i,"TUTU"
       do j=1,nss
         call xdrfint_(ixdrf, ihpb(j), iret)
         call xdrfint_(ixdrf, jhpb(j), iret)
+            write(iout,*), ihpb(j),jhpb(j),"TUTU"
       enddo
       call xdrffloat_(ixdrf,real(eini),iret) 
       call xdrffloat_(ixdrf,real(efree),iret) 
+            write(iout,*) "TUTU", eini
+            write(iout,*) "TUTU", efree
       call xdrffloat_(ixdrf,real(rmsdev),iret) 
       call xdrfint_(ixdrf,iscor,iret) 
 #else
       call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
+            write (iout,*) "iret",iret
+            write (iout,*) "nss",nss,i,"TUTU"
 
       call xdrfint(ixdrf, nss, iret)
       do j=1,nss
         call xdrfint(ixdrf, ihpb(j), iret)
         call xdrfint(ixdrf, jhpb(j), iret)
+            write(iout,*), ihpb(j),jhpb(j),"TUTU"
       enddo
       call xdrffloat(ixdrf,real(eini),iret) 
       call xdrffloat(ixdrf,real(efree),iret) 
+            write(iout,*) "TUTU", eini
+            write(iout,*) "TUTU", efree
       call xdrffloat(ixdrf,real(rmsdev),iret) 
       call xdrfint(ixdrf,iscor,iret) 
 #endif
@@ -1296,7 +1334,7 @@ write(iout,*) "end of read database"
       integer :: islice,iR,ib,iparm
       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
-
+      time_slice=0
       do islice=1,nslice
         if (time_end_collect(iR,ib,iparm).ge.1.0d10) then
           ts(islice)=time_start_collect(iR,ib,iparm)
@@ -1464,6 +1502,11 @@ write(iout,*) "end of read database"
         endif
       enddo
       do j=3,nres
+        mnum=molnum(j)
+        itj=itype(j,mnum)
+        if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
+        if (itype(j-1,mnum).eq.ntyp1_molec(mnum)) cycle
+        if (itype(j-2,mnum).eq.ntyp1_molec(mnum)) cycle
         if (theta(j).le.0.0d0) then
           if (iprint.gt.0) &
           write (iout,*) "Zero theta angle(s) in conformation",ii