cluster average over multiple temperatures
[unres.git] / source / cluster / wham / src-M / main_clust.F
index f01f859..5b3cfe8 100644 (file)
@@ -23,23 +23,25 @@ C
       logical printang(max_cut)
       integer printpdb(max_cut)
       integer printmol2(max_cut)
-      character*240 lineh
+      character*240 lineh,scrachdir2d
       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
+      integer nn,ndis,scount_buf
+      real*4 DISNN, diss_buf(maxdist)
       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
+     & it,ncon_work,ind1,kkk, ijk,is,ie,ilen
       double precision t1,t2,tcpu,difconf
       
       double precision varia(maxvar)
       double precision hrtime,mintime,sectime
       logical eof
+      external ilen
+      integer nTend
 #ifdef MPI
       call MPI_Init( IERROR )
       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
@@ -61,6 +63,7 @@ C
 
       call initialize
       call openunits
+      call cinfo
       call parmread
       call read_control
       call molread
@@ -79,6 +82,12 @@ c      write (iout,*) "after permut"
 c      call flush(iout)
       print *,'MAIN: nnt=',nnt,' nct=',nct
 
+      if (nclust.gt.0) then
+        PRINTANG(1)=.TRUE.
+        PRINTPDB(1)=outpdb
+        printmol2(1)=outmol2
+        ncut=0
+      else
       DO I=1,NCUT
         PRINTANG(I)=.FALSE.
         PRINTPDB(I)=0
@@ -90,12 +99,21 @@ c      call flush(iout)
           printmol2(i)=outmol2
         ENDIF
       ENDDO
+      endif
+      if (ncut.gt.0) then
       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
+      else if (nclust.gt.0) then
+      write (iout,'("Number of clusters requested",i5)') nclust
+      else
+      if (me.eq.Master)
+     & write (iout,*) "ERROR: Either nclust or ncut must be >0"
+      stop
+      endif
       DO I=1,NRES-3  
         MULT(I)=1
       ENDDO
@@ -106,7 +124,12 @@ c      call flush(iout)
       write (iout,*) 'from read_coords: ncon',ncon
       
       write (iout,*) "nT",nT
-      do iT=1,nT
+      if (cumul_prob) then
+        nTend=1
+      else
+        nTend=nT
+      endif
+      do iT=1,nTend
       write (iout,*) "iT",iT
 #ifdef MPI
       call work_partition(.true.,ncon)
@@ -122,7 +145,7 @@ c      call flush(iout)
       ndis=ncon_work*(ncon_work-1)/2
       call work_partition(.true.,ndis)
 #endif
-
+      write(iout,*) "AFTET wort_part",NCON_work
       DO I=1,NCON_work
         ICC(I)=I
       ENDDO
@@ -140,9 +163,10 @@ C
             c(l,k)=allcart(l,k,i)
           enddo 
         enddo
+        kkk=1
         do k=1,nres
           do l=1,3
-            cref(l,k)=c(l,k)
+            cref(l,k,kkk)=c(l,k)
           enddo
         enddo
         DO J=I+1,NCON_work
@@ -164,12 +188,21 @@ c          write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND)
       t1=tcpu()
       PRINT '(a)','End of distance computation'
 
+      scount_buf=scount(me)
+
+      do ijk=1, ndis
+      diss_buf(ijk)=diss(ijk)
+      enddo
+
+
 #ifdef MPI
-      call MPI_Gatherv(diss(1),scount(me),MPI_REAL,diss(1),
+c      WRITE (iout,*) "Wchodze do call MPI_Gatherv"
+      call MPI_Gatherv(diss_buf(1),scount_buf,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')
+      scrachdir2d=scratchdir(:ilen(scratchdir))//'distance'
+      open(80,file=scrachdir2d,form='unformatted')
       do i=1,ndis
         write(80) diss(i)
       enddo
@@ -242,29 +275,39 @@ C
       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 c      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
 
+c 3/3/16 AL: added explicit number of cluters
+      if (nclust.gt.0) then
+        is=nclust-1
+        ie=nclust-1
+        icut=1
+      else
+        is=1
+        ie=lev-1
+      endif
       do i=1,maxgr
         licz(i)=0
       enddo
       icut=1
-      i=1
-      NGR=i+1
+      i=is
+      NGR=is+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
-
+c      do i=1,lev-1
+      do i=is,ie
          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)
+c         write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM),
+c     &    " icut",icut," cutoff",rcutoff(icut)
+         IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
+         if (nclust.le.0) 
+     &    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)
@@ -277,9 +320,10 @@ c     &    nconf(iclass(j,i),licz(iclass(j,i)))
          do l=1,maxgr
           licz(l)=0
          enddo
+         ii=i-is+1
          do j=1,n
-          licz(iclass(j,i))=licz(iclass(j,i))+1
-          nconf(iclass(j,i),licz(iclass(j,i)))=j
+          licz(iclass(j,ii))=licz(iclass(j,ii))+1
+          nconf(iclass(j,ii),licz(iclass(j,ii)))=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),
@@ -305,17 +349,17 @@ C
 C
       close(icbase,status="delete")
 #ifdef MPI
-      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
 #endif
       stop '********** Program terminated normally.'
    20 write (iout,*) "Error reading coordinates"
 #ifdef MPI
-      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
 #endif
       stop
    30 write (iout,*) "Error reading reference structure"
 #ifdef MPI
-      call MPI_Finalize(MPI_COMM_WORLD,IERROR)
+      call MPI_Finalize(IERROR)
 #endif
       stop
       end
@@ -358,7 +402,7 @@ c          write (iout,*) "tutaj",zzz
           ibezperm=(run-1)*chalen+i
           do j=1,3
             xx(j,ii)=allcart(j,iaperm,jcon)
-            yy(j,ii)=cref(j,ibezperm)
+            yy(j,ii)=cref(j,ibezperm,kkk)
           enddo
          enddo
         enddo
@@ -372,7 +416,7 @@ 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)
+              yy(j,ii)=cref(j,ibezperm+nres,kkk)
             enddo
            enddo
 c          endif
@@ -392,7 +436,8 @@ c        do i=nnt,nct
           enddo
          enddo
         enddo
-        call fitsq(rms,c(1,nstart),cref(1,nstart),nend-nstart+1,przes,
+        call fitsq(rms,c(1,nstart),cref(1,nstart,kkk),nend-nstart+1,
+     &       przes,
      &       obrot,non_conv)
       endif
       if (rms.lt.0.0) then