3 ! Program to cluster united-residue MCM results.
12 use io_base, only: permut
13 use geometry_data, only: nres,theta,phi,alph,omeg,&
15 use energy_data, only: nnt,nct
16 use control_data, only: symetr,outpdb,outmol2,titel,&
17 iopt,print_dist,nclust !,MaxProcs
18 use control, only: tcpu,initialize
20 use wham_data, only: punch_dist
21 use io_wham, only: parmread
23 ! include 'DIMENSIONS'
24 ! include 'sizesclu.dat'
29 integer :: IERROR,ERRCODE !STATUS(MPI_STATUS_SIZE)
32 ! include "COMMON.MPI"
34 ! include 'COMMON.TIME1'
35 ! include 'COMMON.INTERACT'
36 ! include 'COMMON.NAMES'
37 ! include 'COMMON.GEO'
38 ! include 'COMMON.HEADER'
39 ! include 'COMMON.CONTROL'
40 ! include 'COMMON.CHAIN'
41 ! include 'COMMON.VAR'
42 ! include 'COMMON.CLUSTER'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.FREE'
46 logical,dimension(:),allocatable :: printang !(max_cut)
47 integer,dimension(:),allocatable :: printpdb !(max_cut)
48 integer,dimension(:),allocatable :: printmol2 !(max_cut)
49 character(len=240) lineh
50 REAL(kind=4),dimension(:),allocatable :: CRIT,MEMBR !(maxconf)
51 REAL(kind=4),dimension(:),allocatable :: CRITVAL !(maxconf-1)
52 INTEGER,dimension(:),allocatable :: IA,IB !(maxconf)
53 INTEGER,dimension(:,:),allocatable :: ICLASS !(maxconf,maxconf-1)
54 INTEGER,dimension(:),allocatable :: HVALS !(maxconf-1)
55 INTEGER,dimension(:),allocatable :: IORDER,HEIGHT !(maxconf-1)
56 integer,dimension(:),allocatable :: nn !(maxconf)
58 real(kind=4),dimension(:),allocatable :: DISNN !(maxconf)
59 LOGICAL,dimension(:),allocatable :: FLAG !(maxconf)
60 integer :: i,j,k,l,m,n,len,lev,idum,ii,ind,jj,icut,ncon,&
62 real(kind=8) :: t1,t2,difconf
64 real(kind=8),dimension(:),allocatable :: varia !(maxvar)
65 real(kind=8),dimension(:),allocatable :: list_conf_ !(maxvar)
66 real(kind=8) :: hrtime,mintime,sectime
70 real(kind=4),dimension(:),allocatable :: diss_ !(maxdist)
71 integer,dimension(:),allocatable :: scount_ !(maxdist)
73 call MPI_Init( IERROR )
74 call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
75 call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR )
78 write(iout,*) "SEVERE ERROR - Can't initialize MPI."
79 call mpi_finalize(ierror)
82 !EL /04/2017 No longer needed, dynamic allocation of arrays
83 ! if (nprocs.gt.MaxProcs+1) then
84 ! write (2,*) "Error - too many processors",&
86 ! write (2,*) "Increase MaxProcs and recompile"
87 ! call MPI_Finalize(IERROR)
91 !elwrite(iout,*) "before parmread"
92 allocate(printang(max_cut))
93 allocate(printpdb(max_cut))
94 allocate(printmol2(max_cut))
96 !elwrite(iout,*) "before parmread"
98 !elwrite(iout,*) "before parmread"
102 !elwrite(iout,*) "after read control"
104 ! if (refstr) call read_ref_structure(*30)
111 ! write (iout,*) "Before permut"
112 ! write (iout,*) "symetr", symetr
115 ! write (iout,*) "after permut"
117 print *,'MAIN: nnt=',nnt,' nct=',nct
118 if (nclust.gt.0) then
128 IF (RCUTOFF(I).LT.0.0) THEN
129 RCUTOFF(I)=ABS(RCUTOFF(I))
137 write (iout,*) 'Number of cutoffs:',NCUT
138 write (iout,*) 'Cutoff values:'
140 WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),&
141 printpdb(icut),printmol2(icut)
143 else if (nclust.gt.0) then
144 write (iout,'("Number of clusters requested",i5)') nclust
147 write (iout,*) "ERROR: Either nclust or ncut must be >0"
153 allocate(list_conf(maxconf))
157 call read_coords(ncon,*20)
159 allocate(list_conf_(maxconf))
161 list_conf_(i)=list_conf(i)
163 deallocate(list_conf)
164 allocate(list_conf(ncon))
166 list_conf(i)=list_conf_(i)
168 deallocate(list_conf_)
170 !el call alloc_clust_arrays(ncon)
172 write (iout,*) 'from read_coords: ncon',ncon
174 write (iout,*) "nT",nT
176 write (iout,*) "iT",iT
178 call work_partition(.true.,ncon)
180 !elwrite(iout,*)"after work partition, ncon_work=", ncon_work,ncon
182 call probabl(iT,ncon_work,ncon,*20)
184 write(iout,*)"after probabl, ncon_work=", ncon_work,ncon
186 if (ncon_work.lt.2) then
187 write (iout,*) "Too few conformations; clustering skipped"
191 ndis=ncon_work*(ncon_work-1)/2
192 call work_partition(.true.,ndis)
194 !el call alloc_clust_arrays(ncon_work)
195 write(iout,*) "before icc allocate",(allocated(ICC)),ncon_work
196 if (allocated(ICC)) then
200 allocate(ICC(ncon_work))
201 allocate(DISS(maxdist))
202 write(iout,*) "after icc allocate",(allocated(ICC)),ncon_work
206 WRITE (iout,'(A80)') TITEL
209 ! CALCULATE DISTANCES
211 call daread_ccoords(1,ncon_work)
212 write(iout,*) "after daread_ccoords"
215 if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
218 c(l,k)=allcart(l,k,i)
228 IND=IOFFSET(NCON_work,I,J)
230 if (ind.ge.indstart(me) .and. ind.le.indend(me)) then
233 DISS(IND1)=DIFCONF(I,J)
234 ! write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
241 WRITE (iout,'(/a,1pe14.5,a/)') &
242 'Time for distance calculation:',T2-T1,' sec.'
244 PRINT '(a)','End of distance computation'
245 if (allocated(diss_)) then
250 allocate(diss_(maxdist))
251 allocate(scount_(0:nprocs))
261 call MPI_Gatherv(diss_(1),scount_(me),MPI_REAL,diss(1),&
262 scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR)
263 if (me.eq.master) then
268 open(80,file='/tmp/distance',form='unformatted')
275 IND=IOFFSET(NCON,I,J)
276 write (jrms,'(2i5,2f10.5)') i,j,diss(IND),&
282 ! Print out the RMS deviation matrix.
284 write(iout,*) "before distout"
285 if (print_dist) CALL DISTOUT(NCON_work)
286 write(iout,*) "after distout"
288 ! call hierarchical clustering HC from F. Murtagh
292 write(iout,*) "-------------------------------------------"
293 write(iout,*) "HIERARCHICAL CLUSTERING using"
295 write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
296 elseif (iopt.eq.2) then
297 write(iout,*) "SINGLE LINK METHOD"
298 elseif (iopt.eq.3) then
299 write(iout,*) "COMPLETE LINK METHOD"
300 elseif (iopt.eq.4) then
301 write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
302 elseif (iopt.eq.5) then
303 write(iout,*) "MCQUITTY'S METHOD"
304 elseif (iopt.eq.6) then
305 write(iout,*) "MEDIAN (GOWER'S) METHOD"
306 elseif (iopt.eq.7) then
307 write(iout,*) "CENTROID METHOD"
309 write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
310 write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
314 write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
315 write(iout,*) "February 1986"
316 write(iout,*) "References:"
317 write(iout,*) "1. Multidimensional clustering algorithms"
318 write(iout,*) " Fionn Murtagh"
319 write(iout,*) " Vienna : Physica-Verlag, 1985."
320 write(iout,*) "2. Multivariate data analysis"
321 write(iout,*) " Fionn Murtagh and Andre Heck"
322 write(iout,*) " Kluwer Academic Publishers, 1987"
323 write(iout,*) "-------------------------------------------"
327 write (iout,*) "The TOTFREE array"
329 write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
332 write (iout,*) "before CRIT", allocated(CRIT)
333 if (allocated(CRIT)) then
347 allocate(CRIT(N),MEMBR(N)) !(maxconf)
348 allocate(CRITVAL(N-1)) !(maxconf-1)
349 allocate(IA(N),IB(N))
350 allocate(ICLASS(N,N-1)) !(maxconf,maxconf-1)
351 allocate(HVALS(N-1)) !(maxconf-1)
352 allocate(IORDER(N-1),HEIGHT(N-1)) !(maxconf-1)
353 allocate(nn(N)) !(maxconf)
354 allocate(DISNN(N)) !(maxconf)
355 allocate(FLAG(N)) !(maxconf)
356 write (iout,*) "after CRIT", allocated(CRIT)
358 CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
360 write (iout,*) "n",n," ncon_work",ncon_work," lev",lev
362 write (iout,*) "Too few conformations to cluster."
365 CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
366 ! CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
367 if (allocated(licz)) then
371 deallocate(totfree_gr)
373 allocate(licz(maxgr))
374 allocate(iass(maxgr))
375 allocate(nconf(maxgr,maxingr))
376 allocate(totfree_gr(maxgr))
377 !c 3/3/16 AL: added explicit number of cluters
378 if (nclust.gt.0) then
393 licz(iclass(j,i))=licz(iclass(j,i))+1
394 nconf(iclass(j,i),licz(iclass(j,i)))=j
395 ! write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
396 ! & nconf(iclass(j,i),licz(iclass(j,i)))
402 IF (HEIGHT(L).EQ.IDUM) GOTO 190
405 write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),&
406 " icut",icut," cutoff",rcutoff(icut)
407 IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
409 WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
410 write (iout,'(a,f8.2)') 'Maximum distance found:',&
412 CALL SRTCLUST(ICUT,ncon_work,iT)
414 CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT)
416 if (icut.gt.ncut) goto 191
424 licz(iclass(j,ii))=licz(iclass(j,ii))+1
425 nconf(iclass(j,ii),licz(iclass(j,ii)))=j
426 !d write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),&
427 !d nconf(iclass(j,i),licz(iclass(j,i)))
428 !d print *,j,iclass(j,i),
429 !d & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
440 WRITE (iout,'(/a,1pe14.5,a/)') &
441 'Total time for clustering:',T2-T1,' sec.'
448 close(icbase,status="delete")
450 !el call MPI_Finalize(MPI_COMM_WORLD,IERROR)
451 call MPI_Finalize(IERROR)
453 stop '********** Program terminated normally.'
454 20 write (iout,*) "Error reading coordinates"
456 !el call MPI_Finalize(MPI_COMM_WORLD,IERROR)
457 call MPI_Finalize(IERROR)
460 30 write (iout,*) "Error reading reference structure"
462 !el call MPI_Finalize(MPI_COMM_WORLD,IERROR)
463 call MPI_Finalize(IERROR)
467 !---------------------------------------------------------------------------
469 !---------------------------------------------------------------------------
470 real(kind=8) function difconf(icon,jcon)
474 use io_units, only: iout
475 use io_base, only: permut
476 use geometry_data, only: nres,c,cref,tabperm
477 use energy_data, only: nct,nnt
478 use control_data, only: symetr,lside,nend,nstart
479 use regularize_, only: fitsq
481 ! include 'DIMENSIONS'
482 ! include 'sizesclu.dat'
483 ! include 'COMMON.CONTROL'
484 ! include 'COMMON.CLUSTER'
485 ! include 'COMMON.CHAIN'
486 ! include 'COMMON.INTERACT'
487 ! include 'COMMON.VAR'
488 ! include 'COMMON.IOUNITS'
490 real(kind=8) :: przes(3),obrot(3,3)
491 real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
492 integer :: i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
493 integer :: iaperm,ibezperm,run
494 real(kind=8) :: rms,rmsmina
495 ! write (iout,*) "tu dochodze"
501 ! write (iout,*) "nperm",nperm
503 ! write (iout,*) "tabperm", tabperm(1,1)
507 chalen=int((nend-nstart+2)/symetr)
509 do i=nstart,(nstart+chalen-1)
511 ! write (iout,*) "tutaj",zzz
513 iaperm=(zzz-1)*chalen+i
514 ibezperm=(run-1)*chalen+i
516 xx(j,ii)=allcart(j,iaperm,jcon)
517 yy(j,ii)=cref(j,ibezperm,kkk)
522 do i=nstart,(nstart+chalen-1)
525 iaperm=(zzz-1)*chalen+i
526 ibezperm=(run-1)*chalen+i
527 ! if (itype(i).ne.10) then
530 xx(j,ii)=allcart(j,iaperm+nres,jcon)
531 yy(j,ii)=cref(j,ibezperm+nres,kkk)
536 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
538 chalen=int((nct-nnt+2)/symetr)
540 do i=nnt,(nnt+chalen-1)
542 ! write (iout,*) "tu szukaj", zzz,run,kkk
543 iaperm=(zzz-1)*chalen+i
544 ibezperm=(run-1)*chalen+i
547 c(j,i)=allcart(j,iaperm,jcon)
551 call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,&
556 print *,'error, rms^2 = ',rms,icon,jcon
559 if (non_conv) print *,non_conv,icon,jcon
560 if (rmsmina.gt.rms) rmsmina=rms
562 difconf=dsqrt(rmsmina)
565 !------------------------------------------------------------------------------
566 subroutine distout(ncon)
569 use hc_, only:ioffset
570 use io_units, only: iout
572 ! include 'DIMENSIONS'
573 ! include 'sizesclu.dat'
575 integer,parameter :: ncol=10
576 ! include 'COMMON.IOUNITS'
577 ! include 'COMMON.CLUSTER'
578 integer :: i,j,k,jlim,jlim1,nlim,ind
579 real(kind=4) :: b(ncol)
581 write (iout,'(a)') 'The distance matrix'
583 nlim=min0(i+ncol-1,ncon)
584 write (iout,1000) (k,k=i,nlim)
585 write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
586 1000 format (/8x,10(i4,3x))
587 1020 format (/1x,80(1h-)/)
598 IND=IOFFSET(NCON,j,k)
600 IND=IOFFSET(NCON,k,j)
603 write (iout,1010) j,(b(k),k=1,jlim-i+1)
606 1010 format (i5,3x,10(f6.2,1x))
608 end subroutine distout
609 !------------------------------------------------------------------------------
611 !------------------------------------------------------------------------------
612 SUBROUTINE SRTCLUST(ICUT,NCON,IB)
615 use io_units, only: iout
616 ! implicit real*8 (a-h,o-z)
617 ! include 'DIMENSIONS'
618 ! include 'sizesclu.dat'
619 ! include 'COMMON.CLUSTER'
620 ! include 'COMMON.FREE'
621 ! include 'COMMON.IOUNITS'
623 real(kind=8),dimension(:),allocatable :: prob !(maxgr)
624 real(kind=8) :: emin,ene,en1,sumprob
625 integer :: igr,i,ii,li1,li2,ligr,ico,jco,ind1,ind2
626 integer :: jgr,li,nco,ib,ncon,icut
628 ! Compute free energies of clusters
630 if (allocated(prob)) deallocate(prob)
631 allocate(prob(maxgr))
634 emin=totfree(nconf(igr,1))
635 totfree_gr(igr)=1.0d0
638 totfree_gr(igr)=totfree_gr(igr)+dexp(-totfree(ii)+emin)
640 ! write (iout,*) "igr",igr," totfree",emin,
641 ! & " totfree_gr",totfree_gr(igr)
642 totfree_gr(igr)=emin-dlog(totfree_gr(igr))
643 ! write (iout,*) igr," efree",totfree_gr(igr)/beta_h(ib)
646 ! SORT CONFORMATIONS IN GROUPS ACC. TO ENERGY
678 NCONF(IGR,I)=NCONF(JGR,I)
689 write (iout,'("Free energies and probabilities of clusters at",f6.1," K")')&
690 1.0d0/(1.987d-3*beta_h(ib)) !'
694 prob(i)=dexp(-(totfree_gr(i)-totfree_gr(1)))
695 sumprob=sumprob+prob(i)
698 prob(i)=prob(i)/sumprob
701 write (iout,'("clust efree prob sumprob")')
703 sumprob=sumprob+prob(i)
704 write (iout,'(i5,f8.1,2f8.5)') i,totfree_gr(i)/beta_h(ib),&
710 82 IASS(NCONF(IGR,I))=IGR
714 iass_tot(i,icut)=iass(i)
715 ! write (iout,*) icut,i,iass(i),iass_tot(i,icut)
719 END SUBROUTINE SRTCLUST
720 !------------------------------------------------------------------------------
721 !------------------------------------------------------------------------------