Adam's corrections to HOMOL_FRAG program

diff --git a/source/fragments/klapaucjusz-longest.f b/source/fragments/klapaucjusz-longest.f
index 41f49a3..8aa3f81 100644 (file)
--- a/source/fragments/klapaucjusz-longest.f
+++ b/source/fragments/klapaucjusz-longest.f
@@ -1,5 +1,6 @@
       parameter (maxmodel=20) 
       parameter (maxres=1000)
+      parameter (maxclust=1000)
       character*32 model(maxmodel)
       integer imodel(maxmodel)
       integer ifragpair(maxres,maxmodel,maxmodel),
@@ -7,14 +8,47 @@
      & iflag(maxres),isumfmax(maxmodel),kmax(maxmodel),
      & isumf(maxmodel),isumfmaxmin(maxmodel),kmaxmin(maxmodel),
      & jmaxmin(maxmodel),isump(maxmodel),ifragpair_temp(maxres),
-     & ifragpair_temp1(maxres),ninclust(maxmodel),nlenclust(maxmodel),
-     & icluster(maxmodel,maxmodel),ifrag(maxres,maxmodel)
+     & ifragpair_temp1(maxres),ninclust(maxclust),nlenclust(maxclust),
+     & icluster(maxclust,maxclust),ifrag(maxres,maxclust),
+     & imodelclust(maxclust)
       logical overlap
+      character*1 charm
+      character*8 ctemp
+      character*128 outfile
  
-      nmodel=5
-      nres=114
+      narg = iargc()
+      if (narg.lt.2) then
+        print '(2a)', "Usage: klapaucjusz nres nmodel",
+     &   " [lencut] [minclust_size] [outfile] < infile > logfile"
+        stop
+      endif
+      call getarg(1,ctemp)
+      read(ctemp,*) nres
+      call getarg(2,ctemp)
+      read(ctemp,*) nmodel
+      if (narg.gt.2) then
+        call getarg(3,ctemp)
+        read(ctemp,*) lencut
+      else
+        lencut=20
+      endif
+      if (narg.gt.3) then
+        call getarg(4,ctemp)
+        read(ctemp,*) minclust_size
+      else
+        minclust_size = 5
+      endif
+      if (narg.gt.4) then
+        call getarg(5,outfile)
+      else
+        outfile = "klapaucjusz.out"
+      endif
+      open(2,file=outfile,status="unknown")
+      print *,"nres",nres," nmodel",nmodel
+      print *,"lencut",lencut," minclust_size",minclust_size
+      print *," outfile ",outfile
       ifragpair_new=0 
-      lencut=7
+      nclust = 0
 
       do i=1,nmodel
         do j=1,nmodel
@@ -67,12 +101,27 @@ c      icut = max0(isummax*2/3,lencut)
       nflag=2
       iflag(imodel(1))=i
       iflag(imodel(2))=j
-c      write(*,*)"Maximum"
-c      write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
-      overlap=.true.
+      write(*,*)"Maximum"
+      write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
       i1 = imodel(1)
       i2 = imodel(2)
-      isumcut=0
+      isumcut=isummax
+      if (isummax.ge.lencut) then
+          kkmod=kfmax
+          kmaxminmax = kfmax
+          ifragpair_temp=0
+          ifragpair_temp1=0
+          do k=1,nres
+            if (ifragpair(k,i1,i2).eq.kfmax) then
+              ifragpair_temp(k)=ifragpair(k,i1,i2)
+              ifragpair_temp1(k)=ifragpair(k,i1,i2)
+            endif
+          enddo
+          overlap=.true.
+      else
+        overlap=.false.
+      endif
+      write (*,*) ">>>>>> dowhile"
       do while (overlap)
         do i=1,nmodel
           if (iflag(i).gt.0) cycle
@@ -81,7 +130,7 @@ c      write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
             isumf=0
             do k=1,nres
               kk = ifragpair(k,i,jj)
-              if (ifragpair(k,i1,i2).ne.kfmax .or. kk.eq.0) cycle
+              if (ifragpair_temp1(k).ne.kfmax .or. kk.eq.0) cycle
               isumf(kk)=isumf(kk)+1
             enddo 
             isumfmax(j)=0
@@ -127,41 +176,62 @@ c      write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
           nflag=nflag+1
           iflag(imaxminmax)=imaxminmax
           imodel(nflag)=imaxminmax 
-c          write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
-c          write(*,'(1000i1)')(ifragpair(k,imodel(1),imodel(2)),k=1,nres)
+          write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
+          write(*,'(1000i1)')(ifragpair(k,imodel(1),imodel(2)),k=1,nres)
           write(*,*) "isumfmaxminmax",isumfmaxminmax
           print *,"kmaxminmax",kmaxminmax
           write(*,'(a32)')(model(imodel(i)),i=1,nflag)
           write(*,'(1000i1)')
      &     (ifragpair(k,imaxminmax,jmaxminmax),k=1,nres)
+          ifragpair_temp=0
+          do k=1,nres
+            if (ifragpair(k,imaxminmax,jmaxminmax).eq.kmaxminmax
+     &       .and. ifragpair_temp1(k).gt.0)
+     &      ifragpair_temp(k)=ifragpair(k,imaxminmax,jmaxminmax)
+          enddo
           do k=1,nres
-            ifragpair_temp(k)=ifragpair(k,imaxminmax,jmaxminmax)
-            ifragpair_temp1(k)=ifragpair(k,i1,i2)
+            ifragpair_temp1(k)=ifragpair_temp(k)
           enddo
+         print *,"ifragpair_temp,ifragpair_temp1"
+         print '(1000i1)',(ifragpair_temp(k),k=1,nres)
+         print '(1000i1)',(ifragpair_temp1(k),k=1,nres)
+         do i=1,nmodel
+           do j=i+1,nmodel
+             write(*,'(2a32,1000i1)')model(i),model(j),
+     &       (ifragpair(k,i,j),k=1,nres)
+           enddo
+         enddo
+        endif
+      enddo ! while
+      print *,"iimod",iimod," jjmod",jjmod," kkmod",kkmod,
+     &   " isumcut",isumcut
+      if (isumcut.ge.icut) then
+          print *,"ifragpair_temp1 before zeroing"
+          print '(64x,1000i1)',(ifragpair_temp1(k),k=1,nres)
           do i=1,nflag
             do j=1,i-1
               ii = imodel(i)
               jj = imodel(j)
               do k=1,nres
-                if (ifragpair_temp(k).eq.kmaxminmax .and. 
-     &           ifragpair_temp1(k).eq.kfmax) then
+c                if (ifragpair_temp(k).eq.kmaxminmax .and. 
+c     &           ifragpair_temp1(k).eq.kfmax) then
+                if (ifragpair_temp1(k).gt.0) then 
                   ifragpair(k,ii,jj)=0
                   ifragpair(k,jj,ii)=0
                 endif
               enddo
             enddo
           enddo
+         print *,">>>>>>>>>>>>>>> IFRAGPAIR AFTER ZEROING <<<<<<<<<<<"
          do i=1,nmodel
            do j=i+1,nmodel
              write(*,'(2a32,1000i1)')model(i),model(j),
      &       (ifragpair(k,i,j),k=1,nres)
            enddo
          enddo
-        endif
-      enddo ! while
-      print *,"iimod",iimod," jjmod",jjmod," kkmod",kkmod,
-     &   " isumcut",isumcut
-      if (isumcut.ge.icut) then
+
+        if (nflag.ge.minclust_size) then
+
         nclust=nclust+1
         print *,"cluster",nclust
         ninclust(nclust)=nflag
@@ -171,10 +241,11 @@ c          write(*,'(1000i1)')(ifragpair(k,imodel(1),imodel(2)),k=1,nres)
         enddo
         ii=0
         do k=1,nres
-          if (ifragpair_temp(k).eq.kkmod .and. 
-     &     ifragpair_temp1(k).eq.kfmax) then
+c          if (ifragpair_temp(k).eq.kkmod .and. 
+c     &     ifragpair_temp1(k).eq.kfmax) then
+          if (ifragpair_temp1(k).gt.0) then
             ii=ii+1
-            print *,k
+c            print *,k
             ifrag(ii,nclust)=k
           endif
         enddo
@@ -183,26 +254,68 @@ c          write(*,'(1000i1)')(ifragpair(k,imodel(1),imodel(2)),k=1,nres)
         do i=1,nflag
           do j=1,i-1
             do k=1,nres
-              if (ifragpair_temp(k).eq.kkmod .and.
-     &           ifragpair_temp1(k).eq.kfmax) then
+              if (ifragpair_temp1(k).gt.0) then
                 ifragpair_new(k,imodel(i),imodel(j))=nclust
                 ifragpair_new(k,imodel(j),imodel(i))=nclust
               endif
             enddo
           enddo
         enddo
+
+        endif ! nflag > 2
+
       endif
       ENDDO ! while
       do i=1,nmodel
         do j=i+1,nmodel
-          write(*,'(2a32,1000i1)')model(i),model(j),
-     &      (ifragpair_new(k,i,j),k=1,nres)
+          write(*,'(2a32,1000a1)')model(i),model(j),
+     &      (charm(ifragpair_new(k,i,j)),k=1,nres)
         enddo
       enddo
-      write (*,*) "nclust",nclust 
+      write (*,*) "nclust",nclust
       do i=1,nclust
-        write(*,*) "Cluster",i," ninclust",ninclust(i)
+        write(*,*) "Cluster",i," ninclust",ninclust(i),
+     &   " length",nlenclust(i)
         write(*,'(a32)') (model(icluster(j,i)),j=1,ninclust(i))
         write(*,'(i5)') (ifrag(j,i),j=1,nlenclust(i))
       enddo
+      imodelclust=0
+      do i=1,nclust
+        do j=1,ninclust(i)
+          imodelclust(icluster(j,i))=1
+        enddo
+      enddo
+      print *,"Models in clusters"
+      nmodel_out = 0
+      do i=1,nmodel
+        if (imodelclust(i).eq.1) then
+          print '(a32)',model(i)
+          nmodel_out = nmodel_out + 1
+        endif
+      enddo
+      print *,"Models excluded"
+      do i=1,nmodel
+        if (imodelclust(i).eq.0) print '(a32)',model(i)
+      enddo
+c Write model information
+      write(2,'(i5)') nmodel_out,nclust
+      do i=1,nmodel
+        if (imodelclust(i).eq.1) write(2,'(a32)') model(i)
+      enddo
+c Write cluster information
+      do i=1,nclust
+        write(2,'(2i5)') ninclust(i),nlenclust(i)
+        write(2,'(16i5)') (icluster(k,i),k=1,ninclust(i))
+        write(2,'(16i5)') (ifrag(k,i),k=1,nlenclust(i))
+      enddo
+      end
+c---------------------------------------------------------------------------------
+      character function charm(i)
+      integer i
+      if (i.eq.0) then
+        charm = "."
+      else
+        charm = char(ichar("A")+i-1)
+      endif
+      return
       end