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