From: Cezary Czaplewski Date: Thu, 5 Apr 2018 13:11:22 +0000 (+0200) Subject: Adam's corrections to HOMOL_FRAG program X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=182795727c2f5bc66b50db6c87693dc4770ec06c 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 --- 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