Added an option to compile unres, wham, and cluster on the homology branch
[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       
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       call initialize
63       call openunits
64       call cinfo
65       call parmread
66       call read_control
67       call molread
68 c      if (refstr) call read_ref_structure(*30)
69       do i=1,nres
70         phi(i)=0.0D0
71         theta(i)=0.0D0
72         alph(i)=0.0D0
73         omeg(i)=0.0D0
74       enddo
75
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       call probabl(iT,ncon_work,ncon,*20)
111
112       if (ncon_work.lt.2) then
113         write (iout,*) "Too few conformations; clustering skipped"
114         exit
115       endif
116 #ifdef MPI
117       ndis=ncon_work*(ncon_work-1)/2
118       call work_partition(.true.,ndis)
119 #endif
120
121       DO I=1,NCON_work
122         ICC(I)=I
123       ENDDO
124       WRITE (iout,'(A80)') TITEL
125       t1=tcpu()
126 C
127 C CALCULATE DISTANCES
128 C
129       call daread_ccoords(1,ncon_work)
130       ind1=0
131       DO I=1,NCON_work-1
132         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
133         do k=1,2*nres
134           do l=1,3
135             c(l,k)=allcart(l,k,i)
136           enddo 
137         enddo
138         do k=1,nres
139           do l=1,3
140             cref(l,k)=c(l,k)
141           enddo
142         enddo
143         DO J=I+1,NCON_work
144           IND=IOFFSET(NCON_work,I,J)
145 #ifdef MPI
146           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
147 #endif
148           ind1=ind1+1
149           DISS(IND1)=DIFCONF(I,J)
150 c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
151 #ifdef MPI
152           endif
153 #endif
154         ENDDO
155       ENDDO
156       t2=tcpu()
157       WRITE (iout,'(/a,1pe14.5,a/)') 
158      & 'Time for distance calculation:',T2-T1,' sec.'
159       t1=tcpu()
160       PRINT '(a)','End of distance computation'
161
162 #ifdef MPI
163       call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
164      &     scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
165       if (me.eq.master) then
166 #endif
167       scrachdir2d=scratchdir(:ilen(scratchdir))//'distance'
168       open(80,file=scrachdir2d,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
333       double precision rms
334       if (lside) then
335         ii=0
336         do i=nstart,nend
337           ii=ii+1
338           do j=1,3
339             xx(j,ii)=allcart(j,i,jcon)
340             yy(j,ii)=cref(j,i)
341           enddo
342         enddo
343         do i=nstart,nend
344 c          if (itype(i).ne.10) then
345             ii=ii+1
346             do j=1,3 
347               xx(j,ii)=allcart(j,i+nres,jcon)
348               yy(j,ii)=cref(j,i+nres)
349             enddo
350 c          endif
351         enddo
352         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
353       else
354         do i=nnt,nct
355           do j=1,3
356             c(j,i)=allcart(j,i,jcon)
357           enddo
358         enddo
359         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
360      &       obrot,non_conv)
361       endif
362       if (rms.lt.0.0) then
363         print *,'error, rms^2 = ',rms,icon,jcon
364         stop
365       endif
366       if (non_conv) print *,non_conv,icon,jcon
367       difconf=dsqrt(rms)
368       RETURN
369       END
370 C------------------------------------------------------------------------------
371       subroutine distout(ncon)
372       implicit none
373       include 'DIMENSIONS'
374       include 'sizesclu.dat'
375       integer ncol,ncon
376       parameter (ncol=10)
377       include 'COMMON.IOUNITS'
378       include 'COMMON.CLUSTER'
379       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
380       real*4 b
381       dimension b(ncol)
382       write (iout,'(a)') 'The distance matrix'
383       do 1 i=1,ncon,ncol
384       nlim=min0(i+ncol-1,ncon)
385       write (iout,1000) (k,k=i,nlim)
386       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
387  1000 format (/8x,10(i4,3x))
388  1020 format (/1x,80(1h-)/)
389       do 2 j=i,ncon
390       jlim=min0(j,nlim)
391       if (jlim.eq.j) then
392         b(jlim-i+1)=0.0d0
393         jlim1=jlim-1
394       else
395         jlim1=jlim
396       endif
397       do 3 k=i,jlim1
398        if (j.lt.k) then 
399           IND=IOFFSET(NCON,j,k)
400        else
401           IND=IOFFSET(NCON,k,j)
402        endif
403     3  b(k-i+1)=diss(IND)
404       write (iout,1010) j,(b(k),k=1,jlim-i+1)
405     2 continue
406     1 continue
407  1010 format (i5,3x,10(f6.2,1x))
408       return
409       end