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