2 C Program to cluster united-residue MCM results.
9 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
12 include 'COMMON.TIME1'
13 include 'COMMON.INTERACT'
14 include 'COMMON.NAMES'
16 include 'COMMON.HEADER'
17 include 'COMMON.CONTROL'
18 include 'COMMON.CHAIN'
20 include 'COMMON.CLUSTER'
21 include 'COMMON.IOUNITS'
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)
34 DIMENSION NN(maxconf),DISNN(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
41 double precision varia(maxvar)
42 double precision hrtime,mintime,sectime
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 )
51 write(iout,*) "SEVERE ERROR - Can't initialize MPI."
52 call mpi_finalize(ierror)
55 if (nprocs.gt.MaxProcs+1) then
56 write (2,*) "Error - too many processors",
58 write (2,*) "Increase MaxProcs and recompile"
59 call MPI_Finalize(IERROR)
69 c if (refstr) call read_ref_structure(*30)
76 c write (iout,*) "Before permut"
77 c write (iout,*) "symetr", symetr
80 c write (iout,*) "after permut"
82 print *,'MAIN: nnt=',nnt,' nct=',nct
94 IF (RCUTOFF(I).LT.0.0) THEN
95 RCUTOFF(I)=ABS(RCUTOFF(I))
103 write (iout,*) 'Number of cutoffs:',NCUT
104 write (iout,*) 'Cutoff values:'
106 WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
107 & printpdb(icut),printmol2(icut)
109 else if (nclust.gt.0) then
110 write (iout,'("Number of clusters requested",i5)') nclust
113 & write (iout,*) "ERROR: Either nclust or ncut must be >0"
122 call read_coords(ncon,*20)
123 write (iout,*) 'from read_coords: ncon',ncon
125 write (iout,*) "nT",nT
127 write (iout,*) "iT",iT
129 call work_partition(.true.,ncon)
131 call probabl(iT,ncon_work,ncon,*20)
133 if (ncon_work.lt.2) then
134 write (iout,*) "Too few conformations; clustering skipped"
138 ndis=ncon_work*(ncon_work-1)/2
139 call work_partition(.true.,ndis)
145 WRITE (iout,'(A80)') TITEL
148 C CALCULATE DISTANCES
150 call daread_ccoords(1,ncon_work)
153 if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
156 c(l,k)=allcart(l,k,i)
165 IND=IOFFSET(NCON_work,I,J)
167 if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
171 DISS_(IND1)=DIFCONF(I,J)
173 DISS(IND1)=DIFCONF(I,J)
175 c write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
182 WRITE (iout,'(/a,1pe14.5,a/)')
183 & 'Time for distance calculation:',T2-T1,' sec.'
185 PRINT '(a)','End of distance computation'
188 call MPI_Gatherv(diss_(1),scount(me),MPI_REAL,diss(1),
189 & scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
190 if (me.eq.master) then
192 scrachdir2d=scratchdir(:ilen(scratchdir))//'distance'
193 open(80,file=scrachdir2d,form='unformatted')
200 IND=IOFFSET(NCON,I,J)
201 write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
202 & energy(j)-energy(i)
207 C Print out the RMS deviation matrix.
209 if (print_dist) CALL DISTOUT(NCON_work)
211 C call hierarchical clustering HC from F. Murtagh
215 write(iout,*) "-------------------------------------------"
216 write(iout,*) "HIERARCHICAL CLUSTERING using"
218 write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
219 elseif (iopt.eq.2) then
220 write(iout,*) "SINGLE LINK METHOD"
221 elseif (iopt.eq.3) then
222 write(iout,*) "COMPLETE LINK METHOD"
223 elseif (iopt.eq.4) then
224 write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
225 elseif (iopt.eq.5) then
226 write(iout,*) "MCQUITTY'S METHOD"
227 elseif (iopt.eq.6) then
228 write(iout,*) "MEDIAN (GOWER'S) METHOD"
229 elseif (iopt.eq.7) then
230 write(iout,*) "CENTROID METHOD"
232 write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
233 write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
237 write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
238 write(iout,*) "February 1986"
239 write(iout,*) "References:"
240 write(iout,*) "1. Multidimensional clustering algorithms"
241 write(iout,*) " Fionn Murtagh"
242 write(iout,*) " Vienna : Physica-Verlag, 1985."
243 write(iout,*) "2. Multivariate data analysis"
244 write(iout,*) " Fionn Murtagh and Andre Heck"
245 write(iout,*) " Kluwer Academic Publishers, 1987"
246 write(iout,*) "-------------------------------------------"
250 write (iout,*) "The TOTFREE array"
252 write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
256 CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
258 write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
260 write (iout,*) "Too few conformations to cluster."
263 CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
264 c CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
266 c 3/3/16 AL: added explicit number of cluters
267 if (nclust.gt.0) then
282 licz(iclass(j,i))=licz(iclass(j,i))+1
283 nconf(iclass(j,i),licz(iclass(j,i)))=j
284 c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
285 c & nconf(iclass(j,i),licz(iclass(j,i)))
291 IF (HEIGHT(L).EQ.IDUM) GOTO 190
294 c write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
295 c & " icut",icut," cutoff",rcutoff(icut)
296 IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
298 & WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
299 write (iout,'(a,f8.2)') 'Maximum distance found:',
301 CALL SRTCLUST(ICUT,ncon_work,iT)
303 CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
305 if (icut.gt.ncut) goto 191
313 licz(iclass(j,ii))=licz(iclass(j,ii))+1
314 nconf(iclass(j,ii),licz(iclass(j,ii)))=j
315 c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
316 c & nconf(iclass(j,i),licz(iclass(j,i)))
317 cd print *,j,iclass(j,i),
318 cd & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
329 WRITE (iout,'(/a,1pe14.5,a/)')
330 & 'Total time for clustering:',T2-T1,' sec.'
338 close(icbase,status="delete")
340 call MPI_Finalize(IERROR)
342 stop '********** Program terminated normally.'
343 20 write (iout,*) "Error reading coordinates"
345 call MPI_Finalize(IERROR)
348 30 write (iout,*) "Error reading reference structure"
350 call MPI_Finalize(IERROR)
354 c---------------------------------------------------------------------------
355 double precision function difconf(icon,jcon)
358 include 'sizesclu.dat'
359 include 'COMMON.CONTROL'
360 include 'COMMON.CLUSTER'
361 include 'COMMON.CHAIN'
362 include 'COMMON.INTERACT'
364 include 'COMMON.IOUNITS'
366 double precision przes(3),obrot(3,3)
367 double precision xx(3,maxres2),yy(3,maxres2)
368 integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
369 integer iaperm,ibezperm,run
370 double precision rms,rmsmina
371 c write (iout,*) "tu dochodze"
377 c write (iout,*) "nperm",nperm
379 c write (iout,*) "tabperm", tabperm(1,1)
383 chalen=int((nend-nstart+2)/symetr)
385 do i=nstart,(nstart+chalen-1)
387 c write (iout,*) "tutaj",zzz
389 iaperm=(zzz-1)*chalen+i
390 ibezperm=(run-1)*chalen+i
392 xx(j,ii)=allcart(j,iaperm,jcon)
393 yy(j,ii)=cref(j,ibezperm)
398 do i=nstart,(nstart+chalen-1)
401 iaperm=(zzz-1)*chalen+i
402 ibezperm=(run-1)*chalen+i
403 c if (itype(i).ne.10) then
406 xx(j,ii)=allcart(j,iaperm+nres,jcon)
407 yy(j,ii)=cref(j,ibezperm+nres)
412 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
414 chalen=int((nct-nnt+2)/symetr)
416 do i=nnt,(nnt+chalen-1)
418 c write (iout,*) "tu szukaj", zzz,run,kkk
419 iaperm=(zzz-1)*chalen+i
420 ibezperm=(run-1)*chalen+i
423 c(j,i)=allcart(j,iaperm,jcon)
427 call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
431 print *,'error, rms^2 = ',rms,icon,jcon
434 if (non_conv) print *,non_conv,icon,jcon
435 if (rmsmina.gt.rms) rmsmina=rms
437 difconf=dsqrt(rmsmina)
440 C------------------------------------------------------------------------------
441 subroutine distout(ncon)
444 include 'sizesclu.dat'
447 include 'COMMON.IOUNITS'
448 include 'COMMON.CLUSTER'
449 integer i,j,k,jlim,jlim1,nlim,ind,ioffset
452 write (iout,'(a)') 'The distance matrix'
454 nlim=min0(i+ncol-1,ncon)
455 write (iout,1000) (k,k=i,nlim)
456 write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
457 1000 format (/8x,10(i4,3x))
458 1020 format (/1x,80(1h-)/)
469 IND=IOFFSET(NCON,j,k)
471 IND=IOFFSET(NCON,k,j)
474 write (iout,1010) j,(b(k),k=1,jlim-i+1)
477 1010 format (i5,3x,10(f6.2,1x))