X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2FTMscore_subroutine.f;fp=source%2Funres%2Fsrc_CSA_DiL%2FTMscore_subroutine.f;h=0000000000000000000000000000000000000000;hp=8e6ee9a24c75e221a025edfd2c13e0132c3b403b;hb=2a226bfc86eabc6e4eae0c3ad1cbc3cb5417a05a;hpb=a0e685f844163003749ba91dfbf4644bcc8cfa30

diff --git a/source/unres/src_CSA_DiL/TMscore_subroutine.f b/source/unres/src_CSA_DiL/TMscore_subroutine.f
deleted file mode 100644
index 8e6ee9a..0000000
--- a/source/unres/src_CSA_DiL/TMscore_subroutine.f
+++ /dev/null
@@ -1,536 +0,0 @@
-*************************************************************************
-*************************************************************************
-*     This is a subroutine to compare two structures and find the 
-*     superposition that has the maximum TM-score.
-*     Reference: Yang Zhang, Jeffrey Skolnick, Proteins 2004 57:702-10.
-*
-*     Explanations:
-*     L1--Length of the first structure
-*     (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
-*     n1(i)--Residue sequence number of i'th residue at the first structure
-*     L2--Length of the second structure
-*     (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
-*     n2(i)--Residue sequence number of i'th residue at the second structure
-*     TM--TM-score of the comparison
-*     Rcomm--RMSD of two structures in the common aligned residues
-*     Lcomm--Length of the common aligned regions
-*
-*     Note: 
-*     1, Always put native as the second structure, by which TM-score
-*        is normalized.
-*     2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after
-*        TM-score superposition.
-*************************************************************************
-*************************************************************************
-      subroutine TMscore(L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,TM,Rcomm,Lcomm)
-      include 'DIMENSIONS'
-      PARAMETER(nmax=maxres)
-      common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
-      common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
-      common/para/d,d0
-      common/align/n_ali,iA(nmax),iB(nmax)
-      common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
-      dimension k_ali(nmax),k_ali0(nmax)
-      dimension L_ini(100),iq(nmax)
-      common/scores/score
-      double precision score,score_max
-      dimension xa(nmax),ya(nmax),za(nmax)
-
-      dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
-      dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
-
-ccc   RMSD:
-      double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
-      double precision u(3,3),t(3),rms,drms !armsd is real
-      data w /nmax*1.0/
-ccc   
-
-********* convert input data ****************
-      nseqA=L1
-      do i=1,nseqA
-         xa(i)=x1(i)
-         ya(i)=y1(i)
-         za(i)=z1(i)
-         nresA(i)=n1(i)
-      enddo
-      nseqB=L2
-      do i=1,L2
-         xb(i)=x2(i)
-         yb(i)=y2(i)
-         zb(i)=z2(i)
-         nresB(i)=n2(i)
-      enddo
-
-******************************************************************
-*     pickup the aligned residues:
-******************************************************************
-      k=0
-      do i=1,nseqA
-         do j=1,nseqB
-            if(nresA(i).eq.nresB(j))then
-               k=k+1
-               iA(k)=i
-               iB(k)=j
-               goto 205
-            endif
-         enddo
- 205     continue
-      enddo
-      n_ali=k                   !number of aligned residues
-      Lcomm=n_ali
-      if(n_ali.lt.1)then
-c         write(*,*)'There is no common residues in the input structures'
-         TM=0
-         Rcomm=0
-         return
-      endif
-
-************/////
-*     parameters:
-*****************
-***   d0------------->
-      d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
-      if(d0.lt.0.5)d0=0.5
-***   d0_search ----->
-      d0_search=d0
-      if(d0_search.gt.8)d0_search=8
-      if(d0_search.lt.4.5)d0_search=4.5
-***   iterative parameters ----->
-      n_it=20                   !maximum number of iterations
-      d_output=5                !for output alignment
-      n_init_max=6              !maximum number of L_init
-      n_init=0
-      L_ini_min=4
-      if(n_ali.lt.4)L_ini_min=n_ali
-      do i=1,n_init_max-1
-         n_init=n_init+1
-         L_ini(n_init)=n_ali/2**(n_init-1)
-         if(L_ini(n_init).le.L_ini_min)then
-            L_ini(n_init)=L_ini_min
-            goto 402
-         endif
-      enddo
-      n_init=n_init+1
-      L_ini(n_init)=L_ini_min
- 402  continue
-
-******************************************************************
-*     find the maximum score starting from local structures superposition
-******************************************************************
-      score_max=-1              !TM-score
-      do 333 i_init=1,n_init
-        L_init=L_ini(i_init)
-        iL_max=n_ali-L_init+1
-        do 300 iL=1,iL_max      !on aligned residues, [1,nseqA]
-           LL=0
-           ka=0
-           do i=1,L_init
-              k=iL+i-1          ![1,n_ali] common aligned
-              r_1(1,i)=xa(iA(k))
-              r_1(2,i)=ya(iA(k))
-              r_1(3,i)=za(iA(k))
-              r_2(1,i)=xb(iB(k))
-              r_2(2,i)=yb(iB(k))
-              r_2(3,i)=zb(iB(k))
-              LL=LL+1
-              ka=ka+1
-              k_ali(ka)=k
-           enddo
-           call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
-           if(i_init.eq.1)then  !global superposition
-              armsd=dsqrt(rms/LL)
-              Rcomm=armsd
-           endif
-           do j=1,nseqA
-              xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
-              yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
-              zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
-           enddo
-           d=d0_search-1
-           call score_fun       !init, get scores, n_cut+i_ali(i) for iteration
-           if(score_max.lt.score)then
-              score_max=score
-              ka0=ka
-              do i=1,ka0
-                 k_ali0(i)=k_ali(i)
-              enddo
-           endif
-***   iteration for extending ---------------------------------->
-           d=d0_search+1
-           do 301 it=1,n_it
-              LL=0
-              ka=0
-              do i=1,n_cut
-                 m=i_ali(i)     ![1,n_ali]
-                 r_1(1,i)=xa(iA(m))
-                 r_1(2,i)=ya(iA(m))
-                 r_1(3,i)=za(iA(m))
-                 r_2(1,i)=xb(iB(m))
-                 r_2(2,i)=yb(iB(m))
-                 r_2(3,i)=zb(iB(m))
-                 ka=ka+1
-                 k_ali(ka)=m
-                 LL=LL+1
-              enddo
-              call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
-              do j=1,nseqA
-                 xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
-                 yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
-                 zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
-              enddo
-              call score_fun    !get scores, n_cut+i_ali(i) for iteration
-              if(score_max.lt.score)then
-                 score_max=score
-                 ka0=ka
-                 do i=1,ka
-                    k_ali0(i)=k_ali(i)
-                 enddo
-              endif
-              if(it.eq.n_it)goto 302
-              if(n_cut.eq.ka)then
-                 neq=0
-                 do i=1,n_cut
-                    if(i_ali(i).eq.k_ali(i))neq=neq+1
-                 enddo
-                 if(n_cut.eq.neq)goto 302
-              endif
- 301       continue             !for iteration
- 302       continue
- 300    continue                !for shift
- 333  continue                  !for initial length, L_ali/M
-
-******** return the final rotation ****************
-      LL=0
-      do i=1,ka0
-         m=k_ali0(i)            !record of the best alignment
-         r_1(1,i)=xa(iA(m))
-         r_1(2,i)=ya(iA(m))
-         r_1(3,i)=za(iA(m))
-         r_2(1,i)=xb(iB(m))
-         r_2(2,i)=yb(iB(m))
-         r_2(3,i)=zb(iB(m))
-         LL=LL+1
-      enddo
-      call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
-      do j=1,nseqA
-         x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
-         y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
-         z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
-      enddo
-      TM=score_max
-
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-      return
-      END
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c     1, collect those residues with dis<d;
-c     2, calculate score_GDT, score_maxsub, score_TM
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      subroutine score_fun
-      include 'DIMENSIONS'
-      PARAMETER(nmax=maxres)
-
-      common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
-      common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
-      common/para/d,d0
-      common/align/n_ali,iA(nmax),iB(nmax)
-      common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
-      common/scores/score
-      double precision score,score_max
-
-      d_tmp=d
- 21   n_cut=0                   !number of residue-pairs dis<d, for iteration
-      score_sum=0               !TMscore
-      do k=1,n_ali
-         i=iA(k)                ![1,nseqA] reoder number of structureA
-         j=iB(k)                ![1,nseqB]
-         dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)
-         if(dis.lt.d_tmp)then
-            n_cut=n_cut+1
-            i_ali(n_cut)=k      ![1,n_ali], mark the residue-pairs in dis<d
-         endif
-         score_sum=score_sum+1/(1+(dis/d0)**2)
-      enddo
-      if(n_cut.lt.3.and.n_ali.gt.3)then
-         d_tmp=d_tmp+.5
-         goto 21
-      endif
-      score=score_sum/float(nseqB) !TM-score
-      
-      return
-      end
-
-cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
-c  w    - w(m) is weight for atom pair  c m           (given)
-c  x    - x(i,m) are coordinates of atom c m in set x       (given)
-c  y    - y(i,m) are coordinates of atom c m in set y       (given)
-c  n    - n is number of atom pairs                         (given)
-c  mode  - 0:calculate rms only                             (given)
-c          1:calculate rms,u,t                              (takes longer)
-c  rms   - sum of w*(ux+t-y)**2 over all atom pairs         (result)
-c  u    - u(i,j) is   rotation  matrix for best superposition  (result)
-c  t    - t(i)   is translation vector for best superposition  (result)
-c  ier  - 0: a unique optimal superposition has been determined(result)
-c       -1: superposition is not unique but optimal
-c       -2: no result obtained because of negative weights w
-c           or all weights equal to zero.
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-      subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
-      double precision w(*), x(3,*), y(3,*)
-      integer n, mode
-      
-      double precision rms, u(3,3), t(3)
-      integer ier
-      
-      integer i, j, k, l, m1, m
-      integer ip(9), ip2312(4)
-      double precision r(3,3), xc(3), yc(3), wc
-      double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
-      double precision e0, d, spur, det, cof, h, g
-      double precision cth, sth, sqrth, p, sigma
-      
-      double precision sqrt3, tol, zero
-      
-      data sqrt3 / 1.73205080756888d+00 /
-      data tol / 1.0d-2 /
-      data zero / 0.0d+00 /
-      data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
-      data ip2312 / 2, 3, 1, 2 /
-      
-      wc = zero
-      rms = zero
-      e0 = zero
-      
-      do i=1, 3
-         xc(i) = zero
-         yc(i) = zero
-         t(i) = zero
-         do j=1, 3
-            r(i,j) = zero
-            u(i,j) = zero
-            a(i,j) = zero
-            if( i .eq. j ) then
-               u(i,j) = 1.0
-               a(i,j) = 1.0
-            end if
-         end do
-      end do
-      
-      ier = -1
-      if( n .lt. 1 ) return
-      ier = -2
-      do m=1, n
-         if( w(m) .lt. 0.0 ) return
-         wc = wc + w(m)
-         do i=1, 3
-            xc(i) = xc(i) + w(m)*x(i,m)
-            yc(i) = yc(i) + w(m)*y(i,m)
-         end do
-      end do
-      if( wc .le. zero ) return
-      do i=1, 3
-         xc(i) = xc(i) / wc
-         yc(i) = yc(i) / wc
-      end do
-      
-      do m=1, n
-         do i=1, 3
-            e0=e0+w(m)*((x(i,m)-xc(i))**2+(y(i,m)-yc(i))**2)
-            d = w(m) * ( y(i,m) - yc(i) )
-            do j=1, 3
-               r(i,j) = r(i,j) + d*( x(j,m) - xc(j) )
-            end do
-         end do
-      end do
-      
-      det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
-     &     - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
-     &     + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
-      
-      sigma = det
-      
-      m = 0;
-      do j=1, 3
-         do i=1, j
-            m = m+1
-            rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
-         end do
-      end do
-      
-      spur = (rr(1)+rr(3)+rr(6)) / 3.0
-      cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
-     &     - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
-      det = det*det
-      
-      do i=1, 3
-         e(i) = spur
-      end do
-      if( spur .le. zero ) goto 40
-      d = spur*spur
-      h = d - cof
-      g = (spur*cof - det)/2.0 - spur*h
-      if( h .le. zero ) then
-         if( mode .eq. 0 ) then
-            goto 50
-         else
-            goto 30
-         end if
-      end if
-      sqrth = dsqrt(h)
-      d = h*h*h - g*g
-      if( d .lt. zero ) d = zero
-      d = datan2( dsqrt(d), -g ) / 3.0
-      cth = sqrth * dcos(d)
-      sth = sqrth*sqrt3*dsin(d)
-      e(1) = (spur + cth) + cth
-      e(2) = (spur - cth) + sth
-      e(3) = (spur - cth) - sth
-	
-      if( mode .eq. 0 ) then
-         goto 50
-      end if
-      
-      do l=1, 3, 2
-         d = e(l)
-         ss(1) = (d-rr(3)) * (d-rr(6))  - rr(5)*rr(5)
-         ss(2) = (d-rr(6)) * rr(2)      + rr(4)*rr(5)
-         ss(3) = (d-rr(1)) * (d-rr(6))  - rr(4)*rr(4)
-         ss(4) = (d-rr(3)) * rr(4)      + rr(2)*rr(5)
-         ss(5) = (d-rr(1)) * rr(5)      + rr(2)*rr(4)
-         ss(6) = (d-rr(1)) * (d-rr(3))  - rr(2)*rr(2)
-         
-         if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
-            j=1
-            if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
-         else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
-            j = 2
-         else
-            j = 3
-         end if
-         
-         d = zero
-         j = 3 * (j - 1)
-         
-         do i=1, 3
-            k = ip(i+j)
-            a(i,l) = ss(k)
-            d = d + ss(k)*ss(k)
-         end do
-         if( d .gt. zero ) d = 1.0 / dsqrt(d)
-         do i=1, 3
-            a(i,l) = a(i,l) * d
-         end do
-      end do
-      
-      d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
-      if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
-         m1 = 3
-         m = 1
-      else
-         m1 = 1
-         m = 3
-      endif
-      
-      p = zero
-      do i=1, 3
-         a(i,m1) = a(i,m1) - d*a(i,m)
-         p = p + a(i,m1)**2
-      end do
-      if( p .le. tol ) then
-         p = 1.0
-         do i=1, 3
-            if (p .lt. dabs(a(i,m))) cycle
-            p = dabs( a(i,m) )
-            j = i
-         end do
-         k = ip2312(j)
-         l = ip2312(j+1)
-         p = dsqrt( a(k,m)**2 + a(l,m)**2 )
-         if( p .le. tol ) goto 40
-         a(j,m1) = zero
-         a(k,m1) = -a(l,m)/p
-         a(l,m1) =  a(k,m)/p
-      else
-         p = 1.0 / dsqrt(p)
-         do i=1, 3
-            a(i,m1) = a(i,m1)*p
-         end do
-      end if
-      
-      a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
-      a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
-      a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
-      
- 30   do l=1, 2
-         d = zero
-         do i=1, 3
-            b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
-            d = d + b(i,l)**2
-         end do
-         if( d .gt. zero ) d = 1.0 / dsqrt(d)
-         do i=1, 3
-            b(i,l) = b(i,l)*d
-         end do
-      end do
-      d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
-      p = zero
-      
-      do i=1, 3
-         b(i,2) = b(i,2) - d*b(i,1)
-         p = p + b(i,2)**2
-      end do
-      if( p .le. tol ) then
-         p = 1.0
-         do i=1, 3
-            if( p .lt. dabs(b(i,1)) ) cycle
-            p = dabs( b(i,1) )
-            j = i
-         end do
-         k = ip2312(j)
-         l = ip2312(j+1)
-         p = dsqrt( b(k,1)**2 + b(l,1)**2 )
-         if( p .le. tol ) goto 40
-         b(j,2) = zero
-         b(k,2) = -b(l,1)/p
-         b(l,2) =  b(k,1)/p
-      else
-         p = 1.0 / dsqrt(p)
-         do i=1, 3
-            b(i,2) = b(i,2)*p
-         end do
-      end if
-      
-      b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
-      b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
-      b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
-      
-      do i=1, 3
-         do j=1, 3
-            u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
-         end do
-      end do
-      
- 40   do i=1, 3
-         t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
-      end do
- 50   do i=1, 3
-         if( e(i) .lt. zero ) e(i) = zero
-         e(i) = dsqrt( e(i) )
-      end do
-      
-      ier = 0
-      if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
-      
-      d = e(3)
-      if( sigma .lt. 0.0 ) then
-         d = - d
-         if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
-      end if
-      d = (d + e(2)) + e(1)
-      
-      rms = (e0 - d) - d
-      if( rms .lt. 0.0 ) rms = 0.0
-      
-      return
-      end