--- /dev/null
+ 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
--- /dev/null
+ 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
--- /dev/null
+#!/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
--- /dev/null
+ 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
+
projects / unres.git / commitdiff