dodanie czytania wagi wsccor do clustrowania-M
[unres.git] / source / cluster / unres / src / main_clust.f
1 C
2 C Program to cluster united-residue MCM results.
3 C
4       include 'DIMENSIONS'
5       include 'sizesclu.dat'
6       include 'COMMON.TIME1'
7       include 'COMMON.INTERACT'
8       include 'COMMON.NAMES'
9       include 'COMMON.GEO'
10       include 'COMMON.HEADER'
11       include 'COMMON.CONTROL'
12       include 'COMMON.CONTACTS'
13       include 'COMMON.CHAIN'
14       include 'COMMON.VAR'
15       include 'COMMON.CLUSTER'
16       include 'COMMON.IOUNITS'
17       logical printang(max_cut)
18       integer printpdb(max_cut)
19       integer printmol2(max_cut)
20       character*240 lineh
21       REAL CRIT(maxconf),MEMBR(maxconf)
22       REAL CRITVAL(maxconf-1)
23       REAL DISS(maxdist)
24       INTEGER IA(maxconf),IB(maxconf)
25       INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
26       INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
27       DIMENSION NN(maxconf),DISNN(maxconf)
28       LOGICAL FLAG(maxconf)
29       
30       double precision varia(maxvar)
31       double precision hrtime,mintime,sectime
32       logical eof
33       double precision eee(n_ene,maxconf)
34
35       call initialize
36       call openunits
37       call read_control
38       call molread
39       do i=1,nres
40         phi(i)=0.0D0
41         theta(i)=0.0D0
42         alph(i)=0.0D0
43         omeg(i)=0.0D0
44       enddo
45
46       print *,'MAIN: nnt=',nnt,' nct=',nct
47
48       DO I=1,NCUT
49         PRINTANG(I)=.FALSE.
50         PRINTPDB(I)=0
51         printmol2(i)=0
52         IF (RCUTOFF(I).LT.0.0) THEN
53           RCUTOFF(I)=ABS(RCUTOFF(I))
54           PRINTANG(I)=.TRUE.
55           PRINTPDB(I)=outpdb
56           printmol2(i)=outmol2
57         ENDIF
58       ENDDO
59       PRINT *,'Number of cutoffs:',NCUT
60       PRINT *,'Cutoff values:'
61       DO ICUT=1,NCUT
62         PRINT '(8HRCUTOFF(,I2,2H)=,F8.1,2i2)',ICUT,RCUTOFF(ICUT),
63      &    printpdb(icut),printmol2(icut)
64       ENDDO
65       DO I=1,NRES-3  
66         MULT(I)=1
67       ENDDO
68       ICON=1
69   123 continue
70       if (from_cart) then
71         if (efree) then
72         read (intin,*,end=13,err=11) energy(icon),totfree(icon),
73      &    rmstab(icon),
74      &    nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
75      &    i=1,nss_all(icon)),iscore(icon)
76         else
77         read (intin,*,end=13,err=11) energy(icon),rmstab(icon),
78      &    nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
79      &    i=1,nss_all(icon)),iscore(icon)
80         endif
81         read (intin,'(8f10.5)',end=13,err=10) 
82      &    ((allcart(j,i,icon),j=1,3),i=1,nres),
83      &    ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
84         print *,icon,energy(icon),nss_all(icon),rmstab(icon)
85       else 
86         read(intin,'(a80)',end=13,err=12) lineh
87         read(lineh(:5),*,err=8) ic
88         if (efree) then
89         read(lineh(6:),*,err=8) energy(icon)
90         else
91         read(lineh(6:),*,err=8) energy(icon)
92         endif
93         goto 9
94     8   ic=1
95         print *,'error, assuming e=1d10',lineh
96         energy(icon)=1d10
97         nss=0
98     9   continue
99 cold        read(lineh(18:),*,end=13,err=11) nss_all(icon)
100         ii = index(lineh(15:)," ")+15
101         read(lineh(ii:),*,end=13,err=11) nss_all(icon)
102         IF (NSS_all(icon).LT.9) THEN
103           read (lineh(20:),*,end=102)
104      &    (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),
105      &    iscore(icon)
106         ELSE
107           read (lineh(20:),*,end=102) 
108      &           (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
109           read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),
110      &      I=9,NSS_all(icon)),iscore(icon)
111         ENDIF
112
113   102   continue  
114
115         PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
116         totfree(icon)=energy(icon)
117         call read_angles(intin,*13)
118         call chainbuild
119         do ii=1,2*nres
120           do jj=1,3
121             allcart(jj,ii,icon)=c(jj,ii)
122           enddo 
123         enddo
124         do i=1,nres
125           phiall(i,icon)=phi(i)
126           thetall(i,icon)=theta(i)
127           alphall(i,icon)=alph(i)
128           omall(i,icon)=omeg(i)
129         enddo
130       endif
131       ICON=ICON+1
132       GOTO 123
133 C      
134 C CALCULATE DISTANCES
135 C
136    10 print *,'something wrong with angles'
137       goto 13
138    11 print *,'something wrong with NSS',nss
139       goto 13
140    12 print *,'something wrong with header'
141
142    13 NCON=ICON-1
143       if (lgrp) then
144         i=0
145         do while (i.lt.ncon)
146 c          read (jstatin,*,err=1111,end=1111) idum,(eee(j,i),j=1,n_ene)
147           read (jstatin,'(a)') lineh
148           if (index(lineh,'#').gt.0) goto 1123
149           i=i+1
150           read (lineh,*) idum,(eee(j,i),j=1,n_ene)
151           print *,i,idum,(eee(j,i),j=1,n_ene)
152           energy(i)=eee(15,i) 
153  1123     continue
154         enddo
155       endif 
156       goto 1112
157  1111 print *,'Error in the statin file; statout will not be produced.'
158       write (iout,*)
159      &  'Error in the statin file; statout will not be produced.'
160       lgrp=.false.
161  1112 continue
162       print *,'ncon',ncon
163       DO I=1,NCON
164         ICC(I)=I
165       ENDDO
166       WRITE (iout,'(A80)') TITEL
167       t1=tcpu()
168       DO I=1,NCON-1
169         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
170         if (from_cart) then
171           do ii=1,2*nres
172             do jj=1,3
173               c(jj,ii)=allcart(jj,ii,i)
174             enddo 
175           enddo
176           call int_from_cart(.true.,.false.)
177           do ii=1,nres
178             phiall(ii,i)=phi(ii)
179             thetall(ii,i)=theta(ii)
180             alphall(ii,i)=alph(ii)
181             omall(ii,i)=omeg(ii)
182           enddo
183         else
184           do ii=1,nres
185             phi(ii)=phiall(ii,i)
186             theta(ii)=thetall(ii,i)
187             alph(ii)=alphall(ii,i)
188             omeg(ii)=omall(ii,i)
189           enddo
190           call chainbuild
191         endif
192         do ii=1,nres
193           do jj=1,3
194             cref(jj,ii)=c(jj,ii)
195           enddo
196         enddo
197         DO J=I+1,NCON
198           IND=IOFFSET(NCON,I,J)
199           ATTALUMS(IND)=DIFCONF(I,J)
200 c          write (iout,'(2i4,f10.5)') i,j,ATTALUMS(IND)
201           DISS(IND)=ATTALUMS(IND)
202         ENDDO
203       ENDDO
204       t2=tcpu()
205       WRITE (iout,'(/a,1pe14.5,a/)') 
206      & 'Time for distance calculation:',T2-T1,' sec.'
207       t1=tcpu()
208       PRINT '(a)','End of distance computation'
209
210       if (punch_dist) then
211         do i=1,ncon-1
212           do j=i+1,ncon
213             IND=IOFFSET(NCON,I,J)
214             write (jrms,'(2i5,2f10.5)') i,j,attalums(IND),
215      &        energy(j)-energy(i)
216           enddo
217         enddo
218       endif
219 C
220 C Print out the RMS deviation matrix.
221 C
222       if (print_dist) CALL DISTOUT(NCON)
223 C
224 C  call hierarchical clustering HC from F. Murtagh
225 C
226       N=NCON
227       LEN = (N*(N-1))/2
228       write(iout,*) "-------------------------------------------"
229       write(iout,*) "HIERARCHICAL CLUSTERING using"
230       if (iopt.eq.1) then
231         write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
232       elseif (iopt.eq.2) then
233         write(iout,*) "SINGLE LINK METHOD"
234       elseif (iopt.eq.3) then
235         write(iout,*) "COMPLETE LINK METHOD"
236       elseif (iopt.eq.4) then
237         write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
238       elseif (iopt.eq.5) then
239         write(iout,*) "MCQUITTY'S METHOD"
240       elseif (iopt.eq.6) then
241         write(iout,*) "MEDIAN (GOWER'S) METHOD"
242       elseif (iopt.eq.7) then
243         write(iout,*) "CENTROID METHOD"
244       else
245         write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
246         write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
247         stop
248       endif
249       write(iout,*)
250       write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
251       write(iout,*) "February 1986"
252       write(iout,*) "References:"
253       write(iout,*) "1. Multidimensional clustering algorithms"
254       write(iout,*) "   Fionn Murtagh"
255       write(iout,*) "   Vienna : Physica-Verlag, 1985."
256       write(iout,*) "2. Multivariate data analysis"
257       write(iout,*) "   Fionn Murtagh and Andre Heck"
258       write(iout,*) "   Kluwer Academic Publishers, 1987"
259       write(iout,*) "-------------------------------------------"
260       write(iout,*)
261
262       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
263       LEV = N-1
264       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
265       CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
266
267       icut=1
268       i=1
269       NGR=i+1
270       do j=1,n
271         licz(iclass(j,i))=licz(iclass(j,i))+1
272         nconf(iclass(j,i),licz(iclass(j,i)))=j
273       enddo        
274       do i=1,lev-1
275
276          idum=lev-i
277          DO L=1,LEV
278             IF (HEIGHT(L).EQ.IDUM) GOTO 190
279          ENDDO
280  190     IDUM=L
281 cd         print *,i+1,CRITVAL(IDUM)
282          IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
283           WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
284           write (iout,'(a,f8.2)') 'Maximum distance found:',
285      &              CRITVAL(IDUM)
286           CALL SRTCLUST(ICUT,ncon)
287           CALL TRACK(ICUT)
288           CALL WRTCLUST(ncon,ICUT,PRINTANG,PRINTPDB,PRINTMOL2)
289           icut=icut+1
290           if (icut.gt.ncut) goto 191
291          ENDIF
292          NGR=i+1
293          do l=1,maxgr
294           licz(l)=0
295          enddo
296          do j=1,n
297           licz(iclass(j,i))=licz(iclass(j,i))+1
298           nconf(iclass(j,i),licz(iclass(j,i)))=j
299 cd          print *,j,iclass(j,i),
300 cd     &     licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
301          enddo
302       enddo
303  191  continue
304 C
305       if (plot_tree) then
306         CALL WRITRACK
307         CALL PLOTREE
308       endif
309 C
310       t2=tcpu()
311       WRITE (iout,'(/a,1pe14.5,a/)') 
312      & 'Total time for clustering:',T2-T1,' sec.'
313       if (lgrp) then
314         write (jstatout,'(a5,18a12,10f4.1)')"#    ",
315      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p",
316      &   "EBE bending","ESC SCloc","ETORS ",
317      &   "ECORR4 ","ECORR5 ","ECORR6 ",
318      &   "EELLO ","ETURN3 ","ETURN4 ","ETURN6 "
319      &   ,"ETOT total","RMSD","nat.contact","nnt.contact",
320      &   "cont.order",(rcutoff(i),i=1,ncut)
321         do i=1,ncon
322           write(jstatout,'(i5,18f12.4,10i6)') i,(eee(j,i),j=1,n_ene),
323      &     (iass_tot(i,icut),icut=1,ncut)
324         enddo
325       endif
326 C
327       stop '********** Program terminated normally.'
328       end
329 c---------------------------------------------------------------------------
330       double precision function difconf(icon,jcon)
331       include 'DIMENSIONS'
332       include 'sizesclu.dat'
333       include 'COMMON.CONTROL'
334       include 'COMMON.CLUSTER'
335       include 'COMMON.CHAIN' 
336       include 'COMMON.INTERACT'
337       include 'COMMON.VAR'
338       logical non_conv
339       double precision przes(3),obrot(3,3)
340       double precision xx(3,maxres2),yy(3,maxres2)
341       do i=1,nres
342         phi(i)=phiall(i,jcon)
343         theta(i)=thetall(i,jcon)
344         alph(i)=alphall(i,jcon)
345         omeg(i)=omall(i,jcon)
346       enddo
347       call chainbuild
348 c     do i=1,nres
349 c       print '(i4,2(3f10.5,5x))',i,(cref(j,i),j=1,3),(c(j,i),j=1,3)
350 c     enddo
351       if (lside) then
352         ii=0
353         do i=nnt,nct
354           ii=ii+1
355           do j=1,3
356             xx(j,ii)=allcart(j,i,jcon)
357             yy(j,ii)=cref(j,i)
358           enddo
359         enddo
360         do i=nnt,nct
361 c          if (itype(i).ne.10) then
362             ii=ii+1
363             do j=1,3 
364               xx(j,ii)=allcart(j,i+nres,jcon)
365               yy(j,ii)=cref(j,i+nres)
366             enddo
367 c          endif
368         enddo
369         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
370       else
371         do i=nnt,nct
372           do j=1,3
373             c(j,i)=allcart(j,i,jcon)
374           enddo
375 c          write (2,'(i5,3f10.5,5x,3f10.5)') i,
376 c     &     (c(j,i),j=1,3),(cref(j,i),j=1,3)
377         enddo
378         call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obrot
379      &       ,non_conv)
380       endif
381       if (rms.lt.0.0) then
382         print *,error,'rms^2 = ',rms,icon,jcon
383         stop
384       endif
385       if (non_conv) print *,non_conv,icon,jcon
386       difconf=dsqrt(rms)
387       RETURN
388       END
389 C------------------------------------------------------------------------------
390       double precision function rmsnat(jcon)
391       include 'DIMENSIONS'
392       include 'sizesclu.dat'
393       include 'COMMON.CONTROL'
394       include 'COMMON.CLUSTER'
395       include 'COMMON.CHAIN' 
396       include 'COMMON.INTERACT'
397       include 'COMMON.VAR'
398       logical non_conv
399       double precision przes(3),obrot(3,3)
400       double precision xx(3,maxres2),yy(3,maxres2)
401       if (lside) then
402         ii=0
403         do i=nnt,nct
404           ii=ii+1
405           do j=1,3
406             xx(j,ii)=allcart(j,i,jcon)
407             yy(j,ii)=cref_pdb(j,i)
408           enddo
409         enddo
410         do i=nnt,nct
411 c          if (itype(i).ne.10) then
412             ii=ii+1
413             do j=1,3 
414               xx(j,ii)=allcart(j,i+nres,jcon)
415               yy(j,ii)=cref_pdb(j,i+nres)
416             enddo
417 c          endif
418         enddo
419         call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
420       else
421         do i=nnt,nct
422           do j=1,3
423             c(j,i)=allcart(j,i,jcon)
424           enddo
425         enddo
426         call fitsq(rms,c(1,nnt),cref_pdb(1,nnt),nct-nnt+1,przes,obrot
427      &       ,non_conv)
428       endif
429       if (rms.lt.0.0) then
430         print *,error,'rms^2 = ',rms,icon,jcon
431         stop
432       endif
433       if (non_conv) print *,non_conv,icon,jcon
434       rmsnat=dsqrt(rms)
435       RETURN
436       END
437 C------------------------------------------------------------------------------
438       subroutine distout(ncon)
439       include 'DIMENSIONS'
440       include 'sizesclu.dat'
441       parameter (ncol=10)
442       include 'COMMON.IOUNITS'
443       include 'COMMON.CLUSTER'
444       dimension b(ncol)
445       write (iout,'(a)') 'The distance matrix'
446       do 1 i=1,ncon,ncol
447       nlim=min0(i+ncol-1,ncon)
448       write (iout,1000) (k,k=i,nlim)
449       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
450  1000 format (/8x,10(i4,3x))
451  1020 format (/1x,80(1h-)/)
452       do 2 j=i,ncon
453       jlim=min0(j,nlim)
454       if (jlim.eq.j) then
455         b(jlim-i+1)=0.0d0
456         jlim1=jlim-1
457       else
458         jlim1=jlim
459       endif
460       do 3 k=i,jlim1
461        if (j.lt.k) then 
462           IND=IOFFSET(NCON,j,k)
463        else
464           IND=IOFFSET(NCON,k,j)
465        endif
466     3  b(k-i+1)=attalums(IND)
467       write (iout,1010) j,(b(k),k=1,jlim-i+1)
468     2 continue
469     1 continue
470  1010 format (i5,3x,10(f6.2,1x))
471       return
472       end
473 C-----------------------------------------------------------------------
474       block data
475       include 'DIMENSIONS'
476       include 'COMMON.LOCAL'
477       data dsc /  1.237,  ! CYS (type  1)
478      &            2.142,  ! MET (type  2)
479      &            2.299,  ! PHE (type  3)
480      &            1.776,  ! ILE (type  4)
481      &            1.939,  ! LEU (type  5)
482      &            1.410,  ! VAL (type  6)
483      &            2.605,  ! TRP (type  7)
484      &            2.484,  ! TYR (type  8)
485      &            0.743,  ! ALA (type  9)
486      &            0.000,  ! GLY (type 10)
487      &            1.393,  ! THR (type 11)
488      &            1.150,  ! SER (type 12)
489      &            2.240,  ! GLN (type 13)
490      &            1.684,  ! ASN (type 14)
491      &            2.254,  ! GLU (type 15)
492      &            1.709,  ! ASP (type 16)
493      &            2.113,  ! HIS (type 17)
494      &            3.020,  ! ARG (type 18)
495      &            2.541,  ! LYS (type 19)
496      &            1.345,  ! PRO (type 20)
497      &            0.000  /! D   (type 21)
498       end