Adding MARTINI

[unres4.git] / source / cluster / cluster.F90

diff --git a/source/cluster/cluster.F90 b/source/cluster/cluster.F90
index 3b75b28..70b685e 100644 (file)
--- a/source/cluster/cluster.F90
+++ b/source/cluster/cluster.F90
@@ -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)
@@ -114,7 +115,12 @@
 !      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
@@ -126,12 +132,21 @@
           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
@@ -327,21 +342,29 @@
       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
@@ -349,7 +372,8 @@
  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)
@@ -363,11 +387,10 @@
          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),