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