Lorentzian restraints on distances in cluster_wham
[unres.git] / source / cluster / wham / src / track.F
1       SUBROUTINE TRACK(ICUT)
2       include 'DIMENSIONS'
3       INCLUDE 'sizesclu.dat'
4       INCLUDE 'COMMON.CLUSTER'
5       COMMON /HISTORY/ NCUR(MAX_CUT),IBACK(MAXGR,MAX_CUT)
6       COMMON /PREVIOUS/ NGRP,LICZP(MAXGR),NCONFP(MAXGR,MAXINGR)
7       IF (ICUT.GT.1) THEN
8 C Find out what of the previous families the current ones came from.        
9         DO IGR=1,NGR
10           NCI1=NCONF(IGR,1) 
11           DO JGR=1,NGRP
12             DO K=1,LICZP(JGR)
13               IF (NCI1.EQ.NCONFP(JGR,K)) THEN
14                 IBACK(IGR,ICUT)=JGR
15                 GOTO 10
16               ENDIF
17             ENDDO ! K
18           ENDDO ! JGR
19   10      CONTINUE
20         ENDDO ! IGR
21       ENDIF ! (ICUT.GT.1)
22 C Save current partition for subsequent backtracking.
23       NCUR(ICUT)=NGR
24       NGRP=NGR
25       DO IGR=1,NGR
26         LICZP(IGR)=LICZ(IGR)
27         DO K=1,LICZ(IGR)
28           NCONFP(IGR,K)=NCONF(IGR,K)
29         ENDDO ! K
30       ENDDO ! IGR
31       RETURN
32       END
33 C------------------------------------------------------------------------------
34       SUBROUTINE WRITRACK
35       include 'DIMENSIONS'
36       INCLUDE 'sizesclu.dat'
37       INCLUDE 'COMMON.CLUSTER'
38       include 'COMMON.IOUNITS'
39       COMMON /HISTORY/ NCUR(MAX_CUT),IBACK(MAXGR,MAX_CUT)
40       COMMON /PREVIOUS/ NGRP,LICZP(MAXGR),NCONFP(MAXGR,MAXINGR)
41       DIMENSION IPART(MAXGR/5,MAXGR/5)
42 c     do icut=2,ncut
43 c       write (iout,'(a,f10.5)') 'Cut-off',rcutoff(icut)
44 c       write (iout,'(16i5)') (iback(k,icut),k=1,ncur(icut))
45 c     enddo
46 C
47 C Print the partition history.
48 C
49       DO ICUT=2,NCUT
50         NCU=NCUR(ICUT)
51         NCUP=NCUR(ICUT-1)
52 cd      print *,'icut=',icut,' ncu=',ncu,' ncur=',ncur
53         WRITE(iout,'(A,f10.5,A,f10.5)') 
54      &  'Partition of families obtained at cut-off',RCUTOFF(ICUT-1),
55      &  ' at cut-off',RCUTOFF(ICUT)
56         DO I=1,NCUP
57           NPART=0
58 cd        print *,'i=',i
59           DO J=1,NCU
60             IF (IBACK(J,ICUT).EQ.I) THEN
61               NPART=NPART+1  
62               IPART(NPART,I)=J
63             ENDIF
64 cd          print *,'j=',j,' iback=',IBACK(J,ICUT),' npart=',npart
65           ENDDO ! J
66           WRITE (iout,'(16I5)') I,(IPART(K,I),K=1,NPART) 
67         ENDDO ! I
68       ENDDO ! ICUT
69       RETURN
70       END
71 C------------------------------------------------------------------------------
72       SUBROUTINE PLOTREE
73       include 'DIMENSIONS'
74       INCLUDE 'sizesclu.dat'
75       INCLUDE 'COMMON.CLUSTER'
76       include 'COMMON.IOUNITS'
77       COMMON /HISTORY/ NCUR(MAX_CUT),IBACK(MAXGR,MAX_CUT)
78       COMMON /PREVIOUS/ NGRP,LICZP(MAXGR),NCONFP(MAXGR,MAXINGR)
79       DIMENSION Y(MAXGR,MAX_CUT)
80       DIMENSION ITREE(MAXGR,MAX_CUT),IFIRST(MAXGR,MAX_CUT),
81      &ILAST(MAXGR,MAX_CUT),IFT(MAXGR),ILT(MAXGR),ITR(MAXGR)
82       CHARACTER*32 FD
83       external ilen
84
85 C Generate the image of the tree (tentatively for LaTeX picture environment).
86 C
87 C
88 C First untangle the branches of the tree
89 C
90       DO I=1,NCUR(1)
91         ITREE(I,1)=I
92       ENDDO
93       DO ICUT=NCUT,2,-1
94 C
95 C Determine the order of families for the (icut)th partition.
96 C
97         NCU=NCUR(ICUT)
98         NCUP=NCUR(ICUT-1)
99         NPART=0
100         DO I=1,NCUP
101           IS=0
102           IF (I.GT.1) ILAST(I-1,ICUT-1)=NPART
103           DO J=1,NCU
104             IF (IBACK(J,ICUT).EQ.I) THEN
105               NPART=NPART+1  
106               IF (IS.EQ.0) THEN
107                 IS=1
108                 IFIRST(I,ICUT-1)=NPART
109               ENDIF
110               ITREE(NPART,ICUT)=J
111             ENDIF
112           ENDDO ! J
113         ENDDO ! I
114         ILAST(NCUP,ICUT-1)=NPART
115 cd      print *,'i=',i,' ncup=',ncup,' ncu=',ncu,' npart=',npart
116       ENDDO ! ICUT
117 c diagnostic printout
118 cd    do icut=1,ncut
119 cd      write (iout,*) 'Cut-off',icut,' = ',rcutoff(icut) 
120 cd      write (iout,*) 'ITREE'
121 cd      write (iout,*) (itree(i,icut),i=1,ncur(icut))
122 cd      write (iout,*) 'IFIRST, ILAST'
123 cd      write (iout,*) (ifirst(i,icut),ilast(i,icut),i=1,ncur(icut))
124 cd    enddo
125 C
126 C Propagate the order of families from cut-off #2 to cut-off #n.
127 C
128       DO ICUT=1,NCUT-1
129         DO J=1,NCUR(ICUT)
130           IFT(J)=IFIRST(J,ICUT)
131           ILT(J)=ILAST(J,ICUT)
132         ENDDO ! J
133         DO J=1,NCUR(ICUT+1)
134           ITR(J)=ITREE(J,ICUT+1)
135         ENDDO
136         DO I=1,NCUR(ICUT)
137           ITI=ITREE(I,ICUT)
138 c         write (iout,*) 'icut=',icut,' i=',i,' iti=',iti
139 C         IF (ITI.NE.I) THEN
140             JF1=IFT(I)
141             JF2=IFT(ITI)
142             JL1=ILT(I)
143             JL2=ILT(ITI)
144             JR1=JL1-JF1+1
145             JR2=JL2-JF2+1
146 Cd          write (iout,*) 'jf1=',jf1,' jl1=',jl1,' jf2=',jf2,
147 Cd   &                     ' jl2=',jl2
148 Cd          write (iout,*) 'jr1=',jr1,' jr2=',jr2
149 C Update IFIRST and ILAST.
150             ILAST(I,ICUT)=IFIRST(I,ICUT)+JR2-1
151             IFIRST(I+1,ICUT)=ILAST(I,ICUT)+1
152 C Update ITREE.
153             JF11=IFIRST(I,ICUT)
154 Cd          write(iout,*) 'jf11=',jf11
155             DO J=JF2,JL2
156 Cd            write (iout,*) j,JF11+J-JF2,ITR(J)
157               ITREE(JF11+J-JF2,ICUT+1)=ITR(J)
158             ENDDO
159 Cd      write (iout,*) (ifirst(k,icut),ilast(k,icut),k=1,i)
160 Cd      write (iout,*) (itree(k,icut+1),k=1,ilast(i,icut))
161 C         ENDIF ! (ITI.NE.I)
162         ENDDO ! I
163       ENDDO ! ICUT
164 c diagnostic printout
165 cd    do icut=1,ncut
166 cd      write (iout,*) 'Cut-off',icut,' = ',rcutoff(icut) 
167 cd      write (iout,*) 'ITREE'
168 cd      write (iout,*) (itree(i,icut),i=1,ncur(icut))
169 cd      write (iout,*) 'IFIRST, ILAST'
170 cd      write (iout,*) (ifirst(i,icut),ilast(i,icut),i=1,ncur(icut))
171 cd    enddo
172 C
173 C Generate the y-coordinates of the branches.
174 C
175       XLEN=400.0/(ncut-1)
176       YLEN=600.0
177       xbox=xlen/4.0
178       deltx=0.5*(xlen-xbox)
179       NNC=NCUR(NCUT)
180       ybox=ylen/(2.0*nnc) 
181       DO J=1,NNC
182         Y(J,NCUT)=J*YLEN/NNC
183       ENDDO
184       DO ICUT=NCUT-1,1,-1
185         NNC=NCUR(ICUT)
186         DO J=1,NNC
187           KF=IFIRST(J,ICUT)
188           KL=ILAST(J,ICUT)
189           YY=0.0
190           DO K=KF,KL
191             YY=YY+Y(K,ICUT+1)
192           ENDDO
193           Y(J,ICUT)=YY/(KL-KF+1)
194         ENDDO ! J 
195       ENDDO ! ICUT
196 c diagnostic output
197 cd    do icut=1,ncut
198 cd      write(iout,*) 'Cut-off=',rcutoff(icut)
199 cd      write(iout,'(8f10.3)') (y(j,icut),j=1,ncur(icut))
200 cd    enddo
201 C
202 C Generate LaTeX script for tree plot
203 C
204       iylen=ylen
205 #ifdef AIX
206       call fdate_(fd)
207 #else
208       call fdate(fd)
209 #endif
210       write(jplot,'(80(1h%))')
211       write(jplot,'(a)')  '% LaTeX code for minimal-tree plotting.'
212       write(jplot,'(3a)') '% Created by UNRES_CLUST on ',
213      &  fd(:ilen(fd)),'.'
214       write(jplot,'(2a)') '% To change the dimensions use the LaTeX',
215      & ' \\unitlength=number command.'
216       write(jplot,'(a)') '% The default dimensions fit an A4 page.'
217       write(jplot,'(80(1h%))')
218       write(jplot,'(a,i5,a)') '\\begin{picture}(1,1)(0,',iylen,')'
219       ycur=ylen+ybox 
220       do icut=ncut,1,-1
221         xcur=xlen*(icut-1)
222         write(jplot,'(a,f6.1,a,f6.1,a,f4.2,a)')
223      &   '  \\put(',xcur,',',ycur,'){',rcutoff(icut),' \\AA}' 
224       enddo ! icut
225       xcur=0.0
226       xdraw=xcur+xbox
227       nnc=ncur(1)
228       write(jplot,'(a,i3,a)') '% Begin cut-off',1,'.'
229       do j=1,nnc
230         ydraw=y(j,1)
231         ycur=ydraw-0.5*ybox
232         ideltx=deltx
233         write(jplot,'(4(a,f6.1),a,i3,a)') 
234      &   '  \\put(',xcur,',',ycur,'){\\framebox(',xbox,',',ybox,'){',
235      &   itree(j,1),'}}'
236         write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)') 
237      &   '  \\put(',xdraw,',',ydraw,'){\\line(',ideltx,
238      &   ',',0,'){',deltx,'}}'
239       enddo ! j
240       do icut=2,ncut
241         write(jplot,'(a,i3,a)') '% Begin cut-off',icut,'.'
242         xcur=xlen*(icut-1)
243         xdraw=xcur-deltx
244 cd      print *,'icut=',icut,' xlen=',xlen,' deltx=',deltx,
245 cd   & ' xcur=',xcur,' xdraw=',xdraw
246         nnc=ncur(icut)
247         do j=1,ncur(icut-1)
248           ydraw=y(ifirst(j,icut-1),icut)
249           delty=y(ilast(j,icut-1),icut)-y(ifirst(j,icut-1),icut)
250           idelty=delty
251           write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)')
252      &   '  \\put(',xdraw,',',ydraw,'){\\line(',0,
253      &   ',',idelty,'){',delty,'}}'
254         enddo
255         do j=1,nnc
256           xcur=xlen*(icut-1)
257           xdraw=xcur-deltx
258           ydraw=y(j,icut)
259           ycur=ydraw-0.5*ybox
260           ideltx=deltx
261           write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)') 
262      &     '  \\put(',xdraw,',',ydraw,'){\\line(',ideltx,
263      &     ',',0,'){',deltx,'}}'
264           write(jplot,'(4(a,f6.1),a,i3,a)') 
265      &     '  \\put(',xcur,',',ycur,'){\\framebox(',xbox,',',ybox,'){',
266      &     itree(j,icut),'}}'
267           if (icut.lt.ncut) then
268             xdraw=xcur+xbox
269             write(jplot,'(2(a,f6.1),2(a,i5),a,f6.1,a)') 
270      &     '  \\put(',xdraw,',',ydraw,'){\\line(',ideltx,
271      &     ',',0,'){',deltx,'}}'
272           endif
273         enddo ! j
274       enddo ! icut
275       write(jplot,'(a)') '\\end{picture}'
276       RETURN
277       END