From: Cezary Czaplewski Date: Thu, 5 Apr 2018 13:09:39 +0000 (+0200) Subject: program genereating HOMOL_FRAG consensus fragments X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=f993ab6738b32897a396baeb156285885fddc002 program genereating HOMOL_FRAG consensus fragments --- diff --git a/source/fragments/fitsq.f b/source/fragments/fitsq.f new file mode 100644 index 0000000..b42473a --- /dev/null +++ b/source/fragments/fitsq.f @@ -0,0 +1,313 @@ + subroutine fitsq(rms,x,y,nn,t,b) +c x and y are the vectors of coordinates (dimensioned (3,n)) of the two +c structures to be superimposed. nn is 3*n, where n is the number of +c points. t and b are respectively the translation vector and the +c rotation matrix that transforms the second set of coordinates to the +c frame of the first set. +c eta = machine-specific variable + +c dimension x(1),y(1),t(3) + dimension x(3*nn),y(3*nn),t(3) + dimension b(3,3),q(3,3),r(3,3),v(3),xav(3),yav(3),e(3),c(3,3) + eta = z00100000 +c small=25.0*rmdcon(3) +c small=25.0*eta +c small=25.0*10.e-10 +c the following is a very lenient value for 'small' + small = 0.0001 + fn=nn + do 10 i=1,3 + xav(i)=0.0 + yav(i)=0.0 + do 10 j=1,3 + 10 b(j,i)=0.0 + nc=0 +c + do 30 n=1,nn + do 20 i=1,3 +c print*,'x = ',x(nc+i),' y = ',y(nc+i) + xav(i)=xav(i)+x(nc+i)/fn + 20 yav(i)=yav(i)+y(nc+i)/fn + 30 nc=nc+3 +c + +* print*,'xav = ',(xav(j),j=1,3) +* print*,'yav = ',(yav(j),j=1,3) + + nc=0 + rms=0.0 + do 50 n=1,nn + do 40 i=1,3 + rms=rms+((x(nc+i)-xav(i))**2+(y(nc+i)-yav(i))**2)/fn + do 40 j=1,3 + b(j,i)=b(j,i)+(x(nc+i)-xav(i))*(y(nc+j)-yav(j))/fn + 40 c(j,i)=b(j,i) + 50 nc=nc+3 + call sivade(b,q,r,d) + sn3=sign(1.0,d) + do 120 i=1,3 + do 120 j=1,3 + 120 b(j,i)=-q(j,1)*r(i,1)-q(j,2)*r(i,2)-sn3*q(j,3)*r(i,3) + call mvvad(b,xav,yav,t) + do 130 i=1,3 + do 130 j=1,3 + rms=rms+2.0*c(j,i)*b(j,i) + 130 b(j,i)=-b(j,i) + if (abs(rms).gt.small) go to 140 +* write (6,301) + return + 140 if (rms.gt.0.0) go to 150 +* write (6,303) rms + rms=0.0 +* stop +* 150 write (6,302) sqrt(rms) + 150 rms=sqrt(rms) + return + 301 format (5x,'rms deviation negligible') + 302 format (5x,'rms deviation ',f14.6) + 303 format (//,5x,'negative ms deviation - ',f14.6) + end + subroutine sivade(x,q,r,dt) +c computes q,e and r such that q(t)xr = diag(e) + dimension x(3,3),q(3,3),r(3,3),e(3) + dimension h(3,3),p(3,3),u(3,3),d(3) + eta = z00100000 + small=25.0*10.e-10 +c small=25.0*eta +c small=2.0*rmdcon(3) + xnrm=0.0 + do 20 i=1,3 + do 10 j=1,3 + xnrm=xnrm+x(j,i)*x(j,i) + u(j,i)=0.0 + r(j,i)=0.0 + 10 h(j,i)=0.0 + u(i,i)=1.0 + 20 r(i,i)=1.0 + xnrm=sqrt(xnrm) + do 110 n=1,2 + xmax=0.0 + do 30 j=n,3 + 30 if (abs(x(j,n)).gt.xmax) xmax=abs(x(j,n)) + a=0.0 + do 40 j=n,3 + h(j,n)=x(j,n)/xmax + 40 a=a+h(j,n)*h(j,n) + a=sqrt(a) + den=a*(a+abs(h(n,n))) + d(n)=1.0/den + h(n,n)=h(n,n)+sign(a,h(n,n)) + do 70 i=n,3 + s=0.0 + do 50 j=n,3 + 50 s=s+h(j,n)*x(j,i) + s=d(n)*s + do 60 j=n,3 + 60 x(j,i)=x(j,i)-s*h(j,n) + 70 continue + if (n.gt.1) go to 110 + xmax=amax1(abs(x(1,2)),abs(x(1,3))) + h(2,3)=x(1,2)/xmax + h(3,3)=x(1,3)/xmax + a=sqrt(h(2,3)*h(2,3)+h(3,3)*h(3,3)) + den=a*(a+abs(h(2,3))) + d(3)=1.0/den + h(2,3)=h(2,3)+sign(a,h(2,3)) + do 100 i=1,3 + s=0.0 + do 80 j=2,3 + 80 s=s+h(j,3)*x(i,j) + s=d(3)*s + do 90 j=2,3 + 90 x(i,j)=x(i,j)-s*h(j,3) + 100 continue + 110 continue + do 130 i=1,3 + do 120 j=1,3 + 120 p(j,i)=-d(1)*h(j,1)*h(i,1) + 130 p(i,i)=1.0+p(i,i) + do 140 i=2,3 + do 140 j=2,3 + u(j,i)=u(j,i)-d(2)*h(j,2)*h(i,2) + 140 r(j,i)=r(j,i)-d(3)*h(j,3)*h(i,3) + call mmmul(p,u,q) + 150 np=1 + nq=1 + if (abs(x(2,3)).gt.small*(abs(x(2,2))+abs(x(3,3)))) go to 160 + x(2,3)=0.0 + nq=nq+1 + 160 if (abs(x(1,2)).gt.small*(abs(x(1,1))+abs(x(2,2)))) go to 180 + x(1,2)=0.0 + if (x(2,3).ne.0.0) go to 170 + nq=nq+1 + go to 180 + 170 np=np+1 + 180 if (nq.eq.3) go to 310 + npq=4-np-nq + if (np.gt.npq) go to 230 + n0=0 + do 220 n=np,npq + nn=n+np-1 + if (abs(x(nn,nn)).gt.small*xnrm) go to 220 + x(nn,nn)=0.0 + if (x(nn,nn+1).eq.0.0) go to 220 + n0=n0+1 + go to (190,210,220),nn + 190 do 200 j=2,3 + 200 call givns(x,q,1,j) + go to 220 + 210 call givns(x,q,2,3) + 220 continue + if (n0.ne.0) go to 150 + 230 nn=3-nq + a=x(nn,nn)*x(nn,nn) + if (nn.gt.1) a=a+x(nn-1,nn)*x(nn-1,nn) + b=x(nn+1,nn+1)*x(nn+1,nn+1)+x(nn,nn+1)*x(nn,nn+1) + c=x(nn,nn)*x(nn,nn+1) + dd=0.5*(a-b) + xn2=c*c + rt=b-xn2/(dd+sign(sqrt(dd*dd+xn2),dd)) + y=x(np,np)*x(np,np)-rt + z=x(np,np)*x(np,np+1) + do 300 n=np,nn + if (abs(y).lt.abs(z)) go to 240 + t=z/y + c=1.0/sqrt(1.0+t*t) + s=c*t + go to 250 + 240 t=y/z + s=1.0/sqrt(1.0+t*t) + c=s*t + 250 do 260 j=1,3 + v=x(j,n) + w=x(j,n+1) + x(j,n)=c*v+s*w + x(j,n+1)=-s*v+c*w + a=r(j,n) + b=r(j,n+1) + r(j,n)=c*a+s*b + 260 r(j,n+1)=-s*a+c*b + y=x(n,n) + z=x(n+1,n) + if (abs(y).lt.abs(z)) go to 270 + t=z/y + c=1.0/sqrt(1.0+t*t) + s=c*t + go to 280 + 270 t=y/z + s=1.0/sqrt(1.0+t*t) + c=s*t + 280 do 290 j=1,3 + v=x(n,j) + w=x(n+1,j) + a=q(j,n) + b=q(j,n+1) + x(n,j)=c*v+s*w + x(n+1,j)=-s*v+c*w + q(j,n)=c*a+s*b + 290 q(j,n+1)=-s*a+c*b + if (n.ge.nn) go to 300 + y=x(n,n+1) + z=x(n,n+2) + 300 continue + go to 150 + 310 do 320 i=1,3 + 320 e(i)=x(i,i) + 330 n0=0 + do 360 i=1,3 + if (e(i).ge.0.0) go to 350 + e(i)=-e(i) + do 340 j=1,3 + 340 q(j,i)=-q(j,i) + 350 if (i.eq.1) go to 360 + if (abs(e(i)).lt.abs(e(i-1))) go to 360 + call switch(i,1,q,r,e) + n0=n0+1 + 360 continue + if (n0.ne.0) go to 330 + if (abs(e(3)).gt.small*xnrm) go to 370 + e(3)=0.0 + if (abs(e(2)).gt.small*xnrm) go to 370 + e(2)=0.0 + 370 dt=det(q(1,1),q(1,2),q(1,3))*det(r(1,1),r(1,2),r(1,3)) +* write (1,501) (e(i),i=1,3) + return + 501 format (/,5x,'singular values - ',3e15.5) + end + subroutine givns(a,b,m,n) + dimension a(3,3),b(3,3) + if (abs(a(m,n)).lt.abs(a(n,n))) go to 10 + t=a(n,n)/a(m,n) + s=1.0/sqrt(1.0+t*t) + c=s*t + go to 20 + 10 t=a(m,n)/a(n,n) + c=1.0/sqrt(1.0+t*t) + s=c*t + 20 do 30 j=1,3 + v=a(m,j) + w=a(n,j) + x=b(j,m) + y=b(j,n) + a(m,j)=c*v-s*w + a(n,j)=s*v+c*w + b(j,m)=c*x-s*y + 30 b(j,n)=s*x+c*y + return + end + subroutine switch(n,m,u,v,d) + dimension u(3,3),v(3,3),d(3) + do 10 i=1,3 + tem=u(i,n) + u(i,n)=u(i,n-1) + u(i,n-1)=tem + if (m.eq.0) go to 10 + tem=v(i,n) + v(i,n)=v(i,n-1) + v(i,n-1)=tem + 10 continue + tem=d(n) + d(n)=d(n-1) + d(n-1)=tem + return + end +c subroutine mvvad(a,b,c,d) + subroutine mvvad(b,xav,yav,t) + dimension b(3,3),xav(3),yav(3),t(3) +c dimension a(3,3),b(3),c(3),d(3) +c do 10 j=1,3 +c d(j)=c(j) +c do 10 i=1,3 +c 10 d(j)=d(j)+a(j,i)*b(i) + do 10 j=1,3 + t(j)=yav(j) + do 10 i=1,3 + 10 t(j)=t(j)+b(j,i)*xav(i) + return + end + real function det (a,b,c) + dimension a(3),b(3),c(3) + det=a(1)*(b(2)*c(3)-b(3)*c(2))+a(2)*(b(3)*c(1)-b(1)*c(3)) + 1 +a(3)*(b(1)*c(2)-b(2)*c(1)) + return + end + subroutine mmmul(a,b,c) + dimension a(3,3),b(3,3),c(3,3) + do 10 i=1,3 + do 10 j=1,3 + c(i,j)=0.0 + do 10 k=1,3 + 10 c(i,j)=c(i,j)+a(i,k)*b(k,j) + return + end + subroutine matvec(uvec,tmat,pvec,nback) + real*4 tmat(3,3),uvec(3,nback), pvec(3,nback) +c + do 2 j=1,nback + do 1 i=1,3 + uvec(i,j) = 0.0 + do 1 k=1,3 + 1 uvec(i,j)=uvec(i,j)+tmat(i,k)*pvec(k,j) + 2 continue + return + end diff --git a/source/fragments/klapaucjusz-longest.f b/source/fragments/klapaucjusz-longest.f new file mode 100644 index 0000000..41f49a3 --- /dev/null +++ b/source/fragments/klapaucjusz-longest.f @@ -0,0 +1,208 @@ + parameter (maxmodel=20) + parameter (maxres=1000) + character*32 model(maxmodel) + integer imodel(maxmodel) + integer ifragpair(maxres,maxmodel,maxmodel), + & ifragpair_new(maxres,maxmodel,maxmodel), + & 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) + logical overlap + + nmodel=5 + nres=114 + ifragpair_new=0 + lencut=7 + + do i=1,nmodel + do j=1,nmodel + if (i.eq.j) cycle + read(*,'(a32,33x,1000i1)')model(i),(ifragpair(k,i,j),k=1,nres) + enddo + enddo + 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 + isumcut=100 + iclust=0 + DO WHILE (isumcut.ge.lencut) +c search the longest fragment (number 1) + iflag=0 + isummax=0 + nfrag=0 + do i=1,nmodel + do j=1,i-1 +c write(*,'(2a32,1000i1)')model(i),model(j),(ifragpair(k,i,j) +c & -ifragpair(k,j,i),k=1,nres) + isump=0 + do k=1,nres + kk = ifragpair(k,i,j) + if (kk.gt.0) isump(kk)=isump(kk)+1 + if (kk.gt.nfrag) nfrag=kk + enddo + isum=0 + do k=1,nfrag + if (isump(k).gt.isum) then + isum=isump(k) + kk=k + endif + enddo +c write(*,'(2a32,i10)')model(i),model(j),isum + if (isum.gt.isummax) then + isummax=isum + imodel(1)=i + imodel(2)=j + kfmax=kk + endif + enddo + enddo + write (*,*) "nfrag",nfrag," kfmax",kfmax +c icut = max0(isummax*2/3,lencut) + icut = 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. + i1 = imodel(1) + i2 = imodel(2) + isumcut=0 + do while (overlap) + do i=1,nmodel + if (iflag(i).gt.0) cycle + do j=1,nflag + jj = imodel(j) + isumf=0 + do k=1,nres + kk = ifragpair(k,i,jj) + if (ifragpair(k,i1,i2).ne.kfmax .or. kk.eq.0) cycle + isumf(kk)=isumf(kk)+1 + enddo + isumfmax(j)=0 + do k=1,nfrag + if (isumf(k).gt.isumfmax(j)) then + isumfmax(j)=isumf(k) + kmax(j)=k + endif + enddo + enddo ! j + isumfmaxmin(i)=1000 + kmaxmin(i)=0 + jmaxmin(i)=0 + do j=1,nflag + if (isumfmax(j).lt.isumfmaxmin(i)) then + isumfmaxmin(i)=isumfmax(j) + kmaxmin(i)=kmax(j) + jmaxmin(i)=imodel(j) + endif + enddo + enddo ! i + isumfmaxminmax=0 + kmaxminmax=0 + imaxminmax=0 + jmaxminmax=0 + do i=1,nmodel + if (iflag(i).gt.0) cycle + if (isumfmaxmin(i).gt.isumfmaxminmax) then + isumfmaxminmax=isumfmaxmin(i) + imaxminmax=i + kmaxminmax=kmaxmin(i) + jmaxminmax=jmaxmin(i) + endif + enddo + print *,"isumfmaxminmax",isumfmaxminmax," icut",icut + if (isumfmaxminmax.lt.icut) then + overlap=.false. + else + isumcut = isumfmaxminmax + iimod=imaxminmax + jjmod=jmaxminmax + kkmod=kmaxminmax + 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(*,*) "isumfmaxminmax",isumfmaxminmax + print *,"kmaxminmax",kmaxminmax + write(*,'(a32)')(model(imodel(i)),i=1,nflag) + write(*,'(1000i1)') + & (ifragpair(k,imaxminmax,jmaxminmax),k=1,nres) + do k=1,nres + ifragpair_temp(k)=ifragpair(k,imaxminmax,jmaxminmax) + ifragpair_temp1(k)=ifragpair(k,i1,i2) + enddo + 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 + ifragpair(k,ii,jj)=0 + ifragpair(k,jj,ii)=0 + endif + enddo + enddo + enddo + 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 + nclust=nclust+1 + print *,"cluster",nclust + ninclust(nclust)=nflag + nlenclust(nclust)=isumcut + do i=1,nflag + icluster(i,nclust)=imodel(i) + enddo + ii=0 + do k=1,nres + if (ifragpair_temp(k).eq.kkmod .and. + & ifragpair_temp1(k).eq.kfmax) then + ii=ii+1 + print *,k + ifrag(ii,nclust)=k + endif + enddo + if (ii.ne.nlenclust(nclust)) write(*,*) "CHUJ NASTAPIL!!!", + & ii,nlenclust(nclust) + 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 + ifragpair_new(k,imodel(i),imodel(j))=nclust + ifragpair_new(k,imodel(j),imodel(i))=nclust + endif + enddo + enddo + enddo + 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) + enddo + enddo + write (*,*) "nclust",nclust + do i=1,nclust + write(*,*) "Cluster",i," ninclust",ninclust(i) + write(*,'(a32)') (model(icluster(j,i)),j=1,ninclust(i)) + write(*,'(i5)') (ifrag(j,i),j=1,nlenclust(i)) + enddo + end diff --git a/source/fragments/run.csh b/source/fragments/run.csh new file mode 100755 index 0000000..2e6d8aa --- /dev/null +++ b/source/fragments/run.csh @@ -0,0 +1,10 @@ +#!/bin/csh +foreach f1 (*[1-5]) + foreach f2 (*[1-5]) + if ($f1 != $f2) then + echo -n $f1 $f2 |awk '{printf "%-32s%-32s",$1,$2}' >> all.txt + ./f $f1 $f2 ooo_${f1}_${f2}_ooo 1 114 >& /dev/null + tail -1 ooo_${f1}_${f2}_ooo >> all.txt + endif +end +end diff --git a/source/fragments/suppdb_ca_dis.f b/source/fragments/suppdb_ca_dis.f new file mode 100644 index 0000000..50b09ce --- /dev/null +++ b/source/fragments/suppdb_ca_dis.f @@ -0,0 +1,397 @@ + dimension s1(3,10000),s2(3,10000),ca1(3,10000),ca2(3,10000) + dimension xyz1(3,5000),xyz2(3,5000),waby(5000),ires(5000), + & atnam(5000),rms(500,500),imin(5000),imax(5000),isel(5000), + & isel_t(5000),issel(5000),diff(5000),iperm(5000) + character*4 atnam + dimension t(3),b(3,3),z(3,5000) + integer resnum,resnum1,resnum2,resnum3,resnum4,frag(5000) + character*80 karta + character*80 fnam1,fnam2,fnam3 + logical insel + ifrag=1 + nnn = iargc() + if (nnn.lt.3) then + stop + & 'Usage: suppdb pdb_file_1 pdb_file_2 out_file [resnum1 resnum2]' + endif + call getarg(1,fnam1) + call getarg(2,fnam2) + call getarg(3,fnam3) + open (2,file=fnam3,status='unknown') + if (nnn.lt.5) then + print *,'resnum1,resnum2' + read (*,*) resnum1,resnum2 + else + call getarg(4,karta) + read (karta,*) resnum1 + call getarg(5,karta) + read (karta,*) resnum2 + endif + if (nnn.gt.5) then + call getarg(6,karta) + read (karta,*) resnum3 + else + resnum3=resnum1 + endif + resnum4=resnum3+resnum2-resnum1 + print *,'resnum1',resnum1,' resnum2',resnum2 + print *,'resnum3',resnum3,' resnum4',resnum4 + print *,'opening file ',fnam1 + open (1,file=fnam1,status='old',err=111) + iat=0 + isup=0 + write (2,'(/2a)') 'Coordinates of the first set from ', + & fnam1(:len_trim(fnam1)) + do irec=1,10000 + read (1,'(a80)',end=10) karta + if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and. + & karta(14:15).eq.'CA') then + read (karta(23:26),*) resnum + if (resnum.ge.resnum1 .and. resnum.le.resnum2) then + iat=iat+1 + read (karta(14:16),'(a3)') atnam(iat) + read (karta(31:54),'(3f8.3)') xyz1(1,iat),xyz1(2,iat), + 1 xyz1(3,iat) + read (karta(23:26),'(i4)') ires(iat) + read (karta(18:20),'(a3)') waby(iat) + if (atnam(iat).eq.'CA ') then + write (2,'(a)') karta + isup=isup+1 + do i=1,3 + ca1(i,isup)=xyz1(i,iat) + s1(i,isup)=xyz1(i,iat) + enddo + endif + endif + endif + enddo + 10 nsup=isup + nat=iat + print *,'from first file: nat=',nat + close(1) + print *,'opening file ',fnam2 + open (1,file=fnam2,status='old',err=112) + iat=0 + isup=0 + write (2,'(/2a)') 'Coordinates of the first set from ', + & fnam2(:len_trim(fnam2)) + do irec=1,10000 + read (1,'(a80)',end=20) karta + if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and. + & karta(14:15).eq.'CA') then + read (karta(23:26),*) resnum + if (resnum.ge.resnum3 .and. resnum.le.resnum4) then + iat=iat+1 + read (karta(31:54),'(3f8.3)') xyz2(1,iat),xyz2(2,iat), + 1 xyz2(3,iat) + if (atnam(iat).eq.'CA ') then + write (2,'(a)') karta + isup=isup+1 + do i=1,3 + ca2(i,isup)=xyz2(i,iat) + s2(i,isup)=xyz2(i,iat) + enddo + endif + endif + endif + enddo + 20 continue + print *,'from second file: nat=',iat +! write (2,'(//4a8)') '# length',' start',' rms' + + + do i=1,nsup + issel(i)=i + enddo + + nsup_t = nsup + + nsup_old=0 + + DO WHILE (NSUP>4 .and. nsup.ne.nsup_old) + + nsup_old = nsup + + do is1=4,nsup +! write (2,'(a)') + do is2=1,nsup-is1+1 + call fitsq(rms(is1,is2),s2(1,is2),s1(1,is2),is1,t,b) +! write (2,'(2(i4,4x),2f8.2)') is1,is2+resnum1-1,rms(is1,is2) + enddo + enddo + write (2,'(//2a8,3a8)') '# length',' start','rmsmin', + & 'rmsmax','rmsave' + do i=4,nsup + imin(i)=1 + imax(i)=1 + rmin=rms(i,1) + rmax=rms(i,1) + rmsave=rms(i,1)**2 + do j=2,nsup-i+1 + if ( rms(i,j).lt.rmin) then + imin(i)=j + rmin=rms(i,j) + else if ( rms(i,j).gt.rmax) then + imax(i)=j + rmax=rms(i,j) + endif + rmsave=rmsave+rms(i,j)*rms(i,j) + enddo + rmsave=sqrt(rmsave/(nsup-i+1)) + write (2,'(2(i4,4x),3(f8.2))') i,imin(i)+resnum1-1, + & rms(i,imin(i)), ! rms(i,imin(i))/sqrt(float(i)), + & rms(i,imax(i)), ! rms(i,imax(i))/sqrt(float(i)), + & rmsave ! ,rmsave/sqrt(float(i)) + enddo + rmscut=4.0d0 + rmscut1=1.0d0*rmscut + do i4=nsup,4,-1 + if (rms(i4,imin(i4)).le.rmscut) exit + enddo + call fitsq(rrms,s2(1,imin(i4)),s1(1,imin(i4)),i4,t,b) + nsup_t=nsup + do i=1,nsup + isel(i)=issel(i) + enddo + rrms=1000.0d0 +c do while (rrms.gt.rmscut) + do while (rrms.gt.rmscut) + nnsup_t=0 + std=0.0d0 + diffmax=-10000.0d0 + do i=1,nsup_t + ii=isel(i) + call mvvad(b,ca2(1,ii),t,z(1,i)) + diffi=0.0d0 + do j=1,3 + diffi=diffi+(z(j,i)-ca1(j,ii))**2 + enddo + std=std+diffi + diffi=sqrt(diffi) +! write (2,*) i,ii,diffi + diff(i)=diffi + if (diffi.gt.diffmax) diffmax=diffi + enddo + do i=1,nsup_t + iperm(i)=i + enddo + write(2,*) 'diff' + write(2,*) (diff(i),i=1,nsup_t) + call mysort(nsup_t,diff,iperm) + write(2,*) 'diff' + write(2,*) (diff(i),i=1,nsup_t) + std=sqrt(std/nsup_t) + write (2,*) 'std=',std + write (*,*) 'std=',std + sstd=dmin1(1.5d0*std,0.9999d0*diffmax) +! do i=1,nsup_t +! ii=isel(i) +! if (diff(i).lt.sstd .or. diff(i).lt.rmscut) then +! nnsup_t=nnsup_t+1 +! isel(nnsup_t)=ii +! endif +! enddo + sumdiff=0.0d0 + do i=1,nsup_t + ii=isel(iperm(i)) +c sumdiff=sumdiff+diff(i)**2 +c if (dsqrt(sumdiff/dfloat(i)).lt.rmscut1) then + if (diff(i).lt.rmscut1) then + nnsup_t=nnsup_t+1 + isel_t(nnsup_t)=ii + endif + enddo + nsup_t=nnsup_t + do i=1,nsup_t + iperm(i)=i + isel(i)=isel_t(i) + enddo + if (nsup_t.eq.0 .or. nnsup_t.eq.0) exit + write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t + write (2,'(20i4)') (isel(i),i=1,nnsup_t) + call mysort(nsup_t,isel,iperm) + print *,'nsup_t',nsup_t,' nnsup_t',nnsup_t + print '(20i4)',(isel(i),i=1,nnsup_t) + write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t + write (2,'(20i4)') (isel(i),i=1,nnsup_t) + len_min=4 + i=1 + len=1 + do i=1,nsup_t + ii=isel(i) + do j=1,3 + s1(j,i)=ca1(j,ii) + s2(j,i)=ca2(j,ii) + enddo + enddo + do i=1,nsup_t + ii=isel(i) + do j=1,3 + s1(j,i)=ca1(j,ii) + s2(j,i)=ca2(j,ii) + enddo + enddo + if (nsup_t.le.2) exit + print *,"before fitsq: nsup_r",nsup_t + call fitsq(rrms,s2(1,1),s1(1,1),nsup_t,t,b) + write (2,*) 'nsup_t=',nsup_t,' rrms=',rrms + write (*,*) 'nsup_t=',nsup_t,' rrms=',rrms + enddo + i=1 + do while (i.le.nnsup_t) + write (*,*) 'i=',i,' isel',isel(i),' isel1',isel(i+1) + if (i.eq.nnsup_t .or. isel(i+1)-isel(i).gt.1) then + if (len.lt.len_min) then + write (*,'(20i4)') (isel(k),k=1,nnsup_t) + write (*,*) 'len=',len,' len_min=',len_min + do j=i+1,nnsup_t + isel(j-len)=isel(j) + enddo + nnsup_t=nnsup_t-len +! write (2,'(20i4)') (isel(k),k=1,nnsup_t) + i=i-len+1 + len=1 + else + len=1 + i=i+1 + endif + else + len=len+1 + i=i+1 + endif +! write (2,*) 'i=',i,' len=',len + enddo + print *,'final nsup_t',nsup_t + print *,"final isel" + print '(20i4)',(isel(i),i=1,nnsup_t) + write (2,*) 'final nsup_t',nsup_t + write (2,*) "final isel" + write (2,'(20i4)') (isel(i),i=1,nnsup_t) + if (nnsup_t.gt.10) then + do i=1,nnsup_t + frag(isel(i))=ifrag + enddo + ifrag=ifrag+1 + endif + print *,"----- end of iteration" + write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t + nsup_t=nnsup_t + write (2,'(/4a)') 'REMARK superimposed coordinates from files ', + &fnam1(:len_trim(fnam1)),' and ',fnam2(:len_trim(fnam2)) + write (2,'(a,f6.2,a)') + & 'REMARK Calpha rms deviation:',rrms, + & ' angstroms' + write (2,'(a,I5)') 'REMARK number of superposed residues',nsup_t + i1=isel(1) + do i=2,nnsup_t + if (isel(i)-isel(i-1).gt.1) then + write (2,'(a,i4,a,i4,a,f5.1,a)') + & 'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1 + i1=isel(i) + endif + enddo + write (2,'(a,i4,a,i4,a,f5.1,a)') + & 'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1 + do i=1,nsup_t + ii=isel(i) + write (2,100) ii+resnum1-1,atnam(ii),waby(ii),ires(ii), + & (xyz1(j,ii),j=1,3) + enddo + write (2,'(a)') 'TER' + do i=1,nnsup_t + ii=isel(i) + call mvvad(b,xyz2(1,ii),t,z(1,ii)) + write (2,100) ii+resnum3-1,atnam(ii),waby(ii),ires(ii), + & (z(j,ii),j=1,3) + enddo + write (2,'(a)') 'TER' + + ind=0 + do i=1,nsup + insel=.false. + do j=1,nsup_t + if (issel(i).eq.isel(j)) then + insel=.true. + exit + endif + enddo + if (.not.insel) then + ind=ind+1 + issel(ind)=issel(i) + endif + enddo + nsup=nsup-nsup_t + print *,"nsup",nsup + if (nsup.gt.0) print *,"issel",(issel(i),i=1,nsup) + do i=1,nsup + do j=1,3 + s1(j,i)=ca1(j,issel(i)) + s2(j,i)=ca2(j,issel(i)) + enddo + enddo + + ENDDO + + write (2,'(5000i1)') (frag(i),i=1,resnum2-resnum1+1) + stop 777 + 111 print '(2a)','Error opening input file ',fnam1(:len_trim(fnam1)) + stop 111 + 112 print '(2a)','Error opening input file ',fnam2(:len_trim(fnam2)) + stop 112 + 100 format ('ATOM',i7,2x,a2,2x,a3,i6,4x,3f8.3) + end +c---------------------------------------------------------------------- + subroutine distout(ncon,rms) + parameter(ncol=10) + dimension rms(200,200) + dimension b(ncol) + iout=2 +c icant(i,j)=((ncon+ncon-j)*(j-1))/2+i-j + write (iout,'(a)') 'The distance matrix' + do 1 i=1,ncon,ncol + nlim=min0(i+ncol-1,ncon) + write (iout,1000) (k,k=i,nlim) + write (iout,'(8h--------,10a)') ('-------',k=i,nlim) + 1000 format (/8x,10(i4,3x)) + 1020 format (/1x,80(1h-)/) + do 2 j=i,ncon + jlim=min0(j,nlim) + if (jlim.eq.j) then + b(jlim-i+1)=0.0d0 + jlim1=jlim-1 + else + jlim1=jlim + endif + do 3 k=i,jlim1 + 3 b(k-i+1)=rms(k,j) + write (iout,1010) j,(b(k),k=1,jlim-i+1) + 2 continue + 1 continue + 1010 format (i5,3x,10(f6.2,1x)) + return + end +!-------------------------------------------------- + subroutine mysort(n, x, ipermut) + implicit none + real x(n),xsort(n) + integer ipermut(n) + integer i,j,imax,ipm,n + real xtemp + do i=1,n + xtemp=x(i) + imax=i + do j=i+1,n + if (x(j).lt.xtemp) then + imax=j + xtemp=x(j) + endif + enddo + x(imax)=x(i) + x(i)=xtemp + ipm=ipermut(imax) + ipermut(imax)=ipermut(i) + ipermut(i)=ipm + enddo + return + end +