program cluster ! ! Program to cluster united-residue MCM results. ! use clust_data use probability use tracking use hc_ use io_clust !#define CLUSTER use io_units use io_base, only: permut use geometry_data, only: nres,theta,phi,alph,omeg,& c,cref use energy_data, only: nnt,nct use control_data, only: symetr,outpdb,outmol2,titel,& iopt,print_dist,MaxProcs use control, only: tcpu,initialize use wham_data, only: punch_dist use io_wham, only: parmread use work_part ! include 'DIMENSIONS' ! include 'sizesclu.dat' #ifdef MPI use mpi_data implicit none include "mpif.h" integer :: IERROR,ERRCODE !STATUS(MPI_STATUS_SIZE) #else implicit none ! include "COMMON.MPI" #endif ! include 'COMMON.TIME1' ! include 'COMMON.INTERACT' ! include 'COMMON.NAMES' ! include 'COMMON.GEO' ! include 'COMMON.HEADER' ! include 'COMMON.CONTROL' ! include 'COMMON.CHAIN' ! include 'COMMON.VAR' ! include 'COMMON.CLUSTER' ! include 'COMMON.IOUNITS' ! include 'COMMON.FREE' logical,dimension(:),allocatable :: printang !(max_cut) integer,dimension(:),allocatable :: printpdb !(max_cut) integer,dimension(:),allocatable :: printmol2 !(max_cut) character(len=240) lineh REAL(kind=4),dimension(:),allocatable :: CRIT,MEMBR !(maxconf) REAL(kind=4),dimension(:),allocatable :: CRITVAL !(maxconf-1) INTEGER,dimension(:),allocatable :: IA,IB !(maxconf) INTEGER,dimension(:,:),allocatable :: ICLASS !(maxconf,maxconf-1) INTEGER,dimension(:),allocatable :: HVALS !(maxconf-1) INTEGER,dimension(:),allocatable :: IORDER,HEIGHT !(maxconf-1) integer,dimension(:),allocatable :: nn !(maxconf) integer :: ndis real(kind=4),dimension(:),allocatable :: DISNN !(maxconf) LOGICAL,dimension(:),allocatable :: FLAG !(maxconf) integer :: i,j,k,l,m,n,len,lev,idum,ii,ind,jj,icut,ncon,& it,ncon_work,ind1,kkk real(kind=8) :: t1,t2,difconf real(kind=8),dimension(:),allocatable :: varia !(maxvar) real(kind=8),dimension(:),allocatable :: list_conf_ !(maxvar) real(kind=8) :: hrtime,mintime,sectime logical :: eof external :: difconf !el real(kind=4),dimension(:),allocatable :: diss_ !(maxdist) integer,dimension(:),allocatable :: scount_ !(maxdist) #ifdef MPI call MPI_Init( IERROR ) call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR ) call MPI_Comm_size( MPI_COMM_WORLD, nprocs, IERROR ) Master = 0 if (ierror.gt.0) then write(iout,*) "SEVERE ERROR - Can't initialize MPI." call mpi_finalize(ierror) stop endif if (nprocs.gt.MaxProcs+1) then write (2,*) "Error - too many processors",& nprocs,MaxProcs+1 write (2,*) "Increase MaxProcs and recompile" call MPI_Finalize(IERROR) stop endif #endif !elwrite(iout,*) "before parmread" allocate(printang(max_cut)) allocate(printpdb(max_cut)) allocate(printmol2(max_cut)) call initialize !elwrite(iout,*) "before parmread" call openunits !elwrite(iout,*) "before parmread" call parmread call read_control !elwrite(iout,*) "after read control" call molread ! if (refstr) call read_ref_structure(*30) do i=1,nres phi(i)=0.0D0 theta(i)=0.0D0 alph(i)=0.0D0 omeg(i)=0.0D0 enddo ! write (iout,*) "Before permut" ! write (iout,*) "symetr", symetr ! call flush(iout) call permut(symetr) ! write (iout,*) "after permut" ! call flush(iout) print *,'MAIN: nnt=',nnt,' nct=',nct DO I=1,NCUT PRINTANG(I)=.FALSE. PRINTPDB(I)=0 printmol2(i)=0 IF (RCUTOFF(I).LT.0.0) THEN RCUTOFF(I)=ABS(RCUTOFF(I)) PRINTANG(I)=.TRUE. PRINTPDB(I)=outpdb printmol2(i)=outmol2 ENDIF ENDDO write (iout,*) 'Number of cutoffs:',NCUT write (iout,*) 'Cutoff values:' DO ICUT=1,NCUT WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT),& printpdb(icut),printmol2(icut) ENDDO DO I=1,NRES-3 MULT(I)=1 ENDDO allocate(list_conf(maxconf)) do i=1,maxconf list_conf(i)=i enddo call read_coords(ncon,*20) allocate(list_conf_(maxconf)) do i=1,maxconf list_conf_(i)=list_conf(i) enddo deallocate(list_conf) allocate(list_conf(ncon)) do i=1,ncon list_conf(i)=list_conf_(i) enddo deallocate(list_conf_) !el call alloc_clust_arrays(ncon) write (iout,*) 'from read_coords: ncon',ncon write (iout,*) "nT",nT do iT=1,nT write (iout,*) "iT",iT #ifdef MPI call work_partition(.true.,ncon) #endif !elwrite(iout,*)"after work partition, ncon_work=", ncon_work,ncon call probabl(iT,ncon_work,ncon,*20) !elwrite(iout,*)"after probabl, ncon_work=", ncon_work,ncon if (ncon_work.lt.2) then write (iout,*) "Too few conformations; clustering skipped" exit endif #ifdef MPI ndis=ncon_work*(ncon_work-1)/2 call work_partition(.true.,ndis) #endif !el call alloc_clust_arrays(ncon_work) allocate(ICC(ncon_work)) allocate(DISS(maxdist)) DO I=1,NCON_work ICC(I)=I ENDDO WRITE (iout,'(A80)') TITEL t1=tcpu() ! ! CALCULATE DISTANCES ! call daread_ccoords(1,ncon_work) ind1=0 DO I=1,NCON_work-1 if (mod(i,100).eq.0) print *,'Calculating RMS i=',i do k=1,2*nres do l=1,3 c(l,k)=allcart(l,k,i) enddo enddo kkk=1 do k=1,nres do l=1,3 cref(l,k,kkk)=c(l,k) enddo enddo DO J=I+1,NCON_work IND=IOFFSET(NCON_work,I,J) #ifdef MPI if (ind.ge.indstart(me) .and. ind.le.indend(me)) then #endif ind1=ind1+1 DISS(IND1)=DIFCONF(I,J) ! write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND) #ifdef MPI endif #endif ENDDO ENDDO t2=tcpu() WRITE (iout,'(/a,1pe14.5,a/)') & 'Time for distance calculation:',T2-T1,' sec.' t1=tcpu() PRINT '(a)','End of distance computation' !el--------------- allocate(diss_(maxdist)) allocate(scount_(0:nprocs)) do i=1,maxdist diss_(i)=diss(i) enddo do i=0,nprocs scount_(i)=scount(i) enddo !el----------- #ifdef MPI call MPI_Gatherv(diss_(1),scount_(me),MPI_REAL,diss(1),& scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR) if (me.eq.master) then #endif deallocate(diss_) deallocate(scount_) open(80,file='/tmp/distance',form='unformatted') do i=1,ndis write(80) diss(i) enddo if (punch_dist) then do i=1,ncon_work-1 do j=i+1,ncon_work IND=IOFFSET(NCON,I,J) write (jrms,'(2i5,2f10.5)') i,j,diss(IND),& energy(j)-energy(i) enddo enddo endif ! ! Print out the RMS deviation matrix. ! if (print_dist) CALL DISTOUT(NCON_work) ! ! call hierarchical clustering HC from F. Murtagh ! N=NCON_work LEN = (N*(N-1))/2 write(iout,*) "-------------------------------------------" write(iout,*) "HIERARCHICAL CLUSTERING using" if (iopt.eq.1) then write(iout,*) "WARD'S MINIMUM VARIANCE METHOD" elseif (iopt.eq.2) then write(iout,*) "SINGLE LINK METHOD" elseif (iopt.eq.3) then write(iout,*) "COMPLETE LINK METHOD" elseif (iopt.eq.4) then write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD" elseif (iopt.eq.5) then write(iout,*) "MCQUITTY'S METHOD" elseif (iopt.eq.6) then write(iout,*) "MEDIAN (GOWER'S) METHOD" elseif (iopt.eq.7) then write(iout,*) "CENTROID METHOD" else write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7" write(*,*) "IOPT=",iopt," IS INVALID, use 1-7" stop endif write(iout,*) write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching" write(iout,*) "February 1986" write(iout,*) "References:" write(iout,*) "1. Multidimensional clustering algorithms" write(iout,*) " Fionn Murtagh" write(iout,*) " Vienna : Physica-Verlag, 1985." write(iout,*) "2. Multivariate data analysis" write(iout,*) " Fionn Murtagh and Andre Heck" write(iout,*) " Kluwer Academic Publishers, 1987" write(iout,*) "-------------------------------------------" write(iout,*) #ifdef DEBUG write (iout,*) "The TOTFREE array" do i=1,ncon_work write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i) enddo #endif allocate(CRIT(N),MEMBR(N)) !(maxconf) allocate(CRITVAL(N-1)) !(maxconf-1) allocate(IA(N),IB(N)) allocate(ICLASS(N,N-1)) !(maxconf,maxconf-1) allocate(HVALS(N-1)) !(maxconf-1) allocate(IORDER(N-1),HEIGHT(N-1)) !(maxconf-1) allocate(nn(N)) !(maxconf) allocate(DISNN(N)) !(maxconf) allocate(FLAG(N)) !(maxconf) call flush(iout) CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS) LEV = N-1 write (iout,*) "n",n," ncon_work",ncon_work," lev",lev if (lev.lt.2) then write (iout,*) "Too few conformations to cluster." goto 192 endif CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT) ! CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL) allocate(licz(maxgr)) allocate(iass(maxgr)) allocate(nconf(maxgr,maxingr)) allocate(totfree_gr(maxgr)) do i=1,maxgr licz(i)=0 enddo icut=1 i=1 NGR=i+1 do j=1,n licz(iclass(j,i))=licz(iclass(j,i))+1 nconf(iclass(j,i),licz(iclass(j,i)))=j ! write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)), ! & nconf(iclass(j,i),licz(iclass(j,i))) enddo do i=1,lev-1 idum=lev-i DO L=1,LEV IF (HEIGHT(L).EQ.IDUM) GOTO 190 ENDDO 190 IDUM=L write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),& " icut",icut," cutoff",rcutoff(icut) IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut) write (iout,'(a,f8.2)') 'Maximum distance found:',& CRITVAL(IDUM) CALL SRTCLUST(ICUT,ncon_work,iT) CALL TRACK(ICUT) CALL WRTCLUST(ncon_work,ICUT,PRINTANG,PRINTPDB,PRINTMOL2,iT) icut=icut+1 if (icut.gt.ncut) goto 191 ENDIF NGR=i+1 do l=1,maxgr licz(l)=0 enddo do j=1,n enddo do j=1,n licz(iclass(j,i))=licz(iclass(j,i))+1 nconf(iclass(j,i),licz(iclass(j,i)))=j !d write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),& !d nconf(iclass(j,i),licz(iclass(j,i))) !d print *,j,iclass(j,i), !d & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i))) enddo enddo 191 continue ! if (plot_tree) then CALL WRITRACK CALL PLOTREE endif ! t2=tcpu() WRITE (iout,'(/a,1pe14.5,a/)') & 'Total time for clustering:',T2-T1,' sec.' #ifdef MPI endif #endif 192 continue enddo ! close(icbase,status="delete") #ifdef MPI !el call MPI_Finalize(MPI_COMM_WORLD,IERROR) call MPI_Finalize(IERROR) #endif stop '********** Program terminated normally.' 20 write (iout,*) "Error reading coordinates" #ifdef MPI !el call MPI_Finalize(MPI_COMM_WORLD,IERROR) call MPI_Finalize(IERROR) #endif stop 30 write (iout,*) "Error reading reference structure" #ifdef MPI !el call MPI_Finalize(MPI_COMM_WORLD,IERROR) call MPI_Finalize(IERROR) #endif stop end program cluster !--------------------------------------------------------------------------- ! !--------------------------------------------------------------------------- real(kind=8) function difconf(icon,jcon) use clust_data use io_units, only: iout use io_base, only: permut use geometry_data, only: nres,c,cref,tabperm use energy_data, only: nct,nnt use control_data, only: symetr,lside,nend,nstart use regularize_, only: fitsq implicit none ! include 'DIMENSIONS' ! include 'sizesclu.dat' ! include 'COMMON.CONTROL' ! include 'COMMON.CLUSTER' ! include 'COMMON.CHAIN' ! include 'COMMON.INTERACT' ! include 'COMMON.VAR' ! include 'COMMON.IOUNITS' logical :: non_conv real(kind=8) :: przes(3),obrot(3,3) real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2) integer :: i,ii,j,icon,jcon,kkk,nperm,chalen,zzz integer :: iaperm,ibezperm,run real(kind=8) :: rms,rmsmina ! write (iout,*) "tu dochodze" rmsmina=10d10 nperm=1 do i=1,symetr nperm=i*nperm enddo ! write (iout,*) "nperm",nperm call permut(symetr) ! write (iout,*) "tabperm", tabperm(1,1) do kkk=1,nperm if (lside) then ii=0 chalen=int((nend-nstart+2)/symetr) do run=1,symetr do i=nstart,(nstart+chalen-1) zzz=tabperm(kkk,run) ! write (iout,*) "tutaj",zzz ii=ii+1 iaperm=(zzz-1)*chalen+i ibezperm=(run-1)*chalen+i do j=1,3 xx(j,ii)=allcart(j,iaperm,jcon) yy(j,ii)=cref(j,ibezperm,kkk) enddo enddo enddo do run=1,symetr do i=nstart,(nstart+chalen-1) zzz=tabperm(kkk,run) ii=ii+1 iaperm=(zzz-1)*chalen+i ibezperm=(run-1)*chalen+i ! if (itype(i).ne.10) then ii=ii+1 do j=1,3 xx(j,ii)=allcart(j,iaperm+nres,jcon) yy(j,ii)=cref(j,ibezperm+nres,kkk) enddo enddo ! endif enddo call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv) else chalen=int((nct-nnt+2)/symetr) do run=1,symetr do i=nnt,(nnt+chalen-1) zzz=tabperm(kkk,run) ! write (iout,*) "tu szukaj", zzz,run,kkk iaperm=(zzz-1)*chalen+i ibezperm=(run-1)*chalen+i ! do i=nnt,nct do j=1,3 c(j,i)=allcart(j,iaperm,jcon) enddo enddo enddo call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,& przes,& obrot,non_conv) endif if (rms.lt.0.0) then print *,'error, rms^2 = ',rms,icon,jcon stop endif if (non_conv) print *,non_conv,icon,jcon if (rmsmina.gt.rms) rmsmina=rms enddo difconf=dsqrt(rmsmina) return end function difconf !------------------------------------------------------------------------------ subroutine distout(ncon) use clust_data use hc_, only:ioffset use io_units, only: iout implicit none ! include 'DIMENSIONS' ! include 'sizesclu.dat' integer :: ncon integer,parameter :: ncol=10 ! include 'COMMON.IOUNITS' ! include 'COMMON.CLUSTER' integer :: i,j,k,jlim,jlim1,nlim,ind real(kind=4) :: b(ncol) write (iout,'(a)') 'The distance matrix' do 1 i=1,ncon,ncol nlim=min0(i+ncol-1,ncon) write (iout,1000) (k,k=i,nlim) write (iout,'(8h--------,10a)') ('-------',k=i,nlim) 1000 format (/8x,10(i4,3x)) 1020 format (/1x,80(1h-)/) do 2 j=i,ncon jlim=min0(j,nlim) if (jlim.eq.j) then b(jlim-i+1)=0.0d0 jlim1=jlim-1 else jlim1=jlim endif do 3 k=i,jlim1 if (j.lt.k) then IND=IOFFSET(NCON,j,k) else IND=IOFFSET(NCON,k,j) endif 3 b(k-i+1)=diss(IND) write (iout,1010) j,(b(k),k=1,jlim-i+1) 2 continue 1 continue 1010 format (i5,3x,10(f6.2,1x)) return end subroutine distout !------------------------------------------------------------------------------ ! srtclust.f !------------------------------------------------------------------------------ SUBROUTINE SRTCLUST(ICUT,NCON,IB) use clust_data use io_units, only: iout ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'sizesclu.dat' ! include 'COMMON.CLUSTER' ! include 'COMMON.FREE' ! include 'COMMON.IOUNITS' implicit none real(kind=8),dimension(:),allocatable :: prob !(maxgr) real(kind=8) :: emin,ene,en1,sumprob integer :: igr,i,ii,li1,li2,ligr,ico,jco,ind1,ind2 integer :: jgr,li,nco,ib,ncon,icut ! ! Compute free energies of clusters ! allocate(prob(maxgr)) do igr=1,ngr emin=totfree(nconf(igr,1)) totfree_gr(igr)=1.0d0 do i=2,licz(igr) ii=nconf(igr,i) totfree_gr(igr)=totfree_gr(igr)+dexp(-totfree(ii)+emin) enddo ! write (iout,*) "igr",igr," totfree",emin, ! & " totfree_gr",totfree_gr(igr) totfree_gr(igr)=emin-dlog(totfree_gr(igr)) ! write (iout,*) igr," efree",totfree_gr(igr)/beta_h(ib) enddo ! ! SORT CONFORMATIONS IN GROUPS ACC. TO ENERGY ! DO 16 IGR=1,NGR LIGR=LICZ(IGR) DO 17 ICO=1,LIGR-1 IND1=NCONF(IGR,ICO) ENE=totfree(IND1) DO 18 JCO=ICO+1,LIGR IND2=NCONF(IGR,JCO) EN1=totfree(IND2) IF (EN1.LT.ENE) THEN NCONF(IGR,ICO)=IND2 NCONF(IGR,JCO)=IND1 IND1=IND2 ENE=EN1 ENDIF 18 CONTINUE 17 CONTINUE 16 CONTINUE ! ! SORT GROUPS ! DO 71 IGR=1,NGR ENE=totfree_gr(IGR) DO 72 JGR=IGR+1,NGR EN1=totfree_gr(JGR) IF (EN1.LT.ENE) THEN LI1=LICZ(IGR) LI2=LICZ(JGR) LI=MAX0(LI1,LI2) DO 73 I=1,LI NCO=NCONF(IGR,I) NCONF(IGR,I)=NCONF(JGR,I) NCONF(JGR,I)=NCO 73 CONTINUE totfree_gr(igr)=en1 totfree_gr(jgr)=ene ENE=EN1 LICZ(IGR)=LI2 LICZ(JGR)=LI1 ENDIF 72 CONTINUE 71 CONTINUE write (iout,'("Free energies and probabilities of clusters at",f6.1," K")')& 1.0d0/(1.987d-3*beta_h(ib)) !' prob(1)=1.0d0 sumprob=1.0d0 do i=2,ngr prob(i)=dexp(-(totfree_gr(i)-totfree_gr(1))) sumprob=sumprob+prob(i) enddo do i=1,ngr prob(i)=prob(i)/sumprob enddo sumprob=0.0d0 write (iout,'("clust efree prob sumprob")') do i=1,ngr sumprob=sumprob+prob(i) write (iout,'(i5,f8.1,2f8.5)') i,totfree_gr(i)/beta_h(ib),& prob(i),sumprob enddo DO 81 IGR=1,NGR LI=LICZ(IGR) DO 82 I=1,LI 82 IASS(NCONF(IGR,I))=IGR 81 CONTINUE if (lgrp) then do i=1,ncon iass_tot(i,icut)=iass(i) ! write (iout,*) icut,i,iass(i),iass_tot(i,icut) enddo endif RETURN END SUBROUTINE SRTCLUST !------------------------------------------------------------------------------ !------------------------------------------------------------------------------