source/cluster/wham/src-HCD-5D/TMscore.F

   1 *************************************************************************
   2 *     This program is to compare two protein structures and identify the
   3 *     best superposition that has the highest TM-score. Input structures
   4 *     must be in the PDB format. By default, TM-score is normalized by
   5 *     the second protein. Users can obtain a brief instruction by simply
   6 *     running the program without arguments. For comments/suggestions,
   7 *     please contact email: zhng@umich.edu.
   8 *
   9 *     Reference:
  10 *     Yang Zhang, Jeffrey Skolnick, Proteins, 2004 57:702-10.
  11 *
  12 *     Permission to use, copy, modify, and distribute this program for
  13 *     any purpose, with or without fee, is hereby granted, provided that
  14 *     the notices on the head, the reference information, and this
  15 *     copyright notice appear in all copies or substantial portions of
  16 *     the Software. It is provided "as is" without express or implied
  17 *     warranty.
  18 ******************* Updating history ************************************
  19 *     2005/10/19: the program was reformed so that the score values.
  20 *                 are not dependent on the specific compilers.
  21 *     2006/06/20: selected 'A' if there is altLoc when reading PDB file.
  22 *     2007/02/05: fixed a bug with length<15 in TMscore_32.
  23 *     2007/02/27: rotation matrix from Chain-1 to Chain-2 was added.
  24 *     2007/12/06: GDT-HA score was added, fixed a bug for reading PDB.
  25 *     2010/08/02: A new RMSD matrix was used and obsolete statement removed.
  26 *     2011/01/03: The length of pdb file names were extended to 500.
  27 *     2011/01/30: An open source license is attached to the program.
  28 *     2012/05/07: Improved RMSD calculation subroutine which speeds up
  29 *                 TM-score program by 30%.
  30 *     2012/06/05: Added option '-l L' which calculates TM-score (and maxsub
  31 *                 and GDT scores) normalized by a specific length 'L'.
  32 *************************************************************************
  33
  34 c      program TMscore
  35       subroutine TMscore_sub(rmsd,gdt_ts,gdt_ha,tmscore,cfname,lprint)
  36       include 'DIMENSIONS'
  37       PARAMETER(nmax=5000)
  38       include 'COMMON.CONTROL'
  39       include 'COMMON.IOUNITS'
  40       include 'COMMON.CHAIN'
  41       include 'COMMON.INTERACT'
  42
  43       real*8 rmsd,gdt_ts,gdt_ha,tmscore
  44       common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
  45       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
  46       common/para/d,d0,d0_fix
  47       common/align/n_ali,iA(nmax),iB(nmax)
  48       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
  49       dimension k_ali(nmax),k_ali0(nmax)
  50
  51       character*500 fnam,pdb(100)!,outname
  52       character*80 cfname
  53       character*3 aa(-1:20),seqA(nmax),seqB(nmax)
  54       character*500 s,du
  55       character seq1A(nmax),seq1B(nmax),ali(nmax)
  56       character sequenceA(nmax),sequenceB(nmax),sequenceM(nmax)
  57
  58       dimension L_ini(100),iq(nmax)
  59       common/scores/score,score_maxsub,score_fix,score10
  60       common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
  61       double precision score,score_max,score_fix,score_fix_max
  62       double precision score_maxsub,score10
  63       dimension xa(nmax),ya(nmax),za(nmax)
  64
  65 ccc   RMSD:
  66       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
  67       double precision u(3,3),tt(3),rms,drms !armsd is real
  68       data w /nmax*1.0/
  69       integer ii,ipermmin,iperm
  70
  71       logical lprint
  72 ccc
  73
  74       data aa/ 'BCK','GLY','ALA','SER','CYS',
  75      &     'VAL','THR','ILE','PRO','MET',
  76      &     'ASP','ASN','LEU','LYS','GLU',
  77      &     'GLN','ARG','HIS','PHE','TYR',
  78      &     'TRP','CYX'/
  79       character*1 slc(-1:20)
  80       data slc/'X','G','A','S','C',
  81      &     'V','T','I','P','M',
  82      &     'D','N','L','K','E',
  83      &     'Q','R','H','F','Y',
  84      &     'W','C'/
  85
  86 *****instructions ----------------->
  87 c      call getarg(1,fnam)
  88 c      if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then
  89 c         write(*,*)
  90 c         write(*,*)'Brief instruction for running TM-score program:'
  91 c         write(*,*)'(For detail: Zhang & Skolnick,  Proteins, 2004',
  92 c     &        ' 57:702-10)'
  93 c         write(*,*)
  94 c         write(*,*)'1. Run TM-score to compare ''model'' and ',
  95 c     &        '''native'':'
  96 c         write(*,*)'   >TMscore model native'
  97 c         write(*,*)
  98 c         write(*,*)'2. TM-score normalized with an assigned scale d0',
  99 c     &        ' e.g. 5 A:'
 100 c         write(*,*)'   >TMscore model native -d 5'
 101 c         write(*,*)
 102 c         write(*,*)'3. TM-score normalized by a specific length, ',
 103 c     &        'e.g. 120 AA:'
 104 c         write(*,*)'   >TMscore model native -l 120'
 105 c         write(*,*)
 106 c         write(*,*)'4. TM-score with superposition output, e.g. ',
 107 c     &        '''TM.sup'':'
 108 c         write(*,*)'   >TMscore model native -o TM.sup'
 109 c         write(*,*)'   To view the superimposed structures by rasmol:'
 110 c         write(*,*)'      >rasmol -script TM.sup'
 111 c         write(*,*)
 112 c         goto 9999
 113 c      endif
 114
 115       pdb(1)=cfname
 116       pdb(2)=pdbfile
 117 ******* options ----------->
 118       m_out=-1
 119       m_fix=-1
 120       m_len=-1
 121 c      narg=iargc()
 122 c      i=0
 123 c      j=0
 124 c 115  continue
 125 c      i=i+1
 126 c      call getarg(i,fnam)
 127 c      if(fnam.eq.'-o')then
 128 c         m_out=1
 129 c         i=i+1
 130 c         call getarg(i,outname)
 131 c      elseif(fnam.eq.'-d')then
 132 c         m_fix=1
 133 c         i=i+1
 134 c         call getarg(i,fnam)
 135 c         read(fnam,*)d0_fix
 136 c      elseif(fnam.eq.'-l')then
 137 c         m_len=1
 138 c         i=i+1
 139 c         call getarg(i,fnam)
 140 c         read(fnam,*)l0_fix
 141 c      else
 142 c         j=j+1
 143 c         pdb(j)=fnam
 144 c      endif
 145 c      if(i.lt.narg)goto 115
 146 c
 147 ccccccccc read data from first CA file:
 148 c      open(unit=10,file=pdb(1),status='old')
 149 c      i=0
 150 c 101  read(10,104,end=102) s
 151 c      if(s(1:3).eq.'TER') goto 102
 152 c      if(s(1:4).eq.'ATOM')then
 153 c         if(s(13:16).eq.'CA  '.or.s(13:16).eq.' CA '.or.s(13:16).
 154 c     &        eq.'  CA')then
 155 c         if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
 156 c            i=i+1
 157 c            read(s,103)du,seqA(i),du,nresA(i),du,xa(i),ya(i),za(i)
 158 c            do j=-1,20
 159 c               if(seqA(i).eq.aa(j))then
 160 c                  seq1A(i)=slc(j)
 161 c                  goto 21
 162 c               endif
 163 c            enddo
 164 c            seq1A(i)=slc(-1)
 165 c 21         continue
 166 c         endif
 167 c         endif
 168 c      endif
 169 c      goto 101
 170 c 102  continue
 171 c 103  format(A17,A3,A2,i4,A4,3F8.3)
 172 c 104  format(A100)
 173 c      close(10)
 174 c      nseqA=i
 175 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 176 c
 177 ccccccccc read data from first CA file:
 178 c      open(unit=10,file=pdb(2),status='old')
 179 c      i=0
 180 c 201  read(10,204,end=202) s
 181 c      if(s(1:3).eq.'TER') goto 202
 182 c      if(s(1:4).eq.'ATOM')then
 183 c         if(s(13:16).eq.'CA  '.or.s(13:16).eq.' CA '.or.s(13:16).
 184 c     &        eq.'  CA')then
 185 c         if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
 186 c            i=i+1
 187 c            read(s,203)du,seqB(i),du,nresB(i),du,xb(i),yb(i),zb(i)
 188 c            do j=-1,20
 189 c               if(seqB(i).eq.aa(j))then
 190 c                  seq1B(i)=slc(j)
 191 c                  goto 22
 192 c               endif
 193 c            enddo
 194 c            seq1B(i)=slc(-1)
 195 c 22         continue
 196 c         endif
 197 c         endif
 198 c      endif
 199 c      goto 201
 200 c 202  continue
 201 c 203  format(A17,A3,A2,i4,A4,3F8.3)
 202 c 204  format(A100)
 203 c      close(10)
 204 c      nseqB=i
 205 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 206
 207 ******************************************************************
 208 *     pickup the aligned residues:
 209 ******************************************************************
 210 c      k=0
 211 c      do i=1,nseqA
 212 c         do j=1,nseqB
 213 c            if(nresA(i).eq.nresB(j))then
 214 c               k=k+1
 215 c               iA(k)=i
 216 c               iB(k)=j
 217 c               goto 205
 218 c            endif
 219 c         enddo
 220 c 205     continue
 221 c      enddo
 222 c      n_ali=k                   !number of aligned residues
 223 c      if(n_ali.lt.1)then
 224 c        write(*,*)'There is no common residues in the input structures'
 225 c        goto 9999
 226 c      endif
 227 c
 228 ************/////
 229 *     parameters:
 230 *****************
 231
 232       DO II=1,NPERMCHAIN
 233
 234       noverlap=nres
 235       if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
 236       nnsup=0
 237       do i=1,noverlap
 238         if (itype(i).ne.ntyp1) then
 239           nnsup=nnsup+1
 240           iA(nnsup)=nnsup
 241           iB(nnsup)=nnsup
 242         endif
 243       enddo
 244       nseqA=nnsup
 245       nseqB=nnsup
 246       n_ali=nnsup
 247 ***   d0------------->
 248       if(nseqB.gt.15)then
 249          d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
 250       else
 251          d0=0.5
 252       endif
 253       if(m_len.eq.1)then
 254          d0=1.24*(l0_fix-15)**(1.0/3.0)-1.8
 255       endif
 256       if(d0.lt.0.5)d0=0.5
 257       if(m_fix.eq.1)d0=d0_fix
 258 ***   d0_search ----->
 259       d0_search=d0
 260       if(d0_search.gt.8)d0_search=8
 261       if(d0_search.lt.4.5)d0_search=4.5
 262 ***   iterative parameters ----->
 263       n_it=20                   !maximum number of iterations
 264       d_output=5                !for output alignment
 265       if(m_fix.eq.1)d_output=d0_fix
 266       n_init_max=6              !maximum number of L_init
 267       n_init=0
 268       L_ini_min=4
 269       if(n_ali.lt.4)L_ini_min=n_ali
 270       do i=1,n_init_max-1
 271          n_init=n_init+1
 272          L_ini(n_init)=n_ali/2**(n_init-1)
 273          if(L_ini(n_init).le.L_ini_min)then
 274             L_ini(n_init)=L_ini_min
 275             goto 402
 276          endif
 277       enddo
 278       n_init=n_init+1
 279       L_ini(n_init)=L_ini_min
 280  402  continue
 281
 282 ******************************************************************
 283 *     find the maximum score starting from local structures superposition
 284 ******************************************************************
 285       score_max=-1              !TM-score
 286       score_maxsub_max=-1       !MaxSub-score
 287       score10_max=-1            !TM-score10
 288       n_GDT05_max=-1            !number of residues<0.5
 289       n_GDT1_max=-1             !number of residues<1
 290       n_GDT2_max=-1             !number of residues<2
 291       n_GDT4_max=-1             !number of residues<4
 292       n_GDT8_max=-1             !number of residues<8
 293
 294 #ifdef DEBUG
 295       write (iout,*) "cref and ccref"
 296 #endif
 297       noverlap=nres
 298       if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
 299       nnsup=0
 300       do i=1,noverlap
 301         if (itype(i).ne.ntyp1) then
 302           nnsup=nnsup+1
 303           xa(nnsup)=c(1,iperm(i,ii))
 304           ya(nnsup)=c(2,iperm(i,ii))
 305           za(nnsup)=c(3,iperm(i,ii))
 306           xb(nnsup)=cref_pdb(1,i)
 307           yb(nnsup)=cref_pdb(2,i)
 308           zb(nnsup)=cref_pdb(3,i)
 309 c          do j=1,3
 310 c            cc(j,nnsup)=c(j,i)
 311 c            ccref(j,nnsup)=cref_pdb(j,i,1)
 312 c          enddo
 313 #ifdef DEBUG
 314           write (iout,'(i5,3f10.5,5x,3f10.5)') nnsup,
 315      &     xa(nnsup),ya(nnsup),za(nnsup),xb(nnsup),yb(nnsup),zb(nnsup)
 316 #endif
 317         endif
 318       enddo
 319
 320       do 333 i_init=1,n_init
 321         L_init=L_ini(i_init)
 322         iL_max=n_ali-L_init+1
 323         do 300 iL=1,iL_max      !on aligned residues, [1,nseqA]
 324            LL=0
 325            ka=0
 326            do i=1,L_init
 327               k=iL+i-1          ![1,n_ali] common aligned
 328               r_1(1,i)=xa(iA(k))
 329               r_1(2,i)=ya(iA(k))
 330               r_1(3,i)=za(iA(k))
 331               r_2(1,i)=xb(iB(k))
 332               r_2(2,i)=yb(iB(k))
 333               r_2(3,i)=zb(iB(k))
 334               ka=ka+1
 335               k_ali(ka)=k
 336               LL=LL+1
 337            enddo
 338            if(i_init.eq.1)then  !global superposition
 339               call u3b(w,r_1,r_2,LL,2,rms,u,tt,ier) !0:rmsd; 1:u,t; 2:rmsd,u,t
 340               armsd=dsqrt(rms/LL)
 341               rmsd_ali=armsd
 342            else
 343               call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
 344            endif
 345            do j=1,nseqA
 346               xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 347               yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 348               zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 349            enddo
 350            d=d0_search-1
 351            call score_fun       !init, get scores, n_cut+i_ali(i) for iteration
 352            if(score_max.lt.score)then
 353               score_max=score
 354               ka0=ka
 355               do i=1,ka0
 356                  k_ali0(i)=k_ali(i)
 357               enddo
 358            endif
 359            if(score10_max.lt.score10)score10_max=score10
 360            if(score_maxsub_max.lt.score_maxsub)score_maxsub_max=
 361      &          score_maxsub
 362            if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
 363            if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
 364            if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
 365            if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
 366            if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
 367 ***   iteration for extending ---------------------------------->
 368            d=d0_search+1
 369            do 301 it=1,n_it
 370               LL=0
 371               ka=0
 372               do i=1,n_cut
 373                  m=i_ali(i)     ![1,n_ali]
 374                  r_1(1,i)=xa(iA(m))
 375                  r_1(2,i)=ya(iA(m))
 376                  r_1(3,i)=za(iA(m))
 377                  r_2(1,i)=xb(iB(m))
 378                  r_2(2,i)=yb(iB(m))
 379                  r_2(3,i)=zb(iB(m))
 380                  ka=ka+1
 381                  k_ali(ka)=m
 382                  LL=LL+1
 383               enddo
 384               call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
 385               do j=1,nseqA
 386                  xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 387                  yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 388                  zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 389               enddo
 390               call score_fun    !get scores, n_cut+i_ali(i) for iteration
 391               if(score_max.lt.score)then
 392                  score_max=score
 393                  ka0=ka
 394                  do i=1,ka
 395                     k_ali0(i)=k_ali(i)
 396                  enddo
 397               endif
 398               if(score10_max.lt.score10)score10_max=score10
 399               if(score_maxsub_max.lt.score_maxsub)score_maxsub_max
 400      &             =score_maxsub
 401               if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
 402               if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
 403               if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
 404               if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
 405               if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
 406               if(it.eq.n_it)goto 302
 407               if(n_cut.eq.ka)then
 408                  neq=0
 409                  do i=1,n_cut
 410                     if(i_ali(i).eq.k_ali(i))neq=neq+1
 411                  enddo
 412                  if(n_cut.eq.neq)goto 302
 413               endif
 414  301       continue             !for iteration
 415  302       continue
 416  300    continue                !for shift
 417  333  continue                  !for initial length, L_ali/M
 418 c
 419       ratio=1
 420       if(m_len.gt.0)then
 421          ratio=float(nseqB)/float(l0_fix)
 422       endif
 423       if(m_len.eq.1)then
 424          score_max=score_max*float(nseqB)/float(l0_fix)
 425       endif
 426       score_GDT=(n_GDT1_max+n_GDT2_max+n_GDT4_max+n_GDT8_max)
 427      &     /float(4*nseqB)
 428       score_GDT_HA=(n_GDT05_max+n_GDT1_max+n_GDT2_max+n_GDT4_max)
 429      &     /float(4*nseqB)
 430       tmscore=score_max
 431       gdt_ts=score_GDT*ratio
 432       gdt_ha=score_GDT_HA*ratio
 433       rmsd=rmsd_ali
 434
 435       if (ii.eq.1 .or. rmsd.lt.rmsd_min) then
 436         rmsd_min=rmsd
 437         tmscore_min=tmscore
 438         gdt_ts_min=gdt_ts
 439         gdt_ha_min=gdt_ha
 440         ipermmin=ii
 441       endif
 442
 443       ENDDO
 444
 445       rmsd=rmsd_min
 446       tmscore=tmscore_min
 447       gdt_ts=gdt_ts_min
 448       gdt_ha=gdt_ha_min
 449
 450 ******************************************************************
 451 *     Output
 452 ******************************************************************
 453 ***   output TM-scale ---------------------------->
 454
 455       if (lprint) then
 456
 457       write(iout,*)
 458       write(iout,*)'**************************************************',
 459      &     '***************************'
 460       write(iout,*)'*                              TM-SCORE           ',
 461      &     '                          *'
 462       write(iout,*)'* A scoring function to assess the similarity of p',
 463      &     'rotein structures         *'
 464       write(iout,*)'* Based on statistics:                            ',
 465      &     '                          *'
 466       write(iout,*)'*       0.0 < TM-score < 0.17, random structural s',
 467      &     'imilarity                 *'
 468       write(iout,*)'*       0.5 < TM-score < 1.00, in about the same f',
 469      &     'old                    *'
 470       write(iout,*)'* Reference: Yang Zhang and Jeffrey Skolnick, ',
 471      &     'Proteins 2004 57: 702-710     *'
 472       write(iout,*)'* For comments, please email to: zhng@umich.edu   ',
 473      &     '                          *'
 474       write(iout,*)'**************************************************',
 475      &     '***************************'
 476       write(iout,*)
 477       write(iout,501)pdb(1),nseqA
 478  501  format('Structure1: ',A10,'  Length= ',I4)
 479       if(m_len.eq.1)then
 480          write(iout,411)pdb(2),nseqB
 481          write(iout,412)l0_fix
 482       else
 483          write(iout,502)pdb(2),nseqB
 484       endif
 485  411  format('Structure2: ',A10,'  Length= ',I4)
 486  412  format('TM-score is notmalized by ',I4)
 487  502  format('Structure2: ',A10,'  Length= ',I4,
 488      &     ' (by which all scores are normalized)')
 489       write(iout,503)n_ali
 490  503  format('Number of residues in common= ',I4)
 491       write(iout,513)rmsd_ali
 492  513  format('RMSD of  the common residues= ',F8.3)
 493       write(iout,*)
 494       write(iout,504)score_max,d0
 495  504  format('TM-score    = ',f6.4,'  (d0=',f5.2,')')
 496       write(iout,505)score_maxsub_max*ratio
 497  505  format('MaxSub-score= ',f6.4,'  (d0= 3.50)')
 498       write(iout,506)score_GDT*ratio,n_GDT1_max/float(nseqB)*ratio,
 499      &     n_GDT2_max/float(nseqB)*ratio,n_GDT4_max/float(nseqB)*ratio,
 500      &     n_GDT8_max/float(nseqB)*ratio
 501  506  format('GDT-TS-score= ',f6.4,' %(d<1)=',f6.4,' %(d<2)=',f6.4,
 502      $     ' %(d<4)=',f6.4,' %(d<8)=',f6.4)
 503       write(iout,507)score_GDT_HA*ratio,n_GDT05_max/float(nseqB)*ratio,
 504      &     n_GDT1_max/float(nseqB)*ratio,n_GDT2_max/float(nseqB)*ratio,
 505      &     n_GDT4_max/float(nseqB)*ratio
 506  507  format('GDT-HA-score= ',f6.4,' %(d<0.5)=',f6.4,' %(d<1)=',f6.4,
 507      $     ' %(d<2)=',f6.4,' %(d<4)=',f6.4)
 508       write (iout,*) "Permutation",ipermmin
 509       write(iout,*)
 510
 511       endif
 512
 513       return
 514       end
 515 c------------------------------------------------------------------------
 516 ***   recall and output the superposition of maxiumum TM-score:
 517 c      LL=0
 518 c      do i=1,ka0
 519 c         m=k_ali0(i)            !record of the best alignment
 520 c         r_1(1,i)=xa(iA(m))
 521 c         r_1(2,i)=ya(iA(m))
 522 c         r_1(3,i)=za(iA(m))
 523 c         r_2(1,i)=xb(iB(m))
 524 c         r_2(2,i)=yb(iB(m))
 525 c         r_2(3,i)=zb(iB(m))
 526 c         LL=LL+1
 527 c      enddo
 528 c      call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
 529 c      do j=1,nseqA
 530 c         xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 531 c         yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 532 c         zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 533 c      enddo
 534 c
 535 c********* extract rotation matrix ------------>
 536 c      write(*,*)'-------- rotation matrix to rotate Chain-1 to ',
 537 c     &     'Chain-2 ------'
 538 c      write(*,*)'i          t(i)         u(i,1)         u(i,2) ',
 539 c     &     '        u(i,3)'
 540 c      do i=1,3
 541 c         write(*,304)i,t(i),u(i,1),u(i,2),u(i,3)
 542 c      enddo
 543 cc      do j=1,nseqA
 544 cc         xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
 545 cc         yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
 546 cc         zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
 547 cc         write(*,*)j,xt(j),yt(j),zt(j)
 548 cc      enddo
 549 c      write(*,*)
 550 c 304  format(I2,f18.10,f15.10,f15.10,f15.10)
 551 c
 552 c********* rmsd in superposed regions --------------->
 553 c      d=d_output                !for output
 554 c      call score_fun()          !give i_ali(i), score_max=score now
 555 c      LL=0
 556 c      do i=1,n_cut
 557 c         m=i_ali(i)             ![1,nseqA]
 558 c         r_1(1,i)=xa(iA(m))
 559 c         r_1(2,i)=ya(iA(m))
 560 c         r_1(3,i)=za(iA(m))
 561 c         r_2(1,i)=xb(iB(m))
 562 c         r_2(2,i)=yb(iB(m))
 563 c         r_2(3,i)=zb(iB(m))
 564 c         LL=LL+1
 565 c      enddo
 566 c      call u3b(w,r_1,r_2,LL,0,rms,u,t,ier)
 567 c      armsd=dsqrt(rms/LL)
 568 c      rmsd=armsd
 569 c
 570 c***   output rotated chain1 + chain2----->
 571 c      if(m_out.ne.1)goto 999
 572 c      OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
 573 c 900  format(A)
 574 c 901  format('select ',I4)
 575 c      write(7,900)'load inline'
 576 c      write(7,900)'select atomno<1000'
 577 cc      write(7,900)'color [255,20,147]'
 578 c      write(7,900)'wireframe .45'
 579 c      write(7,900)'select none'
 580 c      write(7,900)'select atomno>1000'
 581 cc      write(7,900)'color [100,149,237]'
 582 c      write(7,900)'wireframe .15'
 583 c      write(7,900)'color white'
 584 c      do i=1,n_cut
 585 c         write(7,901)nresA(iA(i_ali(i)))
 586 c         write(7,900)'color red'
 587 c      enddo
 588 c      write(7,900)'select all'
 589 c      write(7,900)'exit'
 590 c      write(7,514)rmsd_ali
 591 c 514  format('REMARK  RMSD of the common residues=',F8.3)
 592 c      write(7,515)score_max,d0
 593 c 515  format('REMARK  TM-score=',f6.4,' (d0=',f5.2,')')
 594 c      do i=1,nseqA
 595 c         write(7,1237)nresA(i),seqA(i),nresA(i),xt(i),yt(i),zt(i)
 596 c      enddo
 597 c      write(7,1238)
 598 c      do i=2,nseqA
 599 c         write(7,1239)nresA(i-1),nresA(i)
 600 c      enddo
 601 c      do i=1,nseqB
 602 c         write(7,1237)2000+nresB(i),seqB(i),nresB(i),xb(i),yb(i),zb(i)
 603 c      enddo
 604 c      write(7,1238)
 605 c      do i=2,nseqB
 606 c         write(7,1239)2000+nresB(i-1),2000+nresB(i)
 607 c      enddo
 608 c 1237 format('ATOM  ',i5,'  CA  ',A3,I6,4X,3F8.3)
 609 c 1238 format('TER')
 610 c 1239 format('CONECT',I5,I5)
 611 c 999  continue
 612 c
 613 c***   record aligned residues by i=[1,nseqA], for sequenceM()------------>
 614 c      do i=1,nseqA
 615 c         iq(i)=0
 616 c      enddo
 617 c      do i=1,n_cut
 618 c         j=iA(i_ali(i))         ![1,nseqA]
 619 c         k=iB(i_ali(i))         ![1,nseqB]
 620 c         dis=sqrt((xt(j)-xb(k))**2+(yt(j)-yb(k))**2+(zt(j)-zb(k))**2)
 621 c         if(dis.lt.d_output)then
 622 c            iq(j)=1
 623 c         endif
 624 c      enddo
 625 c*******************************************************************
 626 c***   output aligned sequences
 627 c      k=0
 628 c      i=1
 629 c      j=1
 630 c 800  continue
 631 c      if(i.gt.nseqA.and.j.gt.nseqB)goto 802
 632 c      if(i.gt.nseqA.and.j.le.nseqB)then
 633 c         k=k+1
 634 c         sequenceA(k)='-'
 635 c         sequenceB(k)=seq1B(j)
 636 c         sequenceM(k)=' '
 637 c         j=j+1
 638 c         goto 800
 639 c      endif
 640 c      if(i.le.nseqA.and.j.gt.nseqB)then
 641 c         k=k+1
 642 c         sequenceA(k)=seq1A(i)
 643 c         sequenceB(k)='-'
 644 c         sequenceM(k)=' '
 645 c         i=i+1
 646 c         goto 800
 647 c      endif
 648 c      if(nresA(i).eq.nresB(j))then
 649 c         k=k+1
 650 c         sequenceA(k)=seq1A(i)
 651 c         sequenceB(k)=seq1B(j)
 652 c         if(iq(i).eq.1)then
 653 c            sequenceM(k)=':'
 654 c         else
 655 c            sequenceM(k)=' '
 656 c         endif
 657 c         i=i+1
 658 c         j=j+1
 659 c         goto 800
 660 c      elseif(nresA(i).lt.nresB(j))then
 661 c         k=k+1
 662 c         sequenceA(k)=seq1A(i)
 663 c         sequenceB(k)='-'
 664 c         sequenceM(k)=' '
 665 c         i=i+1
 666 c         goto 800
 667 c      elseif(nresB(j).lt.nresA(i))then
 668 c         k=k+1
 669 c         sequenceA(k)='-'
 670 c         sequenceB(k)=seq1B(j)
 671 c         sequenceM(k)=' '
 672 c         j=j+1
 673 c         goto 800
 674 c      endif
 675 c 802  continue
 676 c
 677 c      write(*,600)d_output,n_cut,rmsd
 678 c 600  format('Superposition in the TM-score: Length(d<',f3.1,
 679 c     $     ')=',i3,'  RMSD=',f6.2)
 680 c      write(*,603)d_output
 681 c 603  format('(":" denotes the residue pairs of distance < ',f3.1,
 682 c     &     ' Angstrom)')
 683 c      write(*,601)(sequenceA(i),i=1,k)
 684 c      write(*,601)(sequenceM(i),i=1,k)
 685 c      write(*,601)(sequenceB(i),i=1,k)
 686 c      write(*,602)(mod(i,10),i=1,k)
 687 c 601  format(2000A1)
 688 c 602  format(2000I1)
 689 c      write(*,*)
 690 c
 691 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 692 c 9999 END
 693
 694 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 695 c     1, collect those residues with dis<d;
 696 c     2, calculate score_GDT, score_maxsub, score_TM
 697 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 698       subroutine score_fun
 699       PARAMETER(nmax=5000)
 700
 701       common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
 702       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
 703       common/para/d,d0,d0_fix
 704       common/align/n_ali,iA(nmax),iB(nmax)
 705       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
 706       common/scores/score,score_maxsub,score_fix,score10
 707       common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
 708       double precision score,score_max,score_fix,score_fix_max
 709       double precision score_maxsub,score10
 710
 711       d_tmp=d
 712  21   n_cut=0                   !number of residue-pairs dis<d, for iteration
 713       n_GDT05=0                 !for GDT-score, # of dis<0.5
 714       n_GDT1=0                  !for GDT-score, # of dis<1
 715       n_GDT2=0                  !for GDT-score, # of dis<2
 716       n_GDT4=0                  !for GDT-score, # of dis<4
 717       n_GDT8=0                  !for GDT-score, # of dis<8
 718       score_maxsub_sum=0        !Maxsub-score
 719       score_sum=0               !TMscore
 720       score_sum10=0             !TMscore10
 721       do k=1,n_ali
 722          i=iA(k)                ![1,nseqA] reoder number of structureA
 723          j=iB(k)                ![1,nseqB]
 724          dis=sqrt((xt(i)-xb(j))**2+(yt(i)-yb(j))**2+(zt(i)-zb(j))**2)
 725 ***   for iteration:
 726          if(dis.lt.d_tmp)then
 727             n_cut=n_cut+1
 728             i_ali(n_cut)=k      ![1,n_ali], mark the residue-pairs in dis<d
 729          endif
 730 ***   for GDT-score:
 731          if(dis.le.8)then
 732             n_GDT8=n_GDT8+1
 733             if(dis.le.4)then
 734                n_GDT4=n_GDT4+1
 735                if(dis.le.2)then
 736                   n_GDT2=n_GDT2+1
 737                   if(dis.le.1)then
 738                      n_GDT1=n_GDT1+1
 739                      if(dis.le.0.5)then
 740                         n_GDT05=n_GDT05+1
 741                      endif
 742                   endif
 743                endif
 744             endif
 745          endif
 746 ***   for MAXsub-score:
 747          if(dis.lt.3.5)then
 748             score_maxsub_sum=score_maxsub_sum+1/(1+(dis/3.5)**2)
 749          endif
 750 ***   for TM-score:
 751          score_sum=score_sum+1/(1+(dis/d0)**2)
 752 ***   for TM-score10:
 753          if(dis.lt.10)then
 754             score_sum10=score_sum10+1/(1+(dis/d0)**2)
 755          endif
 756       enddo
 757       if(n_cut.lt.3.and.n_ali.gt.3)then
 758          d_tmp=d_tmp+.5
 759          goto 21
 760       endif
 761       score_maxsub=score_maxsub_sum/float(nseqB) !MAXsub-score
 762       score=score_sum/float(nseqB) !TM-score
 763       score10=score_sum10/float(nseqB) !TM-score10
 764
 765       return
 766       end
 767
 768 cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
 769 c  w    - w(m) is weight for atom pair  c m                    (given)
 770 c  x    - x(i,m) are coordinates of atom c m in set x          (given)
 771 c  y    - y(i,m) are coordinates of atom c m in set y          (given)
 772 c  n    - n is number of atom pairs                            (given)
 773 c  mode  - 0:calculate rms     only                            (given,short)
 774 c          1:calculate     u,t only                            (given,medium)
 775 c          2:calculate rms,u,t                                 (given,longer)
 776 c  rms   - sum of w*(ux+t-y)**2 over all atom pairs            (result)
 777 c  u    - u(i,j) is   rotation  matrix for best superposition  (result)
 778 c  t    - t(i)   is translation vector for best superposition  (result)
 779 c  ier  - 0: a unique optimal superposition has been determined(result)
 780 c       -1: superposition is not unique but optimal
 781 c       -2: no result obtained because of negative weights w
 782 c           or all weights equal to zero.
 783 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 784       subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
 785       double precision w(*), x(3,*), y(3,*)
 786       integer n, mode
 787
 788       double precision rms, u(3,3), t(3)
 789       integer ier
 790
 791       integer i, j, k, l, m1, m
 792       integer ip(9), ip2312(4)
 793       double precision r(3,3), xc(3), yc(3), wc
 794       double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
 795       double precision e0, d, spur, det, cof, h, g
 796       double precision cth, sth, sqrth, p, sigma
 797       double precision c1x, c1y, c1z, c2x, c2y, c2z
 798       double precision s1x, s1y, s1z, s2x, s2y, s2z
 799       double precision sxx, sxy, sxz, syx, syy, syz, szx, szy, szz
 800
 801       double precision sqrt3, tol, zero
 802
 803       data sqrt3 / 1.73205080756888d+00 /
 804       data tol / 1.0d-2 /
 805       data zero / 0.0d+00 /
 806       data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
 807       data ip2312 / 2, 3, 1, 2 /
 808
 809       wc  = zero
 810       rms = zero
 811       e0  = zero
 812       s1x = zero
 813       s1y = zero
 814       s1z = zero
 815       s2x = zero
 816       s2y = zero
 817       s2z = zero
 818       sxx = zero
 819       sxy = zero
 820       sxz = zero
 821       syx = zero
 822       syy = zero
 823       syz = zero
 824       szx = zero
 825       szy = zero
 826       szz = zero
 827
 828       do i=1, 3
 829          xc(i) = zero
 830          yc(i) = zero
 831          t(i) = zero
 832          do j=1, 3
 833             r(i,j) = zero
 834             u(i,j) = zero
 835             a(i,j) = zero
 836             if( i .eq. j ) then
 837                u(i,j) = 1.0
 838                a(i,j) = 1.0
 839             end if
 840          end do
 841       end do
 842
 843       ier = -1
 844       if( n .lt. 1 ) return
 845       ier = -2
 846
 847       do m=1, n
 848          c1x=x(1, m)
 849          c1y=x(2, m)
 850          c1z=x(3, m)
 851
 852          c2x=y(1, m)
 853          c2y=y(2, m)
 854          c2z=y(3, m)
 855
 856          s1x = s1x + c1x
 857          s1y = s1y + c1y;
 858          s1z = s1z + c1z;
 859
 860          s2x = s2x + c2x;
 861          s2y = s2y + c2y;
 862          s2z = s2z + c2z;
 863
 864          sxx = sxx + c1x*c2x;
 865          sxy = sxy + c1x*c2y;
 866          sxz = sxz + c1x*c2z;
 867
 868          syx = syx + c1y*c2x;
 869          syy = syy + c1y*c2y;
 870          syz = syz + c1y*c2z;
 871
 872          szx = szx + c1z*c2x;
 873          szy = szy + c1z*c2y;
 874          szz = szz + c1z*c2z;
 875       end do
 876
 877       xc(1) = s1x/n;
 878       xc(2) = s1y/n;
 879       xc(3) = s1z/n;
 880
 881       yc(1) = s2x/n;
 882       yc(2) = s2y/n;
 883       yc(3) = s2z/n;
 884       if(mode.eq.2.or.mode.eq.0) then ! need rmsd
 885          do m=1, n
 886             do i=1, 3
 887                e0 = e0+ (x(i, m)-xc(i))**2 + (y(i, m)-yc(i))**2
 888             end do
 889          end do
 890       endif
 891
 892       r(1, 1) = sxx-s1x*s2x/n;
 893       r(2, 1) = sxy-s1x*s2y/n;
 894       r(3, 1) = sxz-s1x*s2z/n;
 895       r(1, 2) = syx-s1y*s2x/n;
 896       r(2, 2) = syy-s1y*s2y/n;
 897       r(3, 2) = syz-s1y*s2z/n;
 898       r(1, 3) = szx-s1z*s2x/n;
 899       r(2, 3) = szy-s1z*s2y/n;
 900       r(3, 3) = szz-s1z*s2z/n;
 901
 902       det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
 903      &     - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
 904      &     + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
 905
 906       sigma = det
 907
 908       m = 0
 909       do j=1, 3
 910          do i=1, j
 911             m = m+1
 912             rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
 913          end do
 914       end do
 915
 916       spur = (rr(1)+rr(3)+rr(6)) / 3.0
 917       cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
 918      &     - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
 919       det = det*det
 920
 921       do i=1, 3
 922          e(i) = spur
 923       end do
 924       if( spur .le. zero ) goto 40
 925       d = spur*spur
 926       h = d - cof
 927       g = (spur*cof - det)/2.0 - spur*h
 928       if( h .le. zero ) then
 929          if( mode .eq. 0 ) then
 930             goto 50
 931          else
 932             goto 30
 933          end if
 934       end if
 935       sqrth = dsqrt(h)
 936       d = h*h*h - g*g
 937       if( d .lt. zero ) d = zero
 938       d = datan2( dsqrt(d), -g ) / 3.0
 939       cth = sqrth * dcos(d)
 940       sth = sqrth*sqrt3*dsin(d)
 941       e(1) = (spur + cth) + cth
 942       e(2) = (spur - cth) + sth
 943       e(3) = (spur - cth) - sth
 944
 945       if( mode .eq. 0 ) then
 946          goto 50
 947       end if
 948
 949       do l=1, 3, 2
 950          d = e(l)
 951          ss(1) = (d-rr(3)) * (d-rr(6))  - rr(5)*rr(5)
 952          ss(2) = (d-rr(6)) * rr(2)      + rr(4)*rr(5)
 953          ss(3) = (d-rr(1)) * (d-rr(6))  - rr(4)*rr(4)
 954          ss(4) = (d-rr(3)) * rr(4)      + rr(2)*rr(5)
 955          ss(5) = (d-rr(1)) * rr(5)      + rr(2)*rr(4)
 956          ss(6) = (d-rr(1)) * (d-rr(3))  - rr(2)*rr(2)
 957
 958          if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
 959             j=1
 960             if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
 961          else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
 962             j = 2
 963          else
 964             j = 3
 965          end if
 966
 967          d = zero
 968          j = 3 * (j - 1)
 969
 970          do i=1, 3
 971             k = ip(i+j)
 972             a(i,l) = ss(k)
 973             d = d + ss(k)*ss(k)
 974          end do
 975          if( d .gt. zero ) d = 1.0 / dsqrt(d)
 976          do i=1, 3
 977             a(i,l) = a(i,l) * d
 978          end do
 979       end do
 980
 981       d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
 982       if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
 983          m1 = 3
 984          m = 1
 985       else
 986          m1 = 1
 987          m = 3
 988       endif
 989
 990       p = zero
 991       do i=1, 3
 992          a(i,m1) = a(i,m1) - d*a(i,m)
 993          p = p + a(i,m1)**2
 994       end do
 995       if( p .le. tol ) then
 996          p = 1.0
 997          do 21 i=1, 3
 998             if (p .lt. dabs(a(i,m))) goto 21
 999             p = dabs( a(i,m) )
1000             j = i
1001  21      continue
1002          k = ip2312(j)
1003          l = ip2312(j+1)
1004          p = dsqrt( a(k,m)**2 + a(l,m)**2 )
1005          if( p .le. tol ) goto 40
1006          a(j,m1) = zero
1007          a(k,m1) = -a(l,m)/p
1008          a(l,m1) =  a(k,m)/p
1009       else
1010          p = 1.0 / dsqrt(p)
1011          do i=1, 3
1012             a(i,m1) = a(i,m1)*p
1013          end do
1014       end if
1015
1016       a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
1017       a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
1018       a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
1019
1020  30   do l=1, 2
1021          d = zero
1022          do i=1, 3
1023             b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
1024             d = d + b(i,l)**2
1025          end do
1026          if( d .gt. zero ) d = 1.0 / dsqrt(d)
1027          do i=1, 3
1028             b(i,l) = b(i,l)*d
1029          end do
1030       end do
1031       d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
1032       p = zero
1033
1034       do i=1, 3
1035          b(i,2) = b(i,2) - d*b(i,1)
1036          p = p + b(i,2)**2
1037       end do
1038       if( p .le. tol ) then
1039          p = 1.0
1040          do 22 i=1, 3
1041             if(p.lt.dabs(b(i,1)))goto 22
1042             p = dabs( b(i,1) )
1043             j = i
1044  22      continue
1045          k = ip2312(j)
1046          l = ip2312(j+1)
1047          p = dsqrt( b(k,1)**2 + b(l,1)**2 )
1048          if( p .le. tol ) goto 40
1049          b(j,2) = zero
1050          b(k,2) = -b(l,1)/p
1051          b(l,2) =  b(k,1)/p
1052       else
1053          p = 1.0 / dsqrt(p)
1054          do i=1, 3
1055             b(i,2) = b(i,2)*p
1056          end do
1057       end if
1058
1059       b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
1060       b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
1061       b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
1062
1063       do i=1, 3
1064          do j=1, 3
1065             u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
1066          end do
1067       end do
1068
1069  40   do i=1, 3
1070          t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
1071       end do
1072  50   do i=1, 3
1073          if( e(i) .lt. zero ) e(i) = zero
1074          e(i) = dsqrt( e(i) )
1075       end do
1076
1077       ier = 0
1078       if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
1079
1080       d = e(3)
1081       if( sigma .lt. 0.0 ) then
1082          d = - d
1083          if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
1084       end if
1085       d = (d + e(2)) + e(1)
1086
1087       if(mode .eq. 2.or.mode.eq.0) then ! need rmsd
1088          rms = (e0 - d) - d
1089          if( rms .lt. 0.0 ) rms = 0.0
1090       endif
1091
1092       return
1093       end
1094
1095