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)
32 integer nn,ndis,scount_buf
33 real*4 DISNN, diss_buf(maxdist)
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, ijk, is,ie
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 write (iout,*) "Main: refstr ",refstr
69 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 c 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,*) "Temperature",1.0d0/(beta_h(iT)*1.987D-3)
128 call work_partition(.true.,ncon)
130 call probabl(iT,ncon_work,ncon,*20)
132 if (ncon_work.lt.2) then
133 write (iout,*) "Too few conformations; clustering skipped"
137 ndis=ncon_work*(ncon_work-1)/2
138 call work_partition(.true.,ndis)
143 WRITE (iout,'(A80)') TITEL
146 C CALCULATE DISTANCES
148 call daread_ccoords(1,ncon_work)
151 c if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
154 c(l,k)=allcart(l,k,i)
164 IND=IOFFSET(NCON_work,I,J)
166 if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
169 DISS(IND1)=DIFCONF(I,J)
170 c write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
177 WRITE (iout,'(/a,1pe14.5,a/)')
178 & 'Time for distance calculation:',T2-T1,' sec.'
180 c PRINT '(a)','End of distance computation'
182 scount_buf=scount(me)
185 diss_buf(ijk)=diss(ijk)
190 WRITE (iout,*) "Wchodze do call MPI_Gatherv"
191 call MPI_Gatherv(diss_buf(1),scount_buf,MPI_REAL,diss(1),
192 & scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
193 if (me.eq.master) then
195 open(80,file='/tmp/distance',form='unformatted')
202 IND=IOFFSET(NCON,I,J)
203 write (jrms,'(2i5,2f10.5)') i,j,diss(IND),
204 & energy(j)-energy(i)
209 C Print out the RMS deviation matrix.
211 if (print_dist) CALL DISTOUT(NCON_work)
213 C call hierarchical clustering HC from F. Murtagh
217 write(iout,*) "-------------------------------------------"
218 write(iout,*) "HIERARCHICAL CLUSTERING using"
220 write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
221 elseif (iopt.eq.2) then
222 write(iout,*) "SINGLE LINK METHOD"
223 elseif (iopt.eq.3) then
224 write(iout,*) "COMPLETE LINK METHOD"
225 elseif (iopt.eq.4) then
226 write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
227 elseif (iopt.eq.5) then
228 write(iout,*) "MCQUITTY'S METHOD"
229 elseif (iopt.eq.6) then
230 write(iout,*) "MEDIAN (GOWER'S) METHOD"
231 elseif (iopt.eq.7) then
232 write(iout,*) "CENTROID METHOD"
234 write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
235 write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
239 write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
240 write(iout,*) "February 1986"
241 write(iout,*) "References:"
242 write(iout,*) "1. Multidimensional clustering algorithms"
243 write(iout,*) " Fionn Murtagh"
244 write(iout,*) " Vienna : Physica-Verlag, 1985."
245 write(iout,*) "2. Multivariate data analysis"
246 write(iout,*) " Fionn Murtagh and Andre Heck"
247 write(iout,*) " Kluwer Academic Publishers, 1987"
248 write(iout,*) "-------------------------------------------"
252 write (iout,*) "The TOTFREE array"
254 write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
258 CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
260 write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
262 write (iout,*) "Too few conformations to cluster."
265 CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
266 c CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
267 c 3/3/16 AL: added explicit number of cluters
268 if (nclust.gt.0) then
283 licz(iclass(j,i))=licz(iclass(j,i))+1
284 nconf(iclass(j,i),licz(iclass(j,i)))=j
285 c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
286 c & nconf(iclass(j,i),licz(iclass(j,i)))
292 IF (HEIGHT(L).EQ.IDUM) GOTO 190
295 write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
296 & " icut",icut," cutoff",rcutoff(icut)
297 IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
299 & WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
300 write (iout,'(a,f8.2)') 'Maximum distance found:',
302 CALL SRTCLUST(ICUT,ncon_work,iT)
304 CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
306 if (icut.gt.ncut) goto 191
313 licz(iclass(j,i))=licz(iclass(j,i))+1
314 nconf(iclass(j,i),licz(iclass(j,i)))=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,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,kkk)
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,kkk)
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,kkk),nend-nstart+1,
432 print *,'error, rms^2 = ',rms,icon,jcon
435 if (non_conv) print *,non_conv,icon,jcon
436 if (rmsmina.gt.rms) rmsmina=rms
438 difconf=dsqrt(rmsmina)
441 C------------------------------------------------------------------------------
442 subroutine distout(ncon)
445 include 'sizesclu.dat'
448 include 'COMMON.IOUNITS'
449 include 'COMMON.CLUSTER'
450 integer i,j,k,jlim,jlim1,nlim,ind,ioffset
453 write (iout,'(a)') 'The distance matrix'
455 nlim=min0(i+ncol-1,ncon)
456 write (iout,1000) (k,k=i,nlim)
457 write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
458 1000 format (/8x,10(i4,3x))
459 1020 format (/1x,80(1h-)/)
470 IND=IOFFSET(NCON,j,k)
472 IND=IOFFSET(NCON,k,j)
475 write (iout,1010) j,(b(k),k=1,jlim-i+1)
478 1010 format (i5,3x,10(f6.2,1x))