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 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 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