arcos = 0.5D0*(PI-DSIGN(1.0D0,X)*PI)
[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
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 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
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       ind1=0
136       DO I=1,NCON_work-1
137         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
138         do k=1,2*nres
139           do l=1,3
140             c(l,k)=allcart(l,k,i)
141           enddo 
142         enddo
143         do k=1,nres
144           do l=1,3
145             cref(l,k)=c(l,k)
146           enddo
147         enddo
148         DO J=I+1,NCON_work
149           IND=IOFFSET(NCON_work,I,J)
150 #ifdef MPI
151           if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
152 #endif
153           ind1=ind1+1
154           DISS(IND1)=DIFCONF(I,J)
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       open(80,file='/tmp/distance',form='unformatted')
173       do i=1,ndis
174         write(80) diss(i)
175       enddo
176       if (punch_dist) then
177         do i=1,ncon_work-1
178           do j=i+1,ncon_work
179             IND=IOFFSET(NCON,I,J)
180             write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
181      &        energy(j)-energy(i)
182           enddo
183         enddo
184       endif
185 C
186 C Print out the RMS deviation matrix.
187 C
188       if (print_dist) CALL DISTOUT(NCON_work)
189 C
190 C  call hierarchical clustering HC from F. Murtagh
191 C
192       N=NCON_work
193       LEN = (N*(N-1))/2
194       write(iout,*) "-------------------------------------------"
195       write(iout,*) "HIERARCHICAL CLUSTERING using"
196       if (iopt.eq.1) then
197         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
198       elseif (iopt.eq.2) then
199         write(iout,*) "SINGLE LINK METHOD"
200       elseif (iopt.eq.3) then
201         write(iout,*) "COMPLETE LINK METHOD"
202       elseif (iopt.eq.4) then
203         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
204       elseif (iopt.eq.5) then
205         write(iout,*) "MCQUITTY'S METHOD"
206       elseif (iopt.eq.6) then
207         write(iout,*) "MEDIAN (GOWER'S) METHOD"
208       elseif (iopt.eq.7) then
209         write(iout,*) "CENTROID METHOD"
210       else
211         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
212         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
213         stop
214       endif
215       write(iout,*)
216       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
217       write(iout,*) "February 1986"
218       write(iout,*) "References:"
219       write(iout,*) "1. Multidimensional clustering algorithms"
220       write(iout,*) "   Fionn Murtagh"
221       write(iout,*) "   Vienna : Physica-Verlag, 1985."
222       write(iout,*) "2. Multivariate data analysis"
223       write(iout,*) "   Fionn Murtagh and Andre Heck"
224       write(iout,*) "   Kluwer Academic Publishers, 1987"
225       write(iout,*) "-------------------------------------------"
226       write(iout,*)
227
228 #ifdef DEBUG
229       write (iout,*) "The TOTFREE array"
230       do i=1,ncon_work
231         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
232       enddo
233 #endif
234       call flush(iout)
235       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
236       LEV = N-1
237       write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
238       if (lev.lt.2) then
239         write (iout,*) "Too few conformations to cluster."
240         goto 192
241       endif
242       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
243 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
244
245       do i=1,maxgr
246         licz(i)=0
247       enddo
248       icut=1
249       i=1
250       NGR=i+1
251       do j=1,n
252         licz(iclass(j,i))=licz(iclass(j,i))+1
253         nconf(iclass(j,i),licz(iclass(j,i)))=j
254 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
255 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
256       enddo        
257       do i=1,lev-1
258
259          idum=lev-i
260          DO L=1,LEV
261             IF (HEIGHT(L).EQ.IDUM) GOTO 190
262          ENDDO
263  190     IDUM=L
264          write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
265      &    " icut",icut," cutoff",rcutoff(icut)
266          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
267           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
268           write (iout,'(a,f8.2)') 'Maximum distance found:',
269      &              CRITVAL(IDUM)
270           CALL SRTCLUST(ICUT,ncon_work,iT)
271           CALL TRACK(ICUT)
272           CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
273           icut=icut+1
274           if (icut.gt.ncut) goto 191
275          ENDIF
276          NGR=i+1
277          do l=1,maxgr
278           licz(l)=0
279          enddo
280          do j=1,n
281           licz(iclass(j,i))=licz(iclass(j,i))+1
282           nconf(iclass(j,i),licz(iclass(j,i)))=j
283 c        write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
284 c     &    nconf(iclass(j,i),licz(iclass(j,i)))
285 cd          print *,j,iclass(j,i),
286 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
287          enddo
288       enddo
289  191  continue
290 C
291       if (plot_tree) then
292         CALL WRITRACK
293         CALL PLOTREE
294       endif
295 C
296       t2=tcpu()
297       WRITE (iout,'(/a,1pe14.5,a/)') 
298      & 'Total time for clustering:',T2-T1,' sec.'
299
300 #ifdef MPI
301       endif
302 #endif
303  192  continue
304       enddo
305 C
306       close(icbase,status="delete")
307 #ifdef MPI
308       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
309 #endif
310       stop '********** Program terminated normally.'
311    20 write (iout,*) "Error reading coordinates"
312 #ifdef MPI
313       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
314 #endif
315       stop
316    30 write (iout,*) "Error reading reference structure"
317 #ifdef MPI
318       call MPI_Finalize(MPI_COMM_WORLD,IERROR)
319 #endif
320       stop
321       end
322 c---------------------------------------------------------------------------
323       double precision function difconf(icon,jcon)
324       implicit none
325       include 'DIMENSIONS'
326       include 'sizesclu.dat'
327       include 'COMMON.CONTROL'
328       include 'COMMON.CLUSTER'
329       include 'COMMON.CHAIN' 
330       include 'COMMON.INTERACT'
331       include 'COMMON.VAR'
332       include 'COMMON.IOUNITS'
333       logical non_conv
334       double precision przes(3),obrot(3,3)
335       double precision xx(3,maxres2),yy(3,maxres2)
336       integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
337       integer iaperm,ibezperm,run
338       double precision rms,rmsmina
339 c      write (iout,*) "tu dochodze"
340       rmsmina=10d10
341       nperm=1
342       do i=1,symetr
343       nperm=i*nperm
344       enddo
345 c      write (iout,*) "nperm",nperm
346       call permut(symetr)
347 c      write (iout,*) "tabperm", tabperm(1,1)
348       do kkk=1,nperm
349       if (lside) then
350         ii=0
351         chalen=int((nend-nstart+2)/symetr)
352         do run=1,symetr
353          do i=nstart,(nstart+chalen-1)
354           zzz=tabperm(kkk,run)
355 c          write (iout,*) "tutaj",zzz
356           ii=ii+1
357           iaperm=(zzz-1)*chalen+i
358           ibezperm=(run-1)*chalen+i
359           do j=1,3
360             xx(j,ii)=allcart(j,iaperm,jcon)
361             yy(j,ii)=cref(j,ibezperm)
362           enddo
363          enddo
364         enddo
365         do run=1,symetr
366          do i=nstart,(nstart+chalen-1)
367           zzz=tabperm(kkk,run)
368           ii=ii+1
369           iaperm=(zzz-1)*chalen+i
370           ibezperm=(run-1)*chalen+i
371 c          if (itype(i).ne.10) then
372             ii=ii+1
373             do j=1,3 
374               xx(j,ii)=allcart(j,iaperm+nres,jcon)
375               yy(j,ii)=cref(j,ibezperm+nres)
376             enddo
377            enddo
378 c          endif
379         enddo
380         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
381       else
382         chalen=int((nct-nnt+2)/symetr)
383         do run=1,symetr
384          do i=nnt,(nnt+chalen-1)
385           zzz=tabperm(kkk,run)
386 c           write (iout,*) "tu szukaj", zzz,run,kkk
387           iaperm=(zzz-1)*chalen+i
388           ibezperm=(run-1)*chalen+i
389 c        do i=nnt,nct
390           do j=1,3
391             c(j,i)=allcart(j,iaperm,jcon)
392           enddo
393          enddo
394         enddo
395         call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
396      &       obrot,non_conv)
397       endif
398       if (rms.lt.0.0) then
399         print *,'error, rms^2 = ',rms,icon,jcon
400         stop
401       endif
402       if (non_conv) print *,non_conv,icon,jcon
403       if (rmsmina.gt.rms) rmsmina=rms
404       enddo
405       difconf=dsqrt(rmsmina)
406       RETURN
407       END
408 C------------------------------------------------------------------------------
409       subroutine distout(ncon)
410       implicit none
411       include 'DIMENSIONS'
412       include 'sizesclu.dat'
413       integer ncol,ncon
414       parameter (ncol=10)
415       include 'COMMON.IOUNITS'
416       include 'COMMON.CLUSTER'
417       integer i,j,k,jlim,jlim1,nlim,ind,ioffset
418       real*4 b
419       dimension b(ncol)
420       write (iout,'(a)') 'The distance matrix'
421       do 1 i=1,ncon,ncol
422       nlim=min0(i+ncol-1,ncon)
423       write (iout,1000) (k,k=i,nlim)
424       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
425  1000 format (/8x,10(i4,3x))
426  1020 format (/1x,80(1h-)/)
427       do 2 j=i,ncon
428       jlim=min0(j,nlim)
429       if (jlim.eq.j) then
430         b(jlim-i+1)=0.0d0
431         jlim1=jlim-1
432       else
433         jlim1=jlim
434       endif
435       do 3 k=i,jlim1
436        if (j.lt.k) then 
437           IND=IOFFSET(NCON,j,k)
438        else
439           IND=IOFFSET(NCON,k,j)
440        endif
441     3  b(k-i+1)=diss(IND)
442       write (iout,1010) j,(b(k),k=1,jlim-i+1)
443     2 continue
444     1 continue
445  1010 format (i5,3x,10(f6.2,1x))
446       return
447       end