2 C Program to cluster united-residue MCM results.
7 include 'COMMON.INTERACT'
10 include 'COMMON.HEADER'
11 include 'COMMON.CONTROL'
12 include 'COMMON.CONTACTS'
13 include 'COMMON.CHAIN'
15 include 'COMMON.CLUSTER'
16 include 'COMMON.IOUNITS'
17 logical printang(max_cut)
18 integer printpdb(max_cut)
19 integer printmol2(max_cut)
21 REAL CRIT(maxconf),MEMBR(maxconf)
22 REAL CRITVAL(maxconf-1)
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)
31 double precision varia(maxvar)
32 double precision hrtime,mintime,sectime
34 double precision eee(n_ene,maxconf)
47 print *,'MAIN: nnt=',nnt,' nct=',nct
53 IF (RCUTOFF(I).LT.0.0) THEN
54 RCUTOFF(I)=ABS(RCUTOFF(I))
60 PRINT *,'Number of cutoffs:',NCUT
61 PRINT *,'Cutoff values:'
63 PRINT '(8HRCUTOFF(,I2,2H)=,F8.1,2i2)',ICUT,RCUTOFF(ICUT),
64 & printpdb(icut),printmol2(icut)
71 call cxread(ncon,*101)
73 101 write (iout,*) "Error in CX file"
83 read (intin,*,end=13,err=11) energy(icon),totfree(icon),
85 & nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
86 & i=1,nss_all(icon)),iscore(icon)
88 read (intin,*,end=13,err=11) energy(icon),rmstab(icon),
89 & nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
90 & i=1,nss_all(icon)),iscore(icon)
92 read (intin,'(8f10.5)',end=13,err=10)
93 & ((allcart(j,i,icon),j=1,3),i=1,nres),
94 & ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
95 print *,icon,energy(icon),nss_all(icon),rmstab(icon)
97 read(intin,'(a80)',end=13,err=12) lineh
98 read(lineh(:5),*,err=8) ic
100 read(lineh(6:),*,err=8) energy(icon)
102 read(lineh(6:),*,err=8) energy(icon)
106 print *,'error, assuming e=1d10',lineh
110 cold read(lineh(18:),*,end=13,err=11) nss_all(icon)
111 ii = index(lineh(15:)," ")+15
112 read(lineh(ii:),*,end=13,err=11) nss_all(icon)
113 IF (NSS_all(icon).LT.9) THEN
114 read (lineh(20:),*,end=102)
115 & (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),
118 read (lineh(20:),*,end=102)
119 & (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
120 read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),
121 & I=9,NSS_all(icon)),iscore(icon)
126 PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
127 totfree(icon)=energy(icon)
128 call read_angles(intin,*13)
132 allcart(jj,ii,icon)=c(jj,ii)
136 phiall(i,icon)=phi(i)
137 thetall(i,icon)=theta(i)
138 alphall(i,icon)=alph(i)
139 omall(i,icon)=omeg(i)
144 10 print *,'something wrong with angles'
146 11 print *,'something wrong with NSS',nss
148 12 print *,'something wrong with header'
156 write (iout,*) "Conformation",icon
158 write (iout,'(i5,3f10.5,5x,3f10.5)') i,
159 & (allcart(j,i,icon),j=1,3),(allcart(j,i+nres,icon),j=1,3)
166 c read (jstatin,*,err=1111,end=1111) idum,(eee(j,i),j=1,n_ene)
167 read (jstatin,'(a)') lineh
168 if (index(lineh,'#').gt.0) goto 1123
170 read (lineh,*) idum,(eee(j,i),j=1,n_ene)
171 print *,i,idum,(eee(j,i),j=1,n_ene)
177 1111 print *,'Error in the statin file; statout will not be produced.'
179 & 'Error in the statin file; statout will not be produced.'
186 WRITE (iout,'(A80)') TITEL
189 C CALCULATE DISTANCES
192 if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
193 if (from_cart .or. from_cx) then
196 c(jj,ii)=allcart(jj,ii,i)
199 call int_from_cart(.true.,.false.)
202 thetall(ii,i)=theta(ii)
203 alphall(ii,i)=alph(ii)
209 theta(ii)=thetall(ii,i)
210 alph(ii)=alphall(ii,i)
221 IND=IOFFSET(NCON,I,J)
222 ATTALUMS(IND)=DIFCONF(I,J)
223 c write (iout,'(2i4,f10.5)') i,j,ATTALUMS(IND)
224 DISS(IND)=ATTALUMS(IND)
228 WRITE (iout,'(/a,1pe14.5,a/)')
229 & 'Time for distance calculation:',T2-T1,' sec.'
231 PRINT '(a)','End of distance computation'
236 IND=IOFFSET(NCON,I,J)
237 write (jrms,'(2i5,2f10.5)') i,j,attalums(IND),
238 & energy(j)-energy(i)
243 C Print out the RMS deviation matrix.
245 if (print_dist) CALL DISTOUT(NCON)
247 C call hierarchical clustering HC from F. Murtagh
251 write(iout,*) "-------------------------------------------"
252 write(iout,*) "HIERARCHICAL CLUSTERING using"
254 write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
255 elseif (iopt.eq.2) then
256 write(iout,*) "SINGLE LINK METHOD"
257 elseif (iopt.eq.3) then
258 write(iout,*) "COMPLETE LINK METHOD"
259 elseif (iopt.eq.4) then
260 write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
261 elseif (iopt.eq.5) then
262 write(iout,*) "MCQUITTY'S METHOD"
263 elseif (iopt.eq.6) then
264 write(iout,*) "MEDIAN (GOWER'S) METHOD"
265 elseif (iopt.eq.7) then
266 write(iout,*) "CENTROID METHOD"
268 write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
269 write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
273 write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
274 write(iout,*) "February 1986"
275 write(iout,*) "References:"
276 write(iout,*) "1. Multidimensional clustering algorithms"
277 write(iout,*) " Fionn Murtagh"
278 write(iout,*) " Vienna : Physica-Verlag, 1985."
279 write(iout,*) "2. Multivariate data analysis"
280 write(iout,*) " Fionn Murtagh and Andre Heck"
281 write(iout,*) " Kluwer Academic Publishers, 1987"
282 write(iout,*) "-------------------------------------------"
284 CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
287 CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
290 c print *,"call HCDEN"
291 CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
298 licz(iclass(j,i))=licz(iclass(j,i))+1
299 nconf(iclass(j,i),licz(iclass(j,i)))=j
305 IF (HEIGHT(L).EQ.IDUM) GOTO 190
308 cd print *,i+1,CRITVAL(IDUM)
309 IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
310 WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
311 write (iout,'(a,f8.2)') 'Maximum distance found:',
313 CALL SRTCLUST(ICUT,ncon)
315 CALL WRTCLUST(ncon,ICUT,PRINTANG,PRINTPDB,PRINTMOL2)
317 if (icut.gt.ncut) goto 191
324 licz(iclass(j,i))=licz(iclass(j,i))+1
325 nconf(iclass(j,i),licz(iclass(j,i)))=j
326 cd print *,j,iclass(j,i),
327 cd & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
338 WRITE (iout,'(/a,1pe14.5,a/)')
339 & 'Total time for clustering:',T2-T1,' sec.'
341 write (jstatout,'(a5,18a12,10f4.1)')"# ",
342 & "EVDW SC-SC","EVDW2 SC-p","EES p-p",
343 & "EBE bending","ESC SCloc","ETORS ",
344 & "ECORR4 ","ECORR5 ","ECORR6 ",
345 & "EELLO ","ETURN3 ","ETURN4 ","ETURN6 "
346 & ,"ETOT total","RMSD","nat.contact","nnt.contact",
347 & "cont.order",(rcutoff(i),i=1,ncut)
349 write(jstatout,'(i5,18f12.4,10i6)') i,(eee(j,i),j=1,n_ene),
350 & (iass_tot(i,icut),icut=1,ncut)
354 stop '********** Program terminated normally.'
356 c---------------------------------------------------------------------------
357 double precision function difconf(icon,jcon)
359 include 'sizesclu.dat'
360 include 'COMMON.CONTROL'
361 include 'COMMON.CLUSTER'
362 include 'COMMON.CHAIN'
363 include 'COMMON.INTERACT'
366 double precision przes(3),obrot(3,3)
367 double precision xx(3,maxres2),yy(3,maxres2)
369 c phi(i)=phiall(i,jcon)
370 c theta(i)=thetall(i,jcon)
371 c alph(i)=alphall(i,jcon)
372 c omeg(i)=omall(i,jcon)
376 c print '(i4,2(3f10.5,5x))',i,(cref(j,i),j=1,3),(c(j,i),j=1,3)
383 xx(j,ii)=allcart(j,i,jcon)
384 yy(j,ii)=allcart(j,i,icon)
388 c if (itype(i).ne.10) then
391 xx(j,ii)=allcart(j,i+nres,jcon)
392 yy(j,ii)=allcart(j,i+nres,icon)
396 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
400 c c(j,i)=allcart(j,i,jcon)
402 c write (2,'(i5,3f10.5,5x,3f10.5)') i,
403 c & (c(j,i),j=1,3),(cref(j,i),j=1,3)
405 call fitsq(rms,allcart(1,nnt,jcon),allcart(1,nnt,icon),
406 & nct-nnt+1,przes,obrot,non_conv)
409 print *,error,'rms^2 = ',rms,icon,jcon
412 if (non_conv) print *,non_conv,icon,jcon
416 C------------------------------------------------------------------------------
417 double precision function rmsnat(jcon)
419 include 'sizesclu.dat'
420 include 'COMMON.CONTROL'
421 include 'COMMON.CLUSTER'
422 include 'COMMON.CHAIN'
423 include 'COMMON.INTERACT'
426 double precision przes(3),obrot(3,3)
427 double precision xx(3,maxres2),yy(3,maxres2)
433 xx(j,ii)=allcart(j,i,jcon)
434 yy(j,ii)=cref_pdb(j,i)
438 c if (itype(i).ne.10) then
441 xx(j,ii)=allcart(j,i+nres,jcon)
442 yy(j,ii)=cref_pdb(j,i+nres)
446 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
450 c(j,i)=allcart(j,i,jcon)
453 call fitsq(rms,c(1,nnt),cref_pdb(1,nnt),nct-nnt+1,przes,obrot
457 print *,error,'rms^2 = ',rms,icon,jcon
460 if (non_conv) print *,non_conv,icon,jcon
464 C------------------------------------------------------------------------------
465 subroutine distout(ncon)
467 include 'sizesclu.dat'
469 include 'COMMON.IOUNITS'
470 include 'COMMON.CLUSTER'
472 write (iout,'(a)') 'The distance matrix'
474 nlim=min0(i+ncol-1,ncon)
475 write (iout,1000) (k,k=i,nlim)
476 write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
477 1000 format (/8x,10(i4,3x))
478 1020 format (/1x,80(1h-)/)
489 IND=IOFFSET(NCON,j,k)
491 IND=IOFFSET(NCON,k,j)
493 3 b(k-i+1)=attalums(IND)
494 write (iout,1010) j,(b(k),k=1,jlim-i+1)
497 1010 format (i5,3x,10(f6.2,1x))
500 C-----------------------------------------------------------------------
503 include 'COMMON.LOCAL'
504 data dsc / 1.237, ! CYS (type 1)
505 & 2.142, ! MET (type 2)
506 & 2.299, ! PHE (type 3)
507 & 1.776, ! ILE (type 4)
508 & 1.939, ! LEU (type 5)
509 & 1.410, ! VAL (type 6)
510 & 2.605, ! TRP (type 7)
511 & 2.484, ! TYR (type 8)
512 & 0.743, ! ALA (type 9)
513 & 0.000, ! GLY (type 10)
514 & 1.393, ! THR (type 11)
515 & 1.150, ! SER (type 12)
516 & 2.240, ! GLN (type 13)
517 & 1.684, ! ASN (type 14)
518 & 2.254, ! GLU (type 15)
519 & 1.709, ! ASP (type 16)
520 & 2.113, ! HIS (type 17)
521 & 3.020, ! ARG (type 18)
522 & 2.541, ! LYS (type 19)
523 & 1.345, ! PRO (type 20)
524 & 0.000 /! D (type 21)