--- /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,kkk
+ 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
+c write (iout,*) "Before permut"
+c write (iout,*) "symetr", symetr
+c call flush(iout)
+ call permut(symetr)
+c write (iout,*) "after permut"
+c 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
+ 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
+ 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)
+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,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
+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,kkk)
+ 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,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
+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 / unres4.git / blobdiff