max_template instead of fix number 19 in wham and cluster
[unres.git] / source / cluster / wham / src / 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
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,ilen
38       double precision t1,t2,tcpu,difconf
39       real diss_(maxdist)
40       
41       double precision varia(maxvar)
42       double precision hrtime,mintime,sectime
43       logical eof
44       external ilen
45 #ifdef MPI
46       call MPI_Init( IERROR )
47       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
48       call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
49       Master = 0
50       if (ierror.gt.0) then
51         write(iout,*) "SEVERE ERROR - Can't initialize MPI."
52         call mpi_finalize(ierror)
53         stop
54       endif
55       if (nprocs.gt.MaxProcs+1) then
56         write (2,*) "Error - too many processors",
57      &   nprocs,MaxProcs+1
58         write (2,*) "Increase MaxProcs and recompile"
59         call MPI_Finalize(IERROR)
60         stop
61       endif
62 #endif
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
77       print *,'MAIN: nnt=',nnt,' nct=',nct
78
79       DO I=1,NCUT
80         PRINTANG(I)=.FALSE.
81         PRINTPDB(I)=0
82         printmol2(i)=0
83         IF (RCUTOFF(I).LT.0.0) THEN
84           RCUTOFF(I)=ABS(RCUTOFF(I))
85           PRINTANG(I)=.TRUE.
86           PRINTPDB(I)=outpdb
87           printmol2(i)=outmol2
88         ENDIF
89       ENDDO
90       write (iout,*) 'Number of cutoffs:',NCUT
91       write (iout,*) 'Cutoff values:'
92       DO ICUT=1,NCUT
93         WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
94      &    printpdb(icut),printmol2(icut)
95       ENDDO
96       DO I=1,NRES-3  
97         MULT(I)=1
98       ENDDO
99       do i=1,maxconf
100         list_conf(i)=i
101       enddo
102       call read_coords(ncon,*20)
103       write (iout,*) 'from read_coords: ncon',ncon
104       
105       write (iout,*) "nT",nT
106       do iT=1,nT
107       write (iout,*) "iT",iT
108 #ifdef MPI
109       call work_partition(.true.,ncon)
110 #endif
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 #ifdef MPI
151           DISS_(IND1)=DIFCONF(I,J)
152 #else
153           DISS(IND1)=DIFCONF(I,J)
154 #endif
155 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
156 #ifdef MPI
157           endif
158 #endif
159         ENDDO
160       ENDDO
161       t2=tcpu()
162       WRITE (iout,'(/a,1pe14.5,a/)') 
163      & 'Time for distance calculation:',T2-T1,' sec.'
164       t1=tcpu()
165       PRINT '(a)','End of distance computation'
166
167 #ifdef MPI
168       call MPI_Gatherv(diss_(1),scount(me),MPI_REAL,diss(1),
169      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
170       if (me.eq.master) then
171 #endif
172       scrachdir2d=scratchdir(:ilen(scratchdir))//'distance'
173       open(80,file=scrachdir2d,form='unformatted')
174       do i=1,ndis
175         write(80) diss(i)
176       enddo
177       if (punch_dist) then
178         do i=1,ncon_work-1
179           do j=i+1,ncon_work
180             IND=IOFFSET(NCON,I,J)
181             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
182      &        energy(j)-energy(i)
183           enddo
184         enddo
185       endif
186 C
187 C Print out the RMS deviation matrix.
188 C
189       if (print_dist) CALL DISTOUT(NCON_work)
190 C
191 C  call hierarchical clustering HC from F. Murtagh
192 C
193       N=NCON_work
194       LEN = (N*(N-1))/2
195       write(iout,*) "-------------------------------------------"
196       write(iout,*) "HIERARCHICAL CLUSTERING using"
197       if (iopt.eq.1) then
198         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
199       elseif (iopt.eq.2) then
200         write(iout,*) "SINGLE LINK METHOD"
201       elseif (iopt.eq.3) then
202         write(iout,*) "COMPLETE LINK METHOD"
203       elseif (iopt.eq.4) then
204         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
205       elseif (iopt.eq.5) then
206         write(iout,*) "MCQUITTY'S METHOD"
207       elseif (iopt.eq.6) then
208         write(iout,*) "MEDIAN (GOWER'S) METHOD"
209       elseif (iopt.eq.7) then
210         write(iout,*) "CENTROID METHOD"
211       else
212         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
213         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
214         stop
215       endif
216       write(iout,*)
217       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
218       write(iout,*) "February 1986"
219       write(iout,*) "References:"
220       write(iout,*) "1. Multidimensional clustering algorithms"
221       write(iout,*) "   Fionn Murtagh"
222       write(iout,*) "   Vienna : Physica-Verlag, 1985."
223       write(iout,*) "2. Multivariate data analysis"
224       write(iout,*) "   Fionn Murtagh and Andre Heck"
225       write(iout,*) "   Kluwer Academic Publishers, 1987"
226       write(iout,*) "-------------------------------------------"
227       write(iout,*)
228
229 #ifdef DEBUG
230       write (iout,*) "The TOTFREE array"
231       do i=1,ncon_work
232         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
233       enddo
234 #endif
235 #ifdef AIX
236       call flush_(iout)
237 #else
238       call flush(iout)
239 #endif
240       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
241       LEV = N-1
242       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
243       if (lev.lt.2) then
244         write (iout,*) "Too few conformations to cluster."
245         goto 192
246       endif
247       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
248 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
249
250       do i=1,maxgr
251         licz(i)=0
252       enddo
253       icut=1
254       i=1
255       NGR=i+1
256       do j=1,n
257         licz(iclass(j,i))=licz(iclass(j,i))+1
258         nconf(iclass(j,i),licz(iclass(j,i)))=j
259 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
260 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
261       enddo        
262       do i=1,lev-1
263
264          idum=lev-i
265          DO L=1,LEV
266             IF (HEIGHT(L).EQ.IDUM) GOTO 190
267          ENDDO
268  190     IDUM=L
269          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
270      &    " icut",icut," cutoff",rcutoff(icut)
271          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
272           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
273           write (iout,'(a,f8.2)') 'Maximum distance found:',
274      &              CRITVAL(IDUM)
275           CALL SRTCLUST(ICUT,ncon_work,iT)
276           CALL TRACK(ICUT)
277           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
278           icut=icut+1
279           if (icut.gt.ncut) goto 191
280          ENDIF
281          NGR=i+1
282          do l=1,maxgr
283           licz(l)=0
284          enddo
285          do j=1,n
286           licz(iclass(j,i))=licz(iclass(j,i))+1
287           nconf(iclass(j,i),licz(iclass(j,i)))=j
288 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
289 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
290 cd          print *,j,iclass(j,i),
291 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
292          enddo
293       enddo
294  191  continue
295 C
296       if (plot_tree) then
297         CALL WRITRACK
298         CALL PLOTREE
299       endif
300 C
301       t2=tcpu()
302       WRITE (iout,'(/a,1pe14.5,a/)') 
303      & 'Total time for clustering:',T2-T1,' sec.'
304
305 #ifdef MPI
306       endif
307 #endif
308  192  continue
309       enddo
310 C
311       close(icbase,status="delete")
312 #ifdef MPI
313       call MPI_Finalize(IERROR)
314 #endif
315       stop '********** Program terminated normally.'
316    20 write (iout,*) "Error reading coordinates"
317 #ifdef MPI
318       call MPI_Finalize(IERROR)
319 #endif
320       stop
321    30 write (iout,*) "Error reading reference structure"
322 #ifdef MPI
323       call MPI_Finalize(IERROR)
324 #endif
325       stop
326       end
327 c---------------------------------------------------------------------------
328       double precision function difconf(icon,jcon)
329       implicit none
330       include 'DIMENSIONS'
331       include 'sizesclu.dat'
332       include 'COMMON.CONTROL'
333       include 'COMMON.CLUSTER'
334       include 'COMMON.CHAIN' 
335       include 'COMMON.INTERACT'
336       include 'COMMON.VAR'
337       include 'COMMON.IOUNITS'
338       logical non_conv
339       double precision przes(3),obrot(3,3)
340       double precision xx(3,maxres2),yy(3,maxres2)
341       integer i,ii,j,icon,jcon
342       double precision rms
343       if (lside) then
344         ii=0
345         do i=nstart,nend
346           ii=ii+1
347           do j=1,3
348             xx(j,ii)=allcart(j,i,jcon)
349             yy(j,ii)=cref(j,i)
350           enddo
351         enddo
352         do i=nstart,nend
353 c          if (itype(i).ne.10) then
354             ii=ii+1
355             do j=1,3 
356               xx(j,ii)=allcart(j,i+nres,jcon)
357               yy(j,ii)=cref(j,i+nres)
358             enddo
359 c          endif
360         enddo
361         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
362       else
363         do i=nnt,nct
364           do j=1,3
365             c(j,i)=allcart(j,i,jcon)
366           enddo
367         enddo
368         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
369      &       obrot,non_conv)
370       endif
371       if (rms.lt.0.0) then
372         print *,'error, rms^2 = ',rms,icon,jcon
373         stop
374       endif
375       if (non_conv) print *,non_conv,icon,jcon
376       difconf=dsqrt(rms)
377       RETURN
378       END
379 C------------------------------------------------------------------------------
380       subroutine distout(ncon)
381       implicit none
382       include 'DIMENSIONS'
383       include 'sizesclu.dat'
384       integer ncol,ncon
385       parameter (ncol=10)
386       include 'COMMON.IOUNITS'
387       include 'COMMON.CLUSTER'
388       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
389       real*4 b
390       dimension b(ncol)
391       write (iout,'(a)') 'The distance matrix'
392       do 1 i=1,ncon,ncol
393       nlim=min0(i+ncol-1,ncon)
394       write (iout,1000) (k,k=i,nlim)
395       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
396  1000 format (/8x,10(i4,3x))
397  1020 format (/1x,80(1h-)/)
398       do 2 j=i,ncon
399       jlim=min0(j,nlim)
400       if (jlim.eq.j) then
401         b(jlim-i+1)=0.0d0
402         jlim1=jlim-1
403       else
404         jlim1=jlim
405       endif
406       do 3 k=i,jlim1
407        if (j.lt.k) then 
408           IND=IOFFSET(NCON,j,k)
409        else
410           IND=IOFFSET(NCON,k,j)
411        endif
412     3  b(k-i+1)=diss(IND)
413       write (iout,1010) j,(b(k),k=1,jlim-i+1)
414     2 continue
415     1 continue
416  1010 format (i5,3x,10(f6.2,1x))
417       return
418       end