external ilen
[unres.git] / source / cluster / wham / src-M / bakup / main_clust.F
1 C
2 C Program to cluster united-residue MCM results.
3 C
4       implicit none
5       include 'DIMENSIONS'
6       include 'sizesclu.dat'
7 #ifdef MPI
8       include "mpif.h"
9       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
10       include "COMMON.MPI"
11 #endif
12       include 'COMMON.TIME1'
13       include 'COMMON.INTERACT'
14       include 'COMMON.NAMES'
15       include 'COMMON.GEO'
16       include 'COMMON.HEADER'
17       include 'COMMON.CONTROL'
18       include 'COMMON.CHAIN'
19       include 'COMMON.VAR'
20       include 'COMMON.CLUSTER'
21       include 'COMMON.IOUNITS'
22       include 'COMMON.FREE'
23       logical printang(max_cut)
24       integer printpdb(max_cut)
25       integer printmol2(max_cut)
26       character*240 lineh
27       REAL CRIT(maxconf),MEMBR(maxconf)
28       REAL CRITVAL(maxconf-1)
29       INTEGER IA(maxconf),IB(maxconf)
30       INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
31       INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
32       integer nn,ndis
33       real*4 DISNN
34       DIMENSION NN(maxconf),DISNN(maxconf)
35       LOGICAL FLAG(maxconf)
36       integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon,
37      & it,ncon_work,ind1
38       double precision t1,t2,tcpu,difconf
39       
40       double precision varia(maxvar)
41       double precision hrtime,mintime,sectime
42       logical eof
43 #ifdef MPI
44       call MPI_Init( IERROR )
45       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
46       call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
47       Master = 0
48       if (ierror.gt.0) then
49         write(iout,*) "SEVERE ERROR - Can't initialize MPI."
50         call mpi_finalize(ierror)
51         stop
52       endif
53       if (nprocs.gt.MaxProcs+1) then
54         write (2,*) "Error - too many processors",
55      &   nprocs,MaxProcs+1
56         write (2,*) "Increase MaxProcs and recompile"
57         call MPI_Finalize(IERROR)
58         stop
59       endif
60 #endif
61
62       call initialize
63       call openunits
64       call parmread
65       call read_control
66       call molread
67 c      if (refstr) call read_ref_structure(*30)
68       do i=1,nres
69         phi(i)=0.0D0
70         theta(i)=0.0D0
71         alph(i)=0.0D0
72         omeg(i)=0.0D0
73       enddo
74        write (iout,*) "symetr", symetr
75       call permut(symetr)
76       print *,'MAIN: nnt=',nnt,' nct=',nct
77
78       DO I=1,NCUT
79         PRINTANG(I)=.FALSE.
80         PRINTPDB(I)=0
81         printmol2(i)=0
82         IF (RCUTOFF(I).LT.0.0) THEN
83           RCUTOFF(I)=ABS(RCUTOFF(I))
84           PRINTANG(I)=.TRUE.
85           PRINTPDB(I)=outpdb
86           printmol2(i)=outmol2
87         ENDIF
88       ENDDO
89       write (iout,*) 'Number of cutoffs:',NCUT
90       write (iout,*) 'Cutoff values:'
91       DO ICUT=1,NCUT
92         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
93      &    printpdb(icut),printmol2(icut)
94       ENDDO
95       DO I=1,NRES-3  
96         MULT(I)=1
97       ENDDO
98       do i=1,maxconf
99         list_conf(i)=i
100       enddo
101       call read_coords(ncon,*20)
102       write (iout,*) 'from read_coords: ncon',ncon
103       
104       write (iout,*) "nT",nT
105       do iT=1,nT
106       write (iout,*) "iT",iT
107 #ifdef MPI
108       call work_partition(.true.,ncon)
109 #endif
110
111       call probabl(iT,ncon_work,ncon,*20)
112
113       if (ncon_work.lt.2) then
114         write (iout,*) "Too few conformations; clustering skipped"
115         exit
116       endif
117 #ifdef MPI
118       ndis=ncon_work*(ncon_work-1)/2
119       call work_partition(.true.,ndis)
120 #endif
121
122       DO I=1,NCON_work
123         ICC(I)=I
124       ENDDO
125       WRITE (iout,'(A80)') TITEL
126       t1=tcpu()
127 C
128 C CALCULATE DISTANCES
129 C
130       call daread_ccoords(1,ncon_work)
131       ind1=0
132       DO I=1,NCON_work-1
133         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
134         do k=1,2*nres
135           do l=1,3
136             c(l,k)=allcart(l,k,i)
137           enddo 
138         enddo
139         do k=1,nres
140           do l=1,3
141             cref(l,k)=c(l,k)
142           enddo
143         enddo
144         DO J=I+1,NCON_work
145           IND=IOFFSET(NCON_work,I,J)
146 #ifdef MPI
147           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
148 #endif
149           ind1=ind1+1
150           DISS(IND1)=DIFCONF(I,J)
151 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
152 #ifdef MPI
153           endif
154 #endif
155         ENDDO
156       ENDDO
157       t2=tcpu()
158       WRITE (iout,'(/a,1pe14.5,a/)') 
159      & 'Time for distance calculation:',T2-T1,' sec.'
160       t1=tcpu()
161       PRINT '(a)','End of distance computation'
162
163 #ifdef MPI
164       call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
165      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
166       if (me.eq.master) then
167 #endif
168       open(80,file='/tmp/distance',form='unformatted')
169       do i=1,ndis
170         write(80) diss(i)
171       enddo
172       if (punch_dist) then
173         do i=1,ncon_work-1
174           do j=i+1,ncon_work
175             IND=IOFFSET(NCON,I,J)
176             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
177      &        energy(j)-energy(i)
178           enddo
179         enddo
180       endif
181 C
182 C Print out the RMS deviation matrix.
183 C
184       if (print_dist) CALL DISTOUT(NCON_work)
185 C
186 C  call hierarchical clustering HC from F. Murtagh
187 C
188       N=NCON_work
189       LEN = (N*(N-1))/2
190       write(iout,*) "-------------------------------------------"
191       write(iout,*) "HIERARCHICAL CLUSTERING using"
192       if (iopt.eq.1) then
193         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
194       elseif (iopt.eq.2) then
195         write(iout,*) "SINGLE LINK METHOD"
196       elseif (iopt.eq.3) then
197         write(iout,*) "COMPLETE LINK METHOD"
198       elseif (iopt.eq.4) then
199         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
200       elseif (iopt.eq.5) then
201         write(iout,*) "MCQUITTY'S METHOD"
202       elseif (iopt.eq.6) then
203         write(iout,*) "MEDIAN (GOWER'S) METHOD"
204       elseif (iopt.eq.7) then
205         write(iout,*) "CENTROID METHOD"
206       else
207         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
208         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
209         stop
210       endif
211       write(iout,*)
212       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
213       write(iout,*) "February 1986"
214       write(iout,*) "References:"
215       write(iout,*) "1. Multidimensional clustering algorithms"
216       write(iout,*) "   Fionn Murtagh"
217       write(iout,*) "   Vienna : Physica-Verlag, 1985."
218       write(iout,*) "2. Multivariate data analysis"
219       write(iout,*) "   Fionn Murtagh and Andre Heck"
220       write(iout,*) "   Kluwer Academic Publishers, 1987"
221       write(iout,*) "-------------------------------------------"
222       write(iout,*)
223
224 #ifdef DEBUG
225       write (iout,*) "The TOTFREE array"
226       do i=1,ncon_work
227         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
228       enddo
229 #endif
230       call flush(iout)
231       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
232       LEV = N-1
233       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
234       if (lev.lt.2) then
235         write (iout,*) "Too few conformations to cluster."
236         goto 192
237       endif
238       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
239 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
240
241       do i=1,maxgr
242         licz(i)=0
243       enddo
244       icut=1
245       i=1
246       NGR=i+1
247       do j=1,n
248         licz(iclass(j,i))=licz(iclass(j,i))+1
249         nconf(iclass(j,i),licz(iclass(j,i)))=j
250 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
251 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
252       enddo        
253       do i=1,lev-1
254
255          idum=lev-i
256          DO L=1,LEV
257             IF (HEIGHT(L).EQ.IDUM) GOTO 190
258          ENDDO
259  190     IDUM=L
260          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
261      &    " icut",icut," cutoff",rcutoff(icut)
262          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
263           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
264           write (iout,'(a,f8.2)') 'Maximum distance found:',
265      &              CRITVAL(IDUM)
266           CALL SRTCLUST(ICUT,ncon_work,iT)
267           CALL TRACK(ICUT)
268           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
269           icut=icut+1
270           if (icut.gt.ncut) goto 191
271          ENDIF
272          NGR=i+1
273          do l=1,maxgr
274           licz(l)=0
275          enddo
276          do j=1,n
277           licz(iclass(j,i))=licz(iclass(j,i))+1
278           nconf(iclass(j,i),licz(iclass(j,i)))=j
279 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
280 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
281 cd          print *,j,iclass(j,i),
282 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
283          enddo
284       enddo
285  191  continue
286 C
287       if (plot_tree) then
288         CALL WRITRACK
289         CALL PLOTREE
290       endif
291 C
292       t2=tcpu()
293       WRITE (iout,'(/a,1pe14.5,a/)') 
294      & 'Total time for clustering:',T2-T1,' sec.'
295
296 #ifdef MPI
297       endif
298 #endif
299  192  continue
300       enddo
301 C
302       close(icbase,status="delete")
303 #ifdef MPI
304       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
305 #endif
306       stop '********** Program terminated normally.'
307    20 write (iout,*) "Error reading coordinates"
308 #ifdef MPI
309       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
310 #endif
311       stop
312    30 write (iout,*) "Error reading reference structure"
313 #ifdef MPI
314       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
315 #endif
316       stop
317       end
318 c---------------------------------------------------------------------------
319       double precision function difconf(icon,jcon)
320       implicit none
321       include 'DIMENSIONS'
322       include 'sizesclu.dat'
323       include 'COMMON.CONTROL'
324       include 'COMMON.CLUSTER'
325       include 'COMMON.CHAIN' 
326       include 'COMMON.INTERACT'
327       include 'COMMON.VAR'
328       include 'COMMON.IOUNITS'
329       logical non_conv
330       double precision przes(3),obrot(3,3)
331       double precision xx(3,maxres2),yy(3,maxres2)
332       integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
333       integer iaperm,ibezperm,run
334       double precision rms,rmsmina
335 c      write (iout,*) "tu dochodze"
336       rmsmina=10d10
337       nperm=1
338       do i=1,symetr
339       nperm=i*nperm
340       enddo
341 c      write (iout,*) "nperm",nperm
342       call permut(symetr)
343 c      write (iout,*) "tabperm", tabperm(1,1)
344       do kkk=1,nperm
345       if (lside) then
346         ii=0
347         chalen=int((nend-nstart+2)/symetr)
348         do run=1,symetr
349          do i=nstart,(nstart+chalen-1)
350           zzz=tabperm(kkk,run)
351 c          write (iout,*) "tutaj",zzz
352           ii=ii+1
353           iaperm=(zzz-1)*chalen+i
354           ibezperm=(run-1)*chalen+i
355           do j=1,3
356             xx(j,ii)=allcart(j,iaperm,jcon)
357             yy(j,ii)=cref(j,ibezperm)
358           enddo
359          enddo
360         enddo
361         do run=1,symetr
362          do i=nstart,(nstart+chalen-1)
363           zzz=tabperm(kkk,run)
364           ii=ii+1
365           iaperm=(zzz-1)*chalen+i
366           ibezperm=(run-1)*chalen+i
367 c          if (itype(i).ne.10) then
368             ii=ii+1
369             do j=1,3 
370               xx(j,ii)=allcart(j,iaperm+nres,jcon)
371               yy(j,ii)=cref(j,ibezperm+nres)
372             enddo
373            enddo
374 c          endif
375         enddo
376         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
377       else
378         chalen=int((nct-nnt+2)/symetr)
379         do run=1,symetr
380          do i=nnt,(nnt+chalen-1)
381           zzz=tabperm(kkk,run)
382 c           write (iout,*) "tu szukaj", zzz,run,kkk
383           iaperm=(zzz-1)*chalen+i
384           ibezperm=(run-1)*chalen+i
385 c        do i=nnt,nct
386           do j=1,3
387             c(j,i)=allcart(j,iaperm,jcon)
388           enddo
389          enddo
390         enddo
391         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
392      &       obrot,non_conv)
393       endif
394       if (rms.lt.0.0) then
395         print *,'error, rms^2 = ',rms,icon,jcon
396         stop
397       endif
398       if (non_conv) print *,non_conv,icon,jcon
399       if (rmsmina.gt.rms) rmsmina=rms
400       enddo
401       difconf=dsqrt(rmsmina)
402       RETURN
403       END
404 C------------------------------------------------------------------------------
405       subroutine distout(ncon)
406       implicit none
407       include 'DIMENSIONS'
408       include 'sizesclu.dat'
409       integer ncol,ncon
410       parameter (ncol=10)
411       include 'COMMON.IOUNITS'
412       include 'COMMON.CLUSTER'
413       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
414       real*4 b
415       dimension b(ncol)
416       write (iout,'(a)') 'The distance matrix'
417       do 1 i=1,ncon,ncol
418       nlim=min0(i+ncol-1,ncon)
419       write (iout,1000) (k,k=i,nlim)
420       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
421  1000 format (/8x,10(i4,3x))
422  1020 format (/1x,80(1h-)/)
423       do 2 j=i,ncon
424       jlim=min0(j,nlim)
425       if (jlim.eq.j) then
426         b(jlim-i+1)=0.0d0
427         jlim1=jlim-1
428       else
429         jlim1=jlim
430       endif
431       do 3 k=i,jlim1
432        if (j.lt.k) then 
433           IND=IOFFSET(NCON,j,k)
434        else
435           IND=IOFFSET(NCON,k,j)
436        endif
437     3  b(k-i+1)=diss(IND)
438       write (iout,1010) j,(b(k),k=1,jlim-i+1)
439     2 continue
440     1 continue
441  1010 format (i5,3x,10(f6.2,1x))
442       return
443       end