+++ /dev/null
- 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
-!------------------------------------------------------------------------------
-!------------------------------------------------------------------------------
projects / unres4.git / blobdiff