rename
[unres4.git] / source / cluster / track.F90
diff --git a/source/cluster/track.F90 b/source/cluster/track.F90
new file mode 100644 (file)
index 0000000..542da6a
--- /dev/null
@@ -0,0 +1,306 @@
+      module tracking
+!------------------------------------------------------------------------------
+      use clust_data
+      implicit none
+!------------------------------------------------------------------------------
+!      COMMON /HISTORY/
+      integer :: NCUR(MAX_CUT),IBACK(MAXGR,MAX_CUT)
+!      COMMON /PREVIOUS/
+      integer :: NGRP,LICZP(MAXGR),NCONFP(MAXGR,MAXINGR) 
+!------------------------------------------------------------------------------
+!
+!
+!------------------------------------------------------------------------------
+      contains
+!------------------------------------------------------------------------------
+      SUBROUTINE TRACK(ICUT)
+!      include 'DIMENSIONS'
+!      INCLUDE 'sizesclu.dat'
+!      INCLUDE 'COMMON.CLUSTER'
+!      COMMON /HISTORY/ NCUR(MAX_CUT),IBACK(MAXGR,MAX_CUT)
+!      COMMON /PREVIOUS/ NGRP,LICZP(MAXGR),NCONFP(MAXGR,MAXINGR)
+      integer :: icut,igr,jgr,k,nci1
+      IF (ICUT.GT.1) THEN
+! Find out what of the previous families the current ones came from.        
+        DO IGR=1,NGR
+          NCI1=NCONF(IGR,1) 
+          DO JGR=1,NGRP
+            DO K=1,LICZP(JGR)
+              IF (NCI1.EQ.NCONFP(JGR,K)) THEN
+                IBACK(IGR,ICUT)=JGR
+                GOTO 10
+              ENDIF
+            ENDDO ! K
+          ENDDO ! JGR
+  10      CONTINUE
+        ENDDO ! IGR
+      ENDIF ! (ICUT.GT.1)
+! Save current partition for subsequent backtracking.
+      NCUR(ICUT)=NGR
+      NGRP=NGR
+      DO IGR=1,NGR
+        LICZP(IGR)=LICZ(IGR)
+        DO K=1,LICZ(IGR)
+          NCONFP(IGR,K)=NCONF(IGR,K)
+        ENDDO ! K
+      ENDDO ! IGR
+      RETURN
+      END SUBROUTINE TRACK
+!------------------------------------------------------------------------------
+      SUBROUTINE WRITRACK
+
+      use io_units, only: iout
+!      include 'DIMENSIONS'
+!      INCLUDE 'sizesclu.dat'
+!      INCLUDE 'COMMON.CLUSTER'
+!      include 'COMMON.IOUNITS'
+!      COMMON /HISTORY/ NCUR(MAX_CUT),IBACK(MAXGR,MAX_CUT)
+!      COMMON /PREVIOUS/ NGRP,LICZP(MAXGR),NCONFP(MAXGR,MAXINGR)
+      integer :: IPART(MAXGR/5,MAXGR/5)
+      integer :: icut,i,j,k,ncu,ncup,npart
+!     do icut=2,ncut
+!       write (iout,'(a,f10.5)') 'Cut-off',rcutoff(icut)
+!       write (iout,'(16i5)') (iback(k,icut),k=1,ncur(icut))
+!     enddo
+!
+! Print the partition history.
+!
+      DO ICUT=2,NCUT
+        NCU=NCUR(ICUT)
+        NCUP=NCUR(ICUT-1)
+!d      print *,'icut=',icut,' ncu=',ncu,' ncur=',ncur
+        WRITE(iout,'(A,f10.5,A,f10.5)') &
+        'Partition of families obtained at cut-off',RCUTOFF(ICUT-1),&
+        ' at cut-off',RCUTOFF(ICUT)
+        DO I=1,NCUP
+          NPART=0
+!d        print *,'i=',i
+          DO J=1,NCU
+            IF (IBACK(J,ICUT).EQ.I) THEN
+              NPART=NPART+1  
+              IPART(NPART,I)=J
+            ENDIF
+!d          print *,'j=',j,' iback=',IBACK(J,ICUT),' npart=',npart
+          ENDDO ! J
+          WRITE (iout,'(16I5)') I,(IPART(K,I),K=1,NPART) 
+        ENDDO ! I
+      ENDDO ! ICUT
+      RETURN
+      END SUBROUTINE WRITRACK
+!------------------------------------------------------------------------------
+      SUBROUTINE PLOTREE
+
+      use io_units, only: jplot
+      use io_base, only: ilen
+!      include 'DIMENSIONS'
+!      INCLUDE 'sizesclu.dat'
+!      INCLUDE 'COMMON.CLUSTER'
+!      include 'COMMON.IOUNITS'
+!      COMMON /HISTORY/ NCUR(MAX_CUT),IBACK(MAXGR,MAX_CUT)
+!      COMMON /PREVIOUS/ NGRP,LICZP(MAXGR),NCONFP(MAXGR,MAXINGR)
+      integer,DIMENSION(MAXGR,MAX_CUT) :: Y
+      integer,DIMENSION(MAXGR,MAX_CUT) :: ITREE,IFIRST,ILAST
+      integer,dimension(MAXGR) :: IFT,ILT,ITR
+      CHARACTER(len=32) :: FD
+      integer :: i,icut,j,k,is,iti,jf1,jf2,jl1,jl2,ncu,ncup,npart
+      integer :: jr1,jr2,jf11,kl,kf,nnc,iylen,ideltx,idelty
+      real(kind=8) :: xlen,ylen,xbox,ybox,deltx,yy
+      real(kind=8) :: ycur,xcur,xdraw,ydraw,delty
+!el      external ilen
+! 
+! Generate the image of the tree (tentatively for LaTeX picture environment).
+!
+!
+! First untangle the branches of the tree
+!
+      DO I=1,NCUR(1)
+        ITREE(I,1)=I
+      ENDDO
+      DO ICUT=NCUT,2,-1
+!
+! Determine the order of families for the (icut)th partition.
+!
+        NCU=NCUR(ICUT)
+        NCUP=NCUR(ICUT-1)
+        NPART=0
+        DO I=1,NCUP
+          IS=0
+          IF (I.GT.1) ILAST(I-1,ICUT-1)=NPART
+          DO J=1,NCU
+            IF (IBACK(J,ICUT).EQ.I) THEN
+              NPART=NPART+1  
+              IF (IS.EQ.0) THEN
+                IS=1
+                IFIRST(I,ICUT-1)=NPART
+              ENDIF
+              ITREE(NPART,ICUT)=J
+            ENDIF
+          ENDDO ! J
+        ENDDO ! I
+        ILAST(NCUP,ICUT-1)=NPART
+!d      print *,'i=',i,' ncup=',ncup,' ncu=',ncu,' npart=',npart
+      ENDDO ! ICUT
+! diagnostic printout
+!d    do icut=1,ncut
+!d      write (iout,*) 'Cut-off',icut,' = ',rcutoff(icut) 
+!d      write (iout,*) 'ITREE'
+!d      write (iout,*) (itree(i,icut),i=1,ncur(icut))
+!d      write (iout,*) 'IFIRST, ILAST'
+!d      write (iout,*) (ifirst(i,icut),ilast(i,icut),i=1,ncur(icut))
+!d    enddo
+!
+! Propagate the order of families from cut-off #2 to cut-off #n.
+!
+      DO ICUT=1,NCUT-1
+        DO J=1,NCUR(ICUT)
+          IFT(J)=IFIRST(J,ICUT)
+          ILT(J)=ILAST(J,ICUT)
+        ENDDO ! J
+        DO J=1,NCUR(ICUT+1)
+          ITR(J)=ITREE(J,ICUT+1)
+        ENDDO
+        DO I=1,NCUR(ICUT)
+          ITI=ITREE(I,ICUT)
+!         write (iout,*) 'icut=',icut,' i=',i,' iti=',iti
+!         IF (ITI.NE.I) THEN
+            JF1=IFT(I)
+            JF2=IFT(ITI)
+            JL1=ILT(I)
+            JL2=ILT(ITI)
+            JR1=JL1-JF1+1
+            JR2=JL2-JF2+1
+!d          write (iout,*) 'jf1=',jf1,' jl1=',jl1,' jf2=',jf2,
+!d   &                     ' jl2=',jl2
+!d          write (iout,*) 'jr1=',jr1,' jr2=',jr2
+! Update IFIRST and ILAST.
+            ILAST(I,ICUT)=IFIRST(I,ICUT)+JR2-1
+            IFIRST(I+1,ICUT)=ILAST(I,ICUT)+1
+! Update ITREE.
+            JF11=IFIRST(I,ICUT)
+!d          write(iout,*) 'jf11=',jf11
+            DO J=JF2,JL2
+!d            write (iout,*) j,JF11+J-JF2,ITR(J)
+              ITREE(JF11+J-JF2,ICUT+1)=ITR(J)
+            ENDDO
+!d      write (iout,*) (ifirst(k,icut),ilast(k,icut),k=1,i)
+!d      write (iout,*) (itree(k,icut+1),k=1,ilast(i,icut))
+!         ENDIF ! (ITI.NE.I)
+        ENDDO ! I
+      ENDDO ! ICUT
+! diagnostic printout
+!d    do icut=1,ncut
+!d      write (iout,*) 'Cut-off',icut,' = ',rcutoff(icut) 
+!d      write (iout,*) 'ITREE'
+!d      write (iout,*) (itree(i,icut),i=1,ncur(icut))
+!d      write (iout,*) 'IFIRST, ILAST'
+!d      write (iout,*) (ifirst(i,icut),ilast(i,icut),i=1,ncur(icut))
+!d    enddo
+!
+! Generate the y-coordinates of the branches.
+!
+      XLEN=400.0/(ncut-1)
+      YLEN=600.0
+      xbox=xlen/4.0
+      deltx=0.5*(xlen-xbox)
+      NNC=NCUR(NCUT)
+      ybox=ylen/(2.0*nnc) 
+      DO J=1,NNC
+        Y(J,NCUT)=J*YLEN/NNC
+      ENDDO
+      DO ICUT=NCUT-1,1,-1
+        NNC=NCUR(ICUT)
+        DO J=1,NNC
+          KF=IFIRST(J,ICUT)
+          KL=ILAST(J,ICUT)
+          YY=0.0
+          DO K=KF,KL
+            YY=YY+Y(K,ICUT+1)
+          ENDDO
+          Y(J,ICUT)=YY/(KL-KF+1)
+        ENDDO ! J 
+      ENDDO ! ICUT
+! diagnostic output
+!d    do icut=1,ncut
+!d      write(iout,*) 'Cut-off=',rcutoff(icut)
+!d      write(iout,'(8f10.3)') (y(j,icut),j=1,ncur(icut))
+!d    enddo
+!
+! Generate LaTeX script for tree plot
+!
+      iylen=ylen
+#ifdef AIX
+      call fdate_(fd)
+#else
+      call fdate(fd)
+#endif
+      write(jplot,'(80(1h%))')
+      write(jplot,'(a)')  '% LaTeX code for minimal-tree plotting.'
+      write(jplot,'(3a)') '% Created by UNRES_CLUST on ',&
+        fd(:ilen(fd)),'.'
+      write(jplot,'(2a)') '% To change the dimensions use the LaTeX',&
+       ' \\unitlength=number command.'
+      write(jplot,'(a)') '% The default dimensions fit an A4 page.'
+      write(jplot,'(80(1h%))')
+      write(jplot,'(a,i5,a)') '\\begin{picture}(1,1)(0,',iylen,')'
+      ycur=ylen+ybox 
+      do icut=ncut,1,-1
+        xcur=xlen*(icut-1)
+        write(jplot,'(a,f6.1,a,f6.1,a,f4.2,a)') &
+         '  \\put(',xcur,',',ycur,'){',rcutoff(icut),' \\AA}' 
+      enddo ! icut
+      xcur=0.0
+      xdraw=xcur+xbox
+      nnc=ncur(1)
+      write(jplot,'(a,i3,a)') '% Begin cut-off',1,'.'
+      do j=1,nnc
+        ydraw=y(j,1)
+        ycur=ydraw-0.5*ybox
+        ideltx=deltx
+        write(jplot,'(4(a,f6.1),a,i3,a)') &
+         '  \\put(',xcur,',',ycur,'){\\framebox(',xbox,',',ybox,'){',&
+         itree(j,1),'}}'
+        write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)') &
+         '  \\put(',xdraw,',',ydraw,'){\\line(',ideltx,&
+         ',',0,'){',deltx,'}}'
+      enddo ! j
+      do icut=2,ncut
+        write(jplot,'(a,i3,a)') '% Begin cut-off',icut,'.'
+        xcur=xlen*(icut-1)
+        xdraw=xcur-deltx
+!d      print *,'icut=',icut,' xlen=',xlen,' deltx=',deltx,
+!d   & ' xcur=',xcur,' xdraw=',xdraw
+        nnc=ncur(icut)
+        do j=1,ncur(icut-1)
+          ydraw=y(ifirst(j,icut-1),icut)
+          delty=y(ilast(j,icut-1),icut)-y(ifirst(j,icut-1),icut)
+          idelty=delty
+          write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)') &
+         '  \\put(',xdraw,',',ydraw,'){\\line(',0,&
+         ',',idelty,'){',delty,'}}'
+        enddo
+        do j=1,nnc
+          xcur=xlen*(icut-1)
+          xdraw=xcur-deltx
+          ydraw=y(j,icut)
+          ycur=ydraw-0.5*ybox
+          ideltx=deltx
+          write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)') &
+           '  \\put(',xdraw,',',ydraw,'){\\line(',ideltx,&
+           ',',0,'){',deltx,'}}'
+          write(jplot,'(4(a,f6.1),a,i3,a)') &
+           '  \\put(',xcur,',',ycur,'){\\framebox(',xbox,',',ybox,'){',&
+           itree(j,icut),'}}'
+          if (icut.lt.ncut) then
+            xdraw=xcur+xbox
+            write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)') & 
+           '  \\put(',xdraw,',',ydraw,'){\\line(',ideltx,&
+           ',',0,'){',deltx,'}}'
+          endif
+        enddo ! j
+      enddo ! icut
+      write(jplot,'(a)') '\\end{picture}'
+      RETURN
+      END SUBROUTINE PLOTREE
+!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
+      end module tracking