source/fragments/klapaucjusz-longest.f

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