+++ /dev/null
-C
-C Program to cluster united-residue MCM results.
-C
- implicit none
- include 'DIMENSIONS'
- include 'sizesclu.dat'
-#ifdef MPI
- include "mpif.h"
- integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
- 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 printang(max_cut)
- integer printpdb(max_cut)
- integer printmol2(max_cut)
- character*240 lineh
- REAL CRIT(maxconf),MEMBR(maxconf)
- REAL CRITVAL(maxconf-1)
- INTEGER IA(maxconf),IB(maxconf)
- INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
- INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
- integer nn,ndis
- real*4 DISNN
- DIMENSION NN(maxconf),DISNN(maxconf)
- LOGICAL FLAG(maxconf)
- integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon,
- & it,ncon_work,ind1
- double precision t1,t2,tcpu,difconf
-
- double precision varia(maxvar)
- double precision hrtime,mintime,sectime
- logical eof
-#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
-
- call initialize
- call openunits
- call parmread
- call read_control
- call molread
-c 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,*) "symetr", symetr
- call permut(symetr)
- 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
- do i=1,maxconf
- list_conf(i)=i
- enddo
- call read_coords(ncon,*20)
- 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
-
- call probabl(iT,ncon_work,ncon,*20)
-
- 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
-
- DO I=1,NCON_work
- ICC(I)=I
- ENDDO
- WRITE (iout,'(A80)') TITEL
- t1=tcpu()
-C
-C CALCULATE DISTANCES
-C
- 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
- do k=1,nres
- do l=1,3
- cref(l,k)=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)
-c 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'
-
-#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
- 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
-C
-C Print out the RMS deviation matrix.
-C
- if (print_dist) CALL DISTOUT(NCON_work)
-C
-C call hierarchical clustering HC from F. Murtagh
-C
- 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
- 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)
-c CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
-
- 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
-c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
-c & 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
- licz(iclass(j,i))=licz(iclass(j,i))+1
- nconf(iclass(j,i),licz(iclass(j,i)))=j
-c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)),
-c & nconf(iclass(j,i),licz(iclass(j,i)))
-cd print *,j,iclass(j,i),
-cd & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
- enddo
- enddo
- 191 continue
-C
- if (plot_tree) then
- CALL WRITRACK
- CALL PLOTREE
- endif
-C
- t2=tcpu()
- WRITE (iout,'(/a,1pe14.5,a/)')
- & 'Total time for clustering:',T2-T1,' sec.'
-
-#ifdef MPI
- endif
-#endif
- 192 continue
- enddo
-C
- close(icbase,status="delete")
-#ifdef MPI
- call MPI_Finalize(MPI_COMM_WORLD,IERROR)
-#endif
- stop '********** Program terminated normally.'
- 20 write (iout,*) "Error reading coordinates"
-#ifdef MPI
- call MPI_Finalize(MPI_COMM_WORLD,IERROR)
-#endif
- stop
- 30 write (iout,*) "Error reading reference structure"
-#ifdef MPI
- call MPI_Finalize(MPI_COMM_WORLD,IERROR)
-#endif
- stop
- end
-c---------------------------------------------------------------------------
- double precision function difconf(icon,jcon)
- 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
- double precision przes(3),obrot(3,3)
- double precision xx(3,maxres2),yy(3,maxres2)
- integer i,ii,j,icon,jcon,kkk,nperm,chalen,zzz
- integer iaperm,ibezperm,run
- double precision rms,rmsmina
-c write (iout,*) "tu dochodze"
- rmsmina=10d10
- nperm=1
- do i=1,symetr
- nperm=i*nperm
- enddo
-c write (iout,*) "nperm",nperm
- call permut(symetr)
-c 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)
-c 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)
- 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
-c 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)
- enddo
- enddo
-c 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)
-c write (iout,*) "tu szukaj", zzz,run,kkk
- iaperm=(zzz-1)*chalen+i
- ibezperm=(run-1)*chalen+i
-c 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),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
-C------------------------------------------------------------------------------
- subroutine distout(ncon)
- implicit none
- include 'DIMENSIONS'
- include 'sizesclu.dat'
- integer ncol,ncon
- parameter (ncol=10)
- include 'COMMON.IOUNITS'
- include 'COMMON.CLUSTER'
- integer i,j,k,jlim,jlim1,nlim,ind,ioffset
- real*4 b
- dimension 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
projects / unres.git / blobdiff