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)
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,kkk
38 double precision t1,t2,tcpu,difconf
40 double precision varia(maxvar)
41 double precision hrtime,mintime,sectime
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 )
49 write(iout,*) "SEVERE ERROR - Can't initialize MPI."
50 call mpi_finalize(ierror)
53 if (nprocs.gt.MaxProcs+1) then
54 write (2,*) "Error - too many processors",
56 write (2,*) "Increase MaxProcs and recompile"
57 call MPI_Finalize(IERROR)
68 c if (refstr) call read_ref_structure(*30)
75 c write (iout,*) "Before permut"
76 c write (iout,*) "symetr", symetr
79 c write (iout,*) "after permut"
81 print *,'MAIN: nnt=',nnt,' nct=',nct
93 IF (RCUTOFF(I).LT.0.0) THEN
94 RCUTOFF(I)=ABS(RCUTOFF(I))
102 write (iout,*) 'Number of cutoffs:',NCUT
103 write (iout,*) 'Cutoff values:'
105 WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),
106 & printpdb(icut),printmol2(icut)
108 else if (nclust.gt.0) then
109 write (iout,'("Number of clusters requested",i5)') nclust
112 & write (iout,*) "ERROR: Either nclust or ncut must be >0"
121 call read_coords(ncon,*20)
122 write (iout,*) 'from read_coords: ncon',ncon
124 write (iout,*) "nT",nT
126 write (iout,*) "iT",iT
128 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)
166 IND=IOFFSET(NCON_work,I,J)
168 if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
171 DISS(IND1)=DIFCONF(I,J)
172 c write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
179 WRITE (iout,'(/a,1pe14.5,a/)')
180 & 'Time for distance calculation:',T2-T1,' sec.'
182 PRINT '(a)','End of distance computation'
185 call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
186 & scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
187 if (me.eq.master) then
189 open(80,file='/tmp/distance',form='unformatted')
196 IND=IOFFSET(NCON,I,J)
197 write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
198 & energy(j)-energy(i)
203 C Print out the RMS deviation matrix.
205 if (print_dist) CALL DISTOUT(NCON_work)
207 C call hierarchical clustering HC from F. Murtagh
211 write(iout,*) "-------------------------------------------"
212 write(iout,*) "HIERARCHICAL CLUSTERING using"
214 write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
215 elseif (iopt.eq.2) then
216 write(iout,*) "SINGLE LINK METHOD"
217 elseif (iopt.eq.3) then
218 write(iout,*) "COMPLETE LINK METHOD"
219 elseif (iopt.eq.4) then
220 write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
221 elseif (iopt.eq.5) then
222 write(iout,*) "MCQUITTY'S METHOD"
223 elseif (iopt.eq.6) then
224 write(iout,*) "MEDIAN (GOWER'S) METHOD"
225 elseif (iopt.eq.7) then
226 write(iout,*) "CENTROID METHOD"
228 write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
229 write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
233 write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
234 write(iout,*) "February 1986"
235 write(iout,*) "References:"
236 write(iout,*) "1. Multidimensional clustering algorithms"
237 write(iout,*) " Fionn Murtagh"
238 write(iout,*) " Vienna : Physica-Verlag, 1985."
239 write(iout,*) "2. Multivariate data analysis"
240 write(iout,*) " Fionn Murtagh and Andre Heck"
241 write(iout,*) " Kluwer Academic Publishers, 1987"
242 write(iout,*) "-------------------------------------------"
246 write (iout,*) "The TOTFREE array"
248 write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
252 CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
254 write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
256 write (iout,*) "Too few conformations to cluster."
259 CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
260 !c CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
261 !c 3/3/16 AL: added explicit number of cluters
262 if (nclust.gt.0) then
278 licz(iclass(j,i))=licz(iclass(j,i))+1
279 nconf(iclass(j,i),licz(iclass(j,i)))=j
280 !c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
281 !c & nconf(iclass(j,i),licz(iclass(j,i)))
287 IF (HEIGHT(L).EQ.IDUM) GOTO 190
290 write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
291 & " icut",icut," cutoff",rcutoff(icut)
292 IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
293 if (nclust.le.0) WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
294 write (iout,'(a,f8.2)') 'Maximum distance found:',
296 CALL SRTCLUST(ICUT,ncon_work,iT)
298 CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
300 if (icut.gt.ncut) goto 191
307 licz(iclass(j,i))=licz(iclass(j,i))+1
308 nconf(iclass(j,i),licz(iclass(j,i)))=j
309 c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
310 c & nconf(iclass(j,i),licz(iclass(j,i)))
311 cd print *,j,iclass(j,i),
312 cd & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
323 WRITE (iout,'(/a,1pe14.5,a/)')
324 & 'Total time for clustering:',T2-T1,' sec.'
332 close(icbase,status="delete")
334 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
336 stop '********** Program terminated normally.'
337 20 write (iout,*) "Error reading coordinates"
339 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
342 30 write (iout,*) "Error reading reference structure"
344 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
348 c---------------------------------------------------------------------------
349 double precision function difconf(icon,jcon)
352 include 'sizesclu.dat'
353 include 'COMMON.CONTROL'
354 include 'COMMON.CLUSTER'
355 include 'COMMON.CHAIN'
356 include 'COMMON.INTERACT'
358 include 'COMMON.IOUNITS'
360 double precision przes(3),obrot(3,3)
361 double precision xx(3,maxres2),yy(3,maxres2)
362 integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
363 integer iaperm,ibezperm,run
364 double precision rms,rmsmina
365 c write (iout,*) "tu dochodze"
371 c write (iout,*) "nperm",nperm
373 c write (iout,*) "tabperm", tabperm(1,1)
377 chalen=int((nend-nstart+2)/symetr)
379 do i=nstart,(nstart+chalen-1)
381 c write (iout,*) "tutaj",zzz
383 iaperm=(zzz-1)*chalen+i
384 ibezperm=(run-1)*chalen+i
386 xx(j,ii)=allcart(j,iaperm,jcon)
387 yy(j,ii)=cref(j,ibezperm,kkk)
392 do i=nstart,(nstart+chalen-1)
395 iaperm=(zzz-1)*chalen+i
396 ibezperm=(run-1)*chalen+i
397 c if (itype(i).ne.10) then
400 xx(j,ii)=allcart(j,iaperm+nres,jcon)
401 yy(j,ii)=cref(j,ibezperm+nres,kkk)
406 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
408 chalen=int((nct-nnt+2)/symetr)
410 do i=nnt,(nnt+chalen-1)
412 c write (iout,*) "tu szukaj", zzz,run,kkk
413 iaperm=(zzz-1)*chalen+i
414 ibezperm=(run-1)*chalen+i
417 c(j,i)=allcart(j,iaperm,jcon)
421 call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,
426 print *,'error, rms^2 = ',rms,icon,jcon
429 if (non_conv) print *,non_conv,icon,jcon
430 if (rmsmina.gt.rms) rmsmina=rms
432 difconf=dsqrt(rmsmina)
435 C------------------------------------------------------------------------------
436 subroutine distout(ncon)
439 include 'sizesclu.dat'
442 include 'COMMON.IOUNITS'
443 include 'COMMON.CLUSTER'
444 integer i,j,k,jlim,jlim1,nlim,ind,ioffset
447 write (iout,'(a)') 'The distance matrix'
449 nlim=min0(i+ncol-1,ncon)
450 write (iout,1000) (k,k=i,nlim)
451 write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
452 1000 format (/8x,10(i4,3x))
453 1020 format (/1x,80(1h-)/)
464 IND=IOFFSET(NCON,j,k)
466 IND=IOFFSET(NCON,k,j)
469 write (iout,1010) j,(b(k),k=1,jlim-i+1)
472 1010 format (i5,3x,10(f6.2,1x))