Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / cluster / cluster.F90
index 3b75b28..a47f418 100644 (file)
@@ -14,7 +14,7 @@
                          c,cref
       use energy_data, only: nnt,nct
       use control_data, only: symetr,outpdb,outmol2,titel,&
-                          iopt,print_dist !,MaxProcs
+                          iopt,print_dist,nclust !,MaxProcs
       use control, only: tcpu,initialize
 
       use wham_data, only: punch_dist
@@ -54,7 +54,7 @@
       INTEGER,dimension(:),allocatable :: HVALS !(maxconf-1)
       INTEGER,dimension(:),allocatable :: IORDER,HEIGHT !(maxconf-1)
       integer,dimension(:),allocatable :: nn !(maxconf)
-      integer :: ndis
+      integer :: ndis,is,ie
       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,&
@@ -96,8 +96,9 @@
 !elwrite(iout,*) "before parmread"
       call openunits
 !elwrite(iout,*) "before parmread"
-      call parmread
       call read_control
+      call parmread
+!      call read_control
 !elwrite(iout,*) "after read control"
       call molread
 !      if (refstr) call read_ref_structure(*30)
 !      write (iout,*) "after permut"
 !      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
           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
 
       call probabl(iT,ncon_work,ncon,*20)
 
-!elwrite(iout,*)"after probabl, ncon_work=", ncon_work,ncon
+       write(iout,*)"after probabl, ncon_work=", ncon_work,ncon
 
       if (ncon_work.lt.2) then
         write (iout,*) "Too few conformations; clustering skipped"
       call work_partition(.true.,ndis)
 #endif
 !el      call alloc_clust_arrays(ncon_work)
+      write(iout,*) "before icc allocate",(allocated(ICC)),ncon_work
+      if (allocated(ICC)) then
+       deallocate(ICC)
+       deallocate(DISS)
+      endif
       allocate(ICC(ncon_work))
       allocate(DISS(maxdist))
-
+      write(iout,*) "after icc allocate",(allocated(ICC)),ncon_work
       DO I=1,NCON_work
         ICC(I)=I
       ENDDO
 ! CALCULATE DISTANCES
 !
       call daread_ccoords(1,ncon_work)
+      write(iout,*) "after daread_ccoords"
       ind1=0
       DO I=1,NCON_work-1
         if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
        'Time for distance calculation:',T2-T1,' sec.'
       t1=tcpu()
       PRINT '(a)','End of distance computation'
+      if (allocated(diss_)) then
+       deallocate(diss_)
+       deallocate(scount_)
+      endif
 !el---------------
       allocate(diss_(maxdist))
       allocate(scount_(0:nprocs))
 !
 ! Print out the RMS deviation matrix.
 !
+      write(iout,*) "before distout"
       if (print_dist) CALL DISTOUT(NCON_work)
+      write(iout,*) "after distout"
 !
 !  call hierarchical clustering HC from F. Murtagh
 !
         write (iout,'(2i5,f10.5)') i,list_conf(i),totfree(i)
       enddo
 #endif
+      write (iout,*) "before CRIT", allocated(CRIT)
+      if (allocated(CRIT)) then
+      deallocate(CRIT)
+      deallocate(MEMBR)
+      deallocate(CRITVAL)
+      deallocate(IA)
+      deallocate(IB)
+      deallocate(IORDER)
+      deallocate(HEIGHT)
+      deallocate(ICLASS)
+      deallocate(HVALS)
+      deallocate(nn)
+      deallocate(DISNN)
+      deallocate(FLAG)
+      endif
       allocate(CRIT(N),MEMBR(N)) !(maxconf)
       allocate(CRITVAL(N-1)) !(maxconf-1)
       allocate(IA(N),IB(N))
       allocate(nn(N)) !(maxconf)
       allocate(DISNN(N)) !(maxconf)
       allocate(FLAG(N)) !(maxconf)
+      write (iout,*) "after CRIT", allocated(CRIT)
       call flush(iout)
       CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
       LEV = N-1
       endif
       CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
 !      CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
-
+      if (allocated(licz)) then
+      deallocate(licz)
+      deallocate(iass)
+      deallocate(nconf)
+      deallocate(totfree_gr)
+      endif
       allocate(licz(maxgr))
       allocate(iass(maxgr))
       allocate(nconf(maxgr,maxingr))
       allocate(totfree_gr(maxgr))
-
+!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
 !        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
-
+!      do i=1,lev-1
+      do i=is,ie
          idum=lev-i
          DO L=1,LEV
             IF (HEIGHT(L).EQ.IDUM) GOTO 190
  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
+         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)
          do l=1,maxgr
           licz(l)=0
          enddo
+         ii=i-is+1
          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
+          licz(iclass(j,ii))=licz(iclass(j,ii))+1
+          nconf(iclass(j,ii),licz(iclass(j,ii)))=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),
 !
 ! Compute free energies of clusters
 !
+      if (allocated(prob)) deallocate(prob)
       allocate(prob(maxgr))
 
       do igr=1,ngr