source/fragments/klapaucjusz-longest.f

   1       parameter (maxmodel=50)
   2       parameter (maxres=1000)
   3       parameter (maxclust=1000)
   4       character*32 model(maxmodel)
   5       integer imodel(maxmodel),iindex(maxmodel)
   6       integer ifragpair(maxres,maxmodel,maxmodel),
   7      &    ifragpair_new(maxres,maxmodel,maxmodel),
   8      & iflag(maxres),isumfmax(maxmodel),kmax(maxmodel),
   9      & isumf(maxmodel),isumfmaxmin(maxmodel),kmaxmin(maxmodel),
  10      & jmaxmin(maxmodel),isump(maxmodel),ifragpair_temp(maxres),
  11      & ifragpair_temp1(maxres),ninclust(maxclust),nlenclust(maxclust),
  12      & icluster(maxclust,maxclust),ifrag(maxres,maxclust),
  13      & imodelclust(maxclust)
  14       logical overlap
  15       character*1 charm
  16       character*8 ctemp
  17       character*128 outfile
  18
  19       narg = iargc()
  20       if (narg.lt.2) then
  21         print '(2a)', "Usage: klapaucjusz nres nmodel",
  22      &   " [lencut] [minclust_size] [outfile] < infile > logfile"
  23         stop
  24       endif
  25       call getarg(1,ctemp)
  26       read(ctemp,*) nres
  27       call getarg(2,ctemp)
  28       read(ctemp,*) nmodel
  29       if (narg.gt.2) then
  30         call getarg(3,ctemp)
  31         read(ctemp,*) lencut
  32       else
  33         lencut=20
  34       endif
  35       if (narg.gt.3) then
  36         call getarg(4,ctemp)
  37         read(ctemp,*) minclust_size
  38       else
  39         minclust_size = 5
  40       endif
  41       if (narg.gt.4) then
  42         call getarg(5,outfile)
  43       else
  44         outfile = "klapaucjusz.out"
  45       endif
  46       open(2,file=outfile,status="unknown")
  47       print *,"nres",nres," nmodel",nmodel
  48       print *,"lencut",lencut," minclust_size",minclust_size
  49       print *," outfile ",outfile
  50       ifragpair_new=0
  51       nclust = 0
  52
  53       do i=1,nmodel
  54         do j=1,nmodel
  55           if (i.eq.j) cycle
  56           print *,"i",i," j",j
  57           read(*,'(a32,32x,1000i2)')model(i),(ifragpair(k,i,j),k=1,nres)
  58         enddo
  59       enddo
  60       do i=1,nmodel
  61         do j=i+1,nmodel
  62           write(*,'(2a32,1000i2)')model(i),model(j),
  63      &    (ifragpair(k,i,j),k=1,nres)
  64         enddo
  65       enddo
  66       isumcut=100
  67       iclust=0
  68       DO WHILE (isumcut.ge.lencut)
  69 c search the longest fragment (number 1)
  70       iflag=0
  71       isummax=0
  72       nfrag=0
  73       do i=1,nmodel
  74         do j=1,i-1
  75 c          write(*,'(2a32,1000i1)')model(i),model(j),(ifragpair(k,i,j)
  76 c     &       -ifragpair(k,j,i),k=1,nres)
  77           isump=0
  78           do k=1,nres
  79             kk = ifragpair(k,i,j)
  80             if (kk.gt.0) isump(kk)=isump(kk)+1
  81             if (kk.gt.nfrag) nfrag=kk
  82           enddo
  83           isum=0
  84           do k=1,nfrag
  85             if (isump(k).gt.isum) then
  86               isum=isump(k)
  87               kk=k
  88             endif
  89           enddo
  90 c          write(*,'(2a32,i10)')model(i),model(j),isum
  91           if (isum.gt.isummax) then
  92             isummax=isum
  93             imodel(1)=i
  94             imodel(2)=j
  95             kfmax=kk
  96           endif
  97         enddo
  98       enddo
  99       write (*,*) "nfrag",nfrag," kfmax",kfmax
 100 c      icut = max0(isummax*2/3,lencut)
 101       icut = lencut
 102       nflag=2
 103       iflag(imodel(1))=i
 104       iflag(imodel(2))=j
 105       write(*,*)"Maximum"
 106       write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
 107       i1 = imodel(1)
 108       i2 = imodel(2)
 109       isumcut=isummax
 110       if (isummax.ge.lencut) then
 111           kkmod=kfmax
 112           kmaxminmax = kfmax
 113           ifragpair_temp=0
 114           ifragpair_temp1=0
 115           do k=1,nres
 116             if (ifragpair(k,i1,i2).eq.kfmax) then
 117               ifragpair_temp(k)=ifragpair(k,i1,i2)
 118               ifragpair_temp1(k)=ifragpair(k,i1,i2)
 119             endif
 120           enddo
 121           overlap=.true.
 122       else
 123         overlap=.false.
 124       endif
 125       write (*,*) ">>>>>> dowhile"
 126       do while (overlap)
 127         do i=1,nmodel
 128           if (iflag(i).gt.0) cycle
 129           do j=1,nflag
 130             jj = imodel(j)
 131             isumf=0
 132             do k=1,nres
 133               kk = ifragpair(k,i,jj)
 134               if (ifragpair_temp1(k).ne.kfmax .or. kk.eq.0) cycle
 135               isumf(kk)=isumf(kk)+1
 136             enddo
 137             isumfmax(j)=0
 138             do k=1,nfrag
 139               if (isumf(k).gt.isumfmax(j)) then
 140                 isumfmax(j)=isumf(k)
 141                 kmax(j)=k
 142               endif
 143             enddo
 144           enddo ! j
 145           isumfmaxmin(i)=1000
 146           kmaxmin(i)=0
 147           jmaxmin(i)=0
 148           do j=1,nflag
 149             if (isumfmax(j).lt.isumfmaxmin(i)) then
 150               isumfmaxmin(i)=isumfmax(j)
 151               kmaxmin(i)=kmax(j)
 152               jmaxmin(i)=imodel(j)
 153             endif
 154           enddo
 155         enddo ! i
 156         isumfmaxminmax=0
 157         kmaxminmax=0
 158         imaxminmax=0
 159         jmaxminmax=0
 160         do i=1,nmodel
 161           if (iflag(i).gt.0) cycle
 162           if (isumfmaxmin(i).gt.isumfmaxminmax) then
 163             isumfmaxminmax=isumfmaxmin(i)
 164             imaxminmax=i
 165             kmaxminmax=kmaxmin(i)
 166             jmaxminmax=jmaxmin(i)
 167           endif
 168         enddo
 169         print *,"isumfmaxminmax",isumfmaxminmax," icut",icut
 170         if (isumfmaxminmax.lt.icut) then
 171           overlap=.false.
 172         else
 173           isumcut = isumfmaxminmax
 174           iimod=imaxminmax
 175           jjmod=jmaxminmax
 176           kkmod=kmaxminmax
 177           nflag=nflag+1
 178           iflag(imaxminmax)=imaxminmax
 179           imodel(nflag)=imaxminmax
 180           write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
 181           write(*,'(1000i1)')(ifragpair(k,imodel(1),imodel(2)),k=1,nres)
 182           write(*,*) "isumfmaxminmax",isumfmaxminmax
 183           print *,"kmaxminmax",kmaxminmax
 184           write(*,'(a32)')(model(imodel(i)),i=1,nflag)
 185           write(*,'(1000i1)')
 186      &     (ifragpair(k,imaxminmax,jmaxminmax),k=1,nres)
 187           ifragpair_temp=0
 188           do k=1,nres
 189             if (ifragpair(k,imaxminmax,jmaxminmax).eq.kmaxminmax
 190      &       .and. ifragpair_temp1(k).gt.0)
 191      &      ifragpair_temp(k)=ifragpair(k,imaxminmax,jmaxminmax)
 192           enddo
 193           do k=1,nres
 194             ifragpair_temp1(k)=ifragpair_temp(k)
 195           enddo
 196          print *,"ifragpair_temp,ifragpair_temp1"
 197          print '(1000i1)',(ifragpair_temp(k),k=1,nres)
 198          print '(1000i1)',(ifragpair_temp1(k),k=1,nres)
 199          do i=1,nmodel
 200            do j=i+1,nmodel
 201              write(*,'(2a32,1000i1)')model(i),model(j),
 202      &       (ifragpair(k,i,j),k=1,nres)
 203            enddo
 204          enddo
 205         endif
 206       enddo ! while
 207       print *,"iimod",iimod," jjmod",jjmod," kkmod",kkmod,
 208      &   " isumcut",isumcut
 209       if (isumcut.ge.icut) then
 210           print *,"ifragpair_temp1 before zeroing"
 211           print '(64x,1000i1)',(ifragpair_temp1(k),k=1,nres)
 212           do i=1,nflag
 213             do j=1,i-1
 214               ii = imodel(i)
 215               jj = imodel(j)
 216               do k=1,nres
 217 c                if (ifragpair_temp(k).eq.kmaxminmax .and.
 218 c     &           ifragpair_temp1(k).eq.kfmax) then
 219                 if (ifragpair_temp1(k).gt.0) then
 220                   ifragpair(k,ii,jj)=0
 221                   ifragpair(k,jj,ii)=0
 222                 endif
 223               enddo
 224             enddo
 225           enddo
 226          print *,">>>>>>>>>>>>>>> IFRAGPAIR AFTER ZEROING <<<<<<<<<<<"
 227          do i=1,nmodel
 228            do j=i+1,nmodel
 229              write(*,'(2a32,1000i1)')model(i),model(j),
 230      &       (ifragpair(k,i,j),k=1,nres)
 231            enddo
 232          enddo
 233
 234         if (nflag.ge.minclust_size) then
 235
 236         nclust=nclust+1
 237         print *,"cluster",nclust
 238         ninclust(nclust)=nflag
 239         nlenclust(nclust)=isumcut
 240         do i=1,nflag
 241           icluster(i,nclust)=imodel(i)
 242         enddo
 243         ii=0
 244         do k=1,nres
 245 c          if (ifragpair_temp(k).eq.kkmod .and.
 246 c     &     ifragpair_temp1(k).eq.kfmax) then
 247           if (ifragpair_temp1(k).gt.0) then
 248             ii=ii+1
 249 c            print *,k
 250             ifrag(ii,nclust)=k
 251           endif
 252         enddo
 253         if (ii.ne.nlenclust(nclust)) write(*,*) "CHUJ NASTAPIL!!!",
 254      &    ii,nlenclust(nclust)
 255         do i=1,nflag
 256           do j=1,i-1
 257             do k=1,nres
 258               if (ifragpair_temp1(k).gt.0) then
 259                 ifragpair_new(k,imodel(i),imodel(j))=nclust
 260                 ifragpair_new(k,imodel(j),imodel(i))=nclust
 261               endif
 262             enddo
 263           enddo
 264         enddo
 265
 266         endif ! nflag > 2
 267
 268       endif
 269       ENDDO ! while
 270       do i=1,nmodel
 271         do j=i+1,nmodel
 272           write(*,'(2a32,1000a1)')model(i),model(j),
 273      &      (charm(ifragpair_new(k,i,j)),k=1,nres)
 274         enddo
 275       enddo
 276       write (*,*) "nclust",nclust
 277       do i=1,nclust
 278         write(*,*) "Cluster",i," ninclust",ninclust(i),
 279      &   " length",nlenclust(i)
 280         write(*,'(a32)') (model(icluster(j,i)),j=1,ninclust(i))
 281         write(*,'(i5)') (ifrag(j,i),j=1,nlenclust(i))
 282       enddo
 283       imodelclust=0
 284       do i=1,nclust
 285         do j=1,ninclust(i)
 286           imodelclust(icluster(j,i))=1
 287         enddo
 288       enddo
 289       print *,"Models in clusters"
 290       nmodel_out = 0
 291       do i=1,nmodel
 292         if (imodelclust(i).eq.1) then
 293           print '(a32)',model(i)
 294           nmodel_out = nmodel_out + 1
 295         endif
 296       enddo
 297       print *,"Models excluded"
 298       do i=1,nmodel
 299         if (imodelclust(i).eq.0) print '(a32)',model(i)
 300       enddo
 301 c Write model information
 302       write(2,'(i5)') nmodel_out,nclust
 303       ii=0
 304       iindex=0
 305       do i=1,nmodel
 306         if (imodelclust(i).eq.1) then
 307           ii=ii+1
 308           write(2,'(a32)') model(i)
 309           iindex(i)=ii
 310         endif
 311       enddo
 312 c Write cluster information
 313       do i=1,nclust
 314         write(2,'(2i5)') ninclust(i),nlenclust(i)
 315         write(2,'(16i5)') (iindex(icluster(k,i)),k=1,ninclust(i))
 316         write(2,'(16i5)') (ifrag(k,i),k=1,nlenclust(i))
 317       enddo
 318       end
 319 c---------------------------------------------------------------------------------
 320       character function charm(i)
 321       integer i
 322       if (i.eq.0) then
 323         charm = "."
 324       else
 325         charm = char(ichar("A")+i-1)
 326       endif
 327       return
 328       end