binaries/src/tcoffee/t_coffee_source/TMalign.f

   1 *************************************************************************\r
   2 *     This program is to identify the best alignment of two protein \r
   3 *     structures to give the best TM-score. By default, TM-score is \r
   4 *     normalized by the second protein. The program can be freely \r
   5 *     copied or modified or redistributed.\r
   6 *     (For comments, please email to: yzhang@ku.edu)\r
   7 *     \r
   8 *     Reference:\r
   9 *     Yang Zhang, Jeffrey Skolnick, Nucl. Acid Res. 2005 33: 2303-9\r
  10 *\r
  11 ************************ updating history *******************************\r
  12 *     2005/06/01: A small bug of two-point superposition was fixed.\r
  13 *     2005/10/19: the program was reformed so that the alignment\r
  14 *                 results are not dependent on the specific compilers.\r
  15 *     2006/06/20: select 'A' if there is altLoc when reading PDB file.\r
  16 *     2007/02/27: rotation matrix from Chain-1 to Chain-2 is added.\r
  17 *     2007/04/18: add options with TM-score normalized by average\r
  18 *                 length, shorter length, or longer length of two \r
  19 *                 structures.\r
  20 *     2007/05/23: add additional output file 'TM.sup_all' for showing\r
  21 *                 all atoms while 'TM.sup' is only for aligned atoms\r
  22 *     2007/09/19: add a new feature alignment to deal with the problem\r
  23 *                 of aligning fractional structures (e.g. protein\r
  24 *                 interfaces).\r
  25 *************************************************************************\r
  26       \r
  27       program compares\r
  28       PARAMETER(nmax=5000)\r
  29       PARAMETER(nmax2=10000)\r
  30 \r
  31       COMMON/BACKBONE/XA(3,nmax,0:1)\r
  32       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
  33       common/alignrst/invmap0(nmax)\r
  34       common/length/nseq1,nseq2\r
  35       common/d0/d0,anseq\r
  36       common/d0min/d0_min\r
  37       common/d00/d00,d002\r
  38 \r
  39       character*100 fnam,pdb(100),outname\r
  40       character*3 aa(-1:20),aanam,ss1(nmax),ss2(nmax)\r
  41       character*100 s,du\r
  42       character*200 outnameall_tmp,outnameall\r
  43       character seq1(0:nmax),seq2(0:nmax)\r
  44       character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2)\r
  45 \r
  46       dimension xx(nmax),yy(nmax),zz(nmax)\r
  47       dimension m1(nmax),m2(nmax)\r
  48       dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
  49       dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
  50       common/init/invmap_i(nmax)\r
  51 \r
  52       common/TM/TM,TMmax\r
  53       common/n1n2/n1(nmax),n2(nmax)\r
  54       common/d8/d8\r
  55       common/initial4/mm1(nmax),mm2(nmax)\r
  56 \r
  57 ccc   RMSD:\r
  58       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
  59       double precision u(3,3),t(3),rms,drms !armsd is real\r
  60       data w /nmax*1.0/\r
  61 ccc   \r
  62 \r
  63       data aa/ 'BCK','GLY','ALA','SER','CYS','VAL','THR','ILE',\r
  64      &     'PRO','MET','ASP','ASN','LEU',\r
  65      &     'LYS','GLU','GLN','ARG',\r
  66      &     'HIS','PHE','TYR','TRP','CYX'/\r
  67       character*1 slc(-1:20)\r
  68       data slc/'X','G','A','S','C','V','T','I',\r
  69      &     'P','M','D','N','L','K','E','Q','R',\r
  70      &     'H','F','Y','W','C'/\r
  71 \r
  72       call getarg(1,fnam)\r
  73       if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then\r
  74          write(*,*)\r
  75          write(*,*)'Brief instruction for running TM-align program:'\r
  76          write(*,*)'(For detail: Zhang & Skolnick, Nucl. Acid Res.',\r
  77      &     '2005 33, 2303)'\r
  78          write(*,*)\r
  79          write(*,*)'1. Align ''structure.pdb'' to ''target.pdb'''\r
  80          write(*,*)'  (By default, TM-score is normalized by the ',\r
  81      &        'length of ''target.pdb'')'\r
  82          write(*,*)'  >TMalign structure.pdb target.pdb'\r
  83          write(*,*)\r
  84          write(*,*)'2. Run TM-align and output the superposition ',\r
  85      &        'to ''TM.sup'' and ''TM.sup_all'':'\r
  86          write(*,*)'  >TMalign structure.pdb target.pdb -o TM.sup'\r
  87          write(*,*)'   To view the superimposed structures of the',\r
  88      &        ' aligned regions by rasmol:'\r
  89          write(*,*)'  >rasmol -script TM.sup)'\r
  90          write(*,*)'   To view the superimposed structures of all',\r
  91      &        ' regions by rasmol:'\r
  92          write(*,*)'  >rasmol -script TM.sup_all)'\r
  93          write(*,*)\r
  94          write(*,*)'3. If you want TM-score normalized by ',\r
  95      &        'an assigned length, e.g. 100 aa:'\r
  96          write(*,*)'  >TMalign structure.pdb target.pdb -L 100'\r
  97          write(*,*)'   If you want TM-score normalized by the ',\r
  98      &        'average length of two structures:'\r
  99          write(*,*)'  >TMalign structure.pdb target.pdb -a'\r
 100          write(*,*)'   If you want TM-score normalized by the ',\r
 101      &        'shorter length of two structures:'\r
 102          write(*,*)'  >TMalign structure.pdb target.pdb -b'\r
 103          write(*,*)'   If you want TM-score normalized by the ',\r
 104      &        'longer length of two structures:'\r
 105          write(*,*)'  >TMalign structure.pdb target.pdb -c'\r
 106          write(*,*)\r
 107 c         write(*,*)'5. If you want to set a minimum cutoff for d0, ',\r
 108 c     &        'e.g. d0>3.0'\r
 109 c         write(*,*)'   (By default d0>0.5, this option need ',\r
 110 c     &        'be considered only when L<35 aa)'\r
 111 c         write(*,*)'  >TMalign structure.pdb target.pdb -dmin 3.0'\r
 112 c         write(*,*)\r
 113          write(*,*)'(All above options does not change the ',\r
 114      &         'final structure alignment result)'\r
 115          write(*,*)\r
 116          goto 9999\r
 117       endif\r
 118 \r
 119 ******* options ----------->\r
 120       m_out=-1                  !decided output\r
 121       m_fix=-1                  !fixed length-scale only for output\r
 122       m_ave=-1                  !using average length\r
 123       m_d0_min=-1               !diminum d0 for search\r
 124       m_d0=-1                   !given d0 for both search and output\r
 125       narg=iargc()\r
 126       i=0\r
 127       j=0\r
 128  115  continue\r
 129       i=i+1\r
 130       call getarg(i,fnam)\r
 131       if(fnam.eq.'-o')then\r
 132          m_out=1\r
 133          i=i+1\r
 134          call getarg(i,outname)\r
 135       elseif(fnam.eq.'-L')then  !change both L_all and d0\r
 136          m_fix=1\r
 137          i=i+1\r
 138          call getarg(i,fnam)\r
 139          read(fnam,*)L_fix\r
 140       elseif(fnam.eq.'-dmin')then\r
 141          m_d0_min=1\r
 142          i=i+1\r
 143          call getarg(i,fnam)\r
 144          read(fnam,*)d0_min_input\r
 145       elseif(fnam.eq.'-d0')then\r
 146          m_d0=1\r
 147          i=i+1\r
 148          call getarg(i,fnam)\r
 149          read(fnam,*)d0_fix\r
 150       elseif(fnam.eq.'-a')then ! this will change the superposed output but not the alignment\r
 151          m_ave=1\r
 152          i=i+1\r
 153       elseif(fnam.eq.'-b')then\r
 154         m_ave=2\r
 155          i=i+1\r
 156       elseif(fnam.eq.'-c')then\r
 157         m_ave=3\r
 158          i=i+1\r
 159       else\r
 160          j=j+1\r
 161          pdb(j)=fnam\r
 162       endif\r
 163       if(i.lt.narg)goto 115\r
 164       \r
 165 ccccccccc read data from first CA file:\r
 166       open(unit=10,file=pdb(1),status='old')\r
 167       i=0\r
 168       do while (.true.)\r
 169          read(10,9001,end=1010) s\r
 170          if(i.gt.0.and.s(1:3).eq.'TER')goto 1010\r
 171          if(s(1:3).eq.'ATO')then\r
 172             if(s(13:16).eq.'CA  '.or.s(13:16).eq.' CA '.or\r
 173      &           .s(13:16).eq.'  CA')then\r
 174                if(s(17:17).eq.' '.or.s(17:17).eq.'A')then\r
 175                   i=i+1\r
 176                   read(s,9000)du,aanam,du,mm1(i),du,\r
 177      $                 xa(1,i,0),xa(2,i,0),xa(3,i,0)\r
 178                   do j=-1,20\r
 179                      if(aanam.eq.aa(j))seq1(i)=slc(j)\r
 180                   enddo\r
 181                   ss1(i)=aanam\r
 182                   if(i.ge.nmax)goto 1010\r
 183                endif\r
 184             endif\r
 185          endif\r
 186       enddo\r
 187  1010 continue\r
 188  9000 format(A17,A3,A2,i4,A4,3F8.3)\r
 189  9001 format(A100)\r
 190       close(10)\r
 191       nseq1=i\r
 192 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
 193 \r
 194 ccccccccc read data from the second CA file:\r
 195       open(unit=10,file=pdb(2),status='old')\r
 196       i=0\r
 197       do while (.true.)\r
 198          read(10,9001,end=1011) s\r
 199          if(i.gt.0.and.s(1:3).eq.'TER')goto 1011\r
 200          if(s(1:3).eq.'ATO')then\r
 201             if(s(13:16).eq.'CA  '.or.s(13:16).eq.' CA '.or.\r
 202      &           s(13:16).eq.'  CA')then\r
 203                if(s(17:17).eq.' '.or.s(17:17).eq.'A')then\r
 204                   i=i+1\r
 205                   read(s,9000)du,aanam,du,mm2(i),du,\r
 206      $                 xa(1,i,1),xa(2,i,1),xa(3,i,1)\r
 207                   do j=-1,20\r
 208                      if(aanam.eq.aa(j))seq2(i)=slc(j)\r
 209                   enddo\r
 210                   ss2(i)=aanam\r
 211                   if(i.ge.nmax)goto 1011\r
 212                endif\r
 213             endif\r
 214          endif\r
 215       enddo\r
 216  1011 continue\r
 217       close(10)\r
 218       nseq2=i\r
 219 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
 220 \r
 221 *!!!  Scale of TM-score in search is based on the smaller protein --------->\r
 222       d0_min=0.5\r
 223       if(m_d0_min.eq.1)then\r
 224          d0_min=d0_min_input    !for search\r
 225       endif\r
 226       anseq_min=min(nseq1,nseq2)\r
 227       anseq=anseq_min           !length for defining TMscore in search\r
 228       d8=1.5*anseq_min**0.3+3.5 !remove pairs with dis>d8 during search & final\r
 229       if(anseq.gt.15)then\r
 230          d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score\r
 231       else\r
 232          d0=d0_min\r
 233       endif\r
 234       if(d0.lt.d0_min)d0=d0_min\r
 235       if(m_d0.eq.1)d0=d0_fix\r
 236       d00=d0                    !for quickly calculate TM-score in searching\r
 237       if(d00.gt.8)d00=8\r
 238       if(d00.lt.4.5)d00=4.5\r
 239       d002=d00**2\r
 240       nseq=max(nseq1,nseq2)\r
 241       do i=1,nseq\r
 242          n1(i)=i\r
 243          n2(i)=i\r
 244       enddo\r
 245       \r
 246 ***** do alignment **************************\r
 247       CALL super_align          !to find invmap(j)\r
 248       \r
 249 ************************************************************\r
 250 ***   resuperpose to find residues of dis<d8 ------------------------>\r
 251       n_al=0\r
 252       do j=1,nseq2\r
 253          if(invmap0(j).gt.0)then\r
 254             i=invmap0(j)\r
 255             n_al=n_al+1\r
 256             xtm1(n_al)=xa(1,i,0)\r
 257             ytm1(n_al)=xa(2,i,0)\r
 258             ztm1(n_al)=xa(3,i,0)\r
 259             xtm2(n_al)=xa(1,j,1)\r
 260             ytm2(n_al)=xa(2,j,1)\r
 261             ztm2(n_al)=xa(3,j,1)\r
 262             m1(n_al)=i          !for recording residue order\r
 263             m2(n_al)=j\r
 264          endif\r
 265       enddo\r
 266       d0_input=d0\r
 267       call TMscore8(d0_input,n_al,xtm1,ytm1,ztm1,n1,n_al,\r
 268      &     xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !TM-score with dis<d8 only\r
 269 \r
 270 *!!!  Output TM-score is based on the second protein------------------>\r
 271       d0_min=0.5                !for output\r
 272       anseq=nseq2               !length for defining final TMscore\r
 273       if(m_ave.eq.1)anseq=(nseq1+nseq2)/2.0 !<L>\r
 274       if(m_ave.eq.2)anseq=min(nseq1,nseq2)\r
 275       if(m_ave.eq.3)anseq=max(nseq1,nseq2)\r
 276       if(anseq.lt.anseq_min)anseq=anseq_min\r
 277       if(m_fix.eq.1)anseq=L_fix !input length\r
 278       if(anseq.gt.15)then\r
 279          d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score\r
 280       else\r
 281          d0=d0_min\r
 282       endif\r
 283       if(d0.lt.d0_min)d0=d0_min\r
 284       if(m_d0.eq.1)d0=d0_fix\r
 285       \r
 286 ***   remove dis>d8 in normal TM-score calculation for final report----->\r
 287       j=0\r
 288       n_eq=0\r
 289       do i=1,n_al\r
 290          dis2=sqrt((xtm1(i)-xtm2(i))**2+(ytm1(i)-ytm2(i))**2+\r
 291      &        (ztm1(i)-ztm2(i))**2)\r
 292          if(dis2.le.d8)then\r
 293             j=j+1\r
 294             xtm1(j)=xtm1(i)\r
 295             ytm1(j)=ytm1(i)\r
 296             ztm1(j)=ztm1(i)\r
 297             xtm2(j)=xtm2(i)\r
 298             ytm2(j)=ytm2(i)\r
 299             ztm2(j)=ztm2(i)\r
 300             m1(j)=m1(i)\r
 301             m2(j)=m2(i)\r
 302             if(ss1(m1(i)).eq.ss2(m2(i)))then\r
 303                n_eq=n_eq+1\r
 304             endif\r
 305          endif\r
 306       enddo\r
 307       seq_id=float(n_eq)/(n_al+0.00000001)\r
 308       n8_al=j\r
 309       d0_input=d0\r
 310       call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n1,n8_al,\r
 311      &     xtm2,ytm2,ztm2,n2,TM8,Rcomm,Lcomm) !normal TMscore\r
 312       rmsd8_al=Rcomm\r
 313       TM8=TM8*n8_al/anseq       !TM-score after cutoff\r
 314       \r
 315 ********* for output summary ******************************\r
 316       write(*,*)\r
 317       write(*,*)'*****************************************************',\r
 318      &     '*********************'\r
 319       write(*,*)'*                               TM-align             ',\r
 320      &     '                    *'\r
 321       write(*,*)'* A protein structural alignment algorithm based on T',\r
 322      &     'M-score             *'\r
 323       write(*,*)'* Reference: Y. Zhang and J. Skolnick, Nucl. Acids Re',\r
 324      &     's. 2005 33, 2302-9  *'\r
 325       write(*,*)'* Comments on the program, please email to: yzhang@ku',\r
 326      &     '.edu                *'\r
 327       write(*,*)'*****************************************************',\r
 328      &     '*********************'\r
 329       write(*,*)\r
 330       write(*,101)pdb(1),nseq1\r
 331  101  format('Chain 1:',A10,'  Size=',I4)\r
 332       write(*,102)pdb(2),nseq2,int(anseq)\r
 333  102  format('Chain 2:',A10,'  Size=',I4,\r
 334      &     ' (TM-score is normalized by ',I4,')')\r
 335       write(*,*)\r
 336       write(*,103)n8_al,rmsd8_al,TM8,seq_id\r
 337  103  format('Aligned length=',I4,', RMSD=',f6.2,\r
 338      &     ', TM-score=',f7.5,', ID=',f5.3)\r
 339       write(*,*)\r
 340 \r
 341 ********* extract rotation matrix ------------>\r
 342       L=0\r
 343       do i=1,n8_al\r
 344          k=m1(i)\r
 345          L=L+1\r
 346          r_1(1,L)=xa(1,k,0)\r
 347          r_1(2,L)=xa(2,k,0)\r
 348          r_1(3,L)=xa(3,k,0)\r
 349          r_2(1,L)=xtm1(i)\r
 350          r_2(2,L)=ytm1(i)\r
 351          r_2(3,L)=ztm1(i)\r
 352        enddo\r
 353        if(L.gt.3)then\r
 354          call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2\r
 355          armsd=dsqrt(rms/L)\r
 356          write(*,*)'-------- rotation matrix to rotate Chain-1 to ',\r
 357      &        'Chain-2 ------'\r
 358          write(*,*)'i          t(i)         u(i,1)         u(i,2) ',\r
 359      &        '        u(i,3)'\r
 360          do i=1,3\r
 361             write(*,204)i,t(i),u(i,1),u(i,2),u(i,3)\r
 362          enddo\r
 363 c         do i=1,nseq1\r
 364 c            ax=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)\r
 365 c            ay=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)\r
 366 c            az=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)\r
 367 c         enddo\r
 368          write(*,*)\r
 369       endif\r
 370  204  format(I2,f18.10,f15.10,f15.10,f15.10)\r
 371       \r
 372 ********* for output superposition ******************************\r
 373       if(m_out.eq.1)then\r
 374  1237    format('ATOM  ',i5,'  CA  ',A3,I6,4X,3F8.3)\r
 375  1238    format('TER')\r
 376  1239    format('CONECT',I5,I5)\r
 377  900     format(A)\r
 378  901     format('select atomno=',I4)\r
 379  104     format('REMARK Chain 1:',A10,'  Size=',I4)\r
 380  105     format('REMARK Chain 2:',A10,'  Size=',I4,\r
 381      &        ' (TM-score is normalized by ',I4,')')\r
 382  106     format('REMARK Aligned length=',I4,', RMSD=',f6.2,\r
 383      &        ', TM-score=',f7.5,', ID=',f5.3)\r
 384          OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln\r
 385 ***   script:\r
 386          write(7,900)'load inline'\r
 387          write(7,900)'select atomno<2000'\r
 388          write(7,900)'wireframe .45'\r
 389          write(7,900)'select none'\r
 390          write(7,900)'select atomno>2000'\r
 391          write(7,900)'wireframe .20'\r
 392          write(7,900)'color white'\r
 393          do i=1,n8_al\r
 394             dis2=sqrt((xtm1(i)-xtm2(i))**2+\r
 395      &           (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)\r
 396             if(dis2.le.5)then\r
 397                write(7,901)m1(i)\r
 398                write(7,900)'color red'\r
 399                write(7,901)2000+m2(i)\r
 400                write(7,900)'color red'\r
 401             endif\r
 402          enddo\r
 403          write(7,900)'select all'\r
 404          write(7,900)'exit'\r
 405          write(7,104)pdb(1),nseq1\r
 406          write(7,105)pdb(2),nseq2,int(anseq)\r
 407          write(7,106)n8_al,rmsd8_al,TM8,seq_id\r
 408 ***   chain1:\r
 409          do i=1,n8_al\r
 410             write(7,1237)m1(i),ss1(m1(i)),mm1(m1(i)),\r
 411      &           xtm1(i),ytm1(i),ztm1(i)\r
 412          enddo\r
 413          write(7,1238)          !TER\r
 414          do i=2,n8_al\r
 415             write(7,1239)m1(i-1),m1(i) !connect atoms\r
 416          enddo\r
 417 ***   chain2:\r
 418          do i=1,n8_al\r
 419             write(7,1237)2000+m2(i),ss2(m2(i)),mm2(m2(i)),\r
 420      $           xtm2(i),ytm2(i),ztm2(i)\r
 421          enddo\r
 422          write(7,1238)\r
 423          do i=2,n8_al\r
 424             write(7,1239)2000+m2(i-1),2000+m2(i)\r
 425          enddo\r
 426          close(7)\r
 427 ccc   \r
 428          k=0\r
 429          outnameall_tmp=outname//'_all'\r
 430          outnameall=''\r
 431          do i=1,200\r
 432             if(outnameall_tmp(i:i).ne.' ')then\r
 433                k=k+1\r
 434                outnameall(k:k)=outnameall_tmp(i:i)\r
 435             endif\r
 436          enddo\r
 437          OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln\r
 438 ***   script:\r
 439          write(8,900)'load inline'\r
 440          write(8,900)'select atomno<2000'\r
 441          write(8,900)'wireframe .45'\r
 442          write(8,900)'select none'\r
 443          write(8,900)'select atomno>2000'\r
 444          write(8,900)'wireframe .20'\r
 445          write(8,900)'color white'\r
 446          do i=1,n8_al\r
 447             dis2=sqrt((xtm1(i)-xtm2(i))**2+\r
 448      &           (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)\r
 449             if(dis2.le.5)then\r
 450                write(8,901)m1(i)\r
 451                write(8,900)'color red'\r
 452                write(8,901)2000+m2(i)\r
 453                write(8,900)'color red'\r
 454             endif\r
 455          enddo\r
 456          write(8,900)'select all'\r
 457          write(8,900)'exit'\r
 458          write(8,104)pdb(1),nseq1\r
 459          write(8,105)pdb(2),nseq2,int(anseq)\r
 460          write(8,106)n8_al,rmsd8_al,TM8,seq_id\r
 461 ***   chain1:\r
 462          do i=1,nseq1\r
 463             ax=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)\r
 464             ay=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)\r
 465             az=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)\r
 466             write(8,1237)i,ss1(i),mm1(i),ax,ay,az\r
 467          enddo\r
 468          write(8,1238)          !TER\r
 469          do i=2,nseq1\r
 470             write(8,1239)i-1,i\r
 471          enddo\r
 472 ***   chain2:\r
 473          do i=1,nseq2\r
 474             write(8,1237)2000+i,ss2(i),mm2(i),\r
 475      $           xa(1,i,1),xa(2,i,1),xa(3,i,1)\r
 476          enddo\r
 477          write(8,1238)\r
 478          do i=2,nseq2\r
 479             write(8,1239)2000+i-1,2000+i\r
 480          enddo\r
 481          close(8)\r
 482       endif\r
 483 *^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
 484 \r
 485 ************  output aligned sequences **************************\r
 486       ii=0\r
 487       i1_old=1\r
 488       i2_old=1\r
 489       do i=1,n8_al\r
 490          do j=i1_old,m1(i)-1\r
 491             ii=ii+1\r
 492             aseq1(ii)=seq1(j)\r
 493             aseq2(ii)='-'\r
 494             aseq3(ii)=' '\r
 495          enddo\r
 496          do j=i2_old,m2(i)-1\r
 497             ii=ii+1\r
 498             aseq1(ii)='-'\r
 499             aseq2(ii)=seq2(j)\r
 500             aseq3(ii)=' '\r
 501          enddo\r
 502          ii=ii+1\r
 503          aseq1(ii)=seq1(m1(i))\r
 504          aseq2(ii)=seq2(m2(i))\r
 505          dis2=sqrt((xtm1(i)-xtm2(i))**2+\r
 506      &     (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)\r
 507          if(dis2.le.5)then\r
 508            aseq3(ii)=':'\r
 509          else\r
 510            aseq3(ii)='.'\r
 511          endif\r
 512          i1_old=m1(i)+1\r
 513          i2_old=m2(i)+1\r
 514       enddo\r
 515       do i=i1_old,nseq1\r
 516          ii=ii+1\r
 517          aseq1(ii)=seq1(i)\r
 518          aseq2(ii)='-'\r
 519          aseq3(ii)=' '\r
 520       enddo\r
 521       do i=i2_old,nseq2\r
 522          ii=ii+1\r
 523          aseq1(ii)='-'\r
 524          aseq2(ii)=seq2(i)\r
 525          aseq3(ii)=' '\r
 526       enddo\r
 527       write(*,50)\r
 528  50   format('(":" denotes the residue pairs of distance < 5.0 ',\r
 529      &     'Angstrom)')\r
 530       write(*,10)(aseq1(i),i=1,ii)\r
 531       write(*,10)(aseq3(i),i=1,ii)\r
 532       write(*,10)(aseq2(i),i=1,ii)\r
 533  10   format(10000A1)\r
 534       write(*,*)\r
 535 \r
 536 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
 537  9999 END\r
 538 \r
 539 ***********************************************************************\r
 540 ***********************************************************************\r
 541 *     Structure superposition\r
 542 ***********************************************************************\r
 543 ***********************************************************************\r
 544 ***********************************************************************\r
 545       SUBROUTINE super_align\r
 546       PARAMETER(nmax=5000)\r
 547       COMMON/BACKBONE/XA(3,nmax,0:1)\r
 548       common/length/nseq1,nseq2\r
 549       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
 550       common/alignrst/invmap0(nmax)\r
 551       common/zscore/zrms,n_al,rmsd_al\r
 552       common/TM/TM,TMmax\r
 553       common/init/invmap_i(nmax)\r
 554       dimension gapp(100)\r
 555 \r
 556       TMmax=0\r
 557       n_gapp=2\r
 558       gapp(1)=-0.6\r
 559       gapp(2)=0\r
 560 \r
 561 c      n_gapp=11\r
 562 c      do i=1,n_gapp\r
 563 c         gapp(i)=-(n_gapp-i)\r
 564 c      enddo\r
 565 \r
 566 *11111111111111111111111111111111111111111111111111111111\r
 567 *     get initial alignment from gapless threading\r
 568 **********************************************************\r
 569       call get_initial          !gapless threading\r
 570       do i=1,nseq2\r
 571          invmap(i)=invmap_i(i)  !with highest zcore\r
 572       enddo\r
 573       call get_score            !TM, matrix score(i,j)\r
 574       if(TM.gt.TMmax)then\r
 575          TMmax=TM\r
 576          do j=1,nseq2\r
 577             invmap0(j)=invmap(j)\r
 578          enddo\r
 579       endif\r
 580 \r
 581 *****************************************************************\r
 582 *       initerative alignment, for different gap_open:\r
 583 *****************************************************************\r
 584       DO 111 i_gapp=1,n_gapp    !different gap panalties\r
 585          GAP_OPEN=gapp(i_gapp)  !gap panalty\r
 586          do 222 id=1,30         !maximum interation is 200\r
 587             call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
 588 *     Input: score(i,j), and gap_open\r
 589 *     Output: invmap(j)\r
 590 \r
 591             call get_score      !calculate TM-score, score(i,j)\r
 592 c     record the best alignment in whole search ---------->\r
 593             if(TM.gt.TMmax)then\r
 594                TMmax=TM\r
 595                do j=1,nseq2\r
 596                   invmap0(j)=invmap(j)\r
 597                enddo\r
 598             endif\r
 599             if(id.gt.1)then\r
 600                diff=abs(TM-TM_old)\r
 601                if(diff.lt.0.000001)goto 33\r
 602             endif\r
 603             TM_old=TM\r
 604  222     continue\r
 605  33      continue\r
 606  111  continue\r
 607 \r
 608 *222222222222222222222222222222222222222222222222222222222\r
 609 *     get initial alignment from secondary structure alignment\r
 610 **********************************************************\r
 611       call get_initial2         !DP for secondary structure\r
 612       do i=1,nseq2\r
 613          invmap(i)=invmap_i(i)  !with highest zcore\r
 614       enddo\r
 615       call get_score            !TM, score(i,j)\r
 616       if(TM.gt.TMmax)then\r
 617          TMmax=TM\r
 618          do j=1,nseq2\r
 619             invmap0(j)=invmap(j)\r
 620          enddo\r
 621       endif\r
 622 \r
 623 *****************************************************************\r
 624 *     initerative alignment, for different gap_open:\r
 625 *****************************************************************\r
 626       DO 1111 i_gapp=1,n_gapp   !different gap panalties\r
 627          GAP_OPEN=gapp(i_gapp)  !gap panalty\r
 628          do 2222 id=1,30        !maximum interation is 200\r
 629             call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
 630 *     Input: score(i,j), and gap_open\r
 631 *     Output: invmap(j)\r
 632 \r
 633             call get_score      !calculate TM-score, score(i,j)\r
 634 c     write(*,21)gap_open,rmsd_al,n_al,TM\r
 635 c     record the best alignment in whole search ---------->\r
 636             if(TM.gt.TMmax)then\r
 637                TMmax=TM\r
 638                do j=1,nseq2\r
 639                   invmap0(j)=invmap(j)\r
 640                enddo\r
 641             endif\r
 642             if(id.gt.1)then\r
 643                diff=abs(TM-TM_old)\r
 644                if(diff.lt.0.000001)goto 333\r
 645             endif\r
 646             TM_old=TM\r
 647  2222    continue\r
 648  333     continue\r
 649  1111 continue\r
 650       \r
 651 *333333333333333333333333333333333333333333333333333333333333\r
 652 *     get initial alignment from invmap0+SS\r
 653 *************************************************************\r
 654       call get_initial3         !invmap0+SS\r
 655       do i=1,nseq2\r
 656          invmap(i)=invmap_i(i)  !with highest zcore\r
 657       enddo\r
 658       call get_score            !TM, score(i,j)\r
 659       if(TM.gt.TMmax)then\r
 660          TMmax=TM\r
 661          do j=1,nseq2\r
 662             invmap0(j)=invmap(j)\r
 663          enddo\r
 664       endif\r
 665 \r
 666 *****************************************************************\r
 667 *     initerative alignment, for different gap_open:\r
 668 *****************************************************************\r
 669       DO 1110 i_gapp=1,n_gapp   !different gap panalties\r
 670          GAP_OPEN=gapp(i_gapp)  !gap panalty\r
 671          do 2220 id=1,30        !maximum interation is 200\r
 672             call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
 673 *     Input: score(i,j), and gap_open\r
 674 *     Output: invmap(j)\r
 675 \r
 676             call get_score      !calculate TM-score, score(i,j)\r
 677 c     write(*,21)gap_open,rmsd_al,n_al,TM\r
 678 c     record the best alignment in whole search ---------->\r
 679             if(TM.gt.TMmax)then\r
 680                TMmax=TM\r
 681                do j=1,nseq2\r
 682                   invmap0(j)=invmap(j)\r
 683                enddo\r
 684             endif\r
 685             if(id.gt.1)then\r
 686                diff=abs(TM-TM_old)\r
 687                if(diff.lt.0.000001)goto 330\r
 688             endif\r
 689             TM_old=TM\r
 690  2220    continue\r
 691  330     continue\r
 692  1110 continue\r
 693 \r
 694 *444444444444444444444444444444444444444444444444444444444\r
 695 *     get initial alignment of pieces from gapless threading\r
 696 **********************************************************\r
 697       call get_initial4         !gapless threading\r
 698       do i=1,nseq2\r
 699          invmap(i)=invmap_i(i)  !with highest zcore\r
 700       enddo\r
 701       call get_score            !TM, matrix score(i,j)\r
 702       if(TM.gt.TMmax)then\r
 703          TMmax=TM\r
 704          do j=1,nseq2\r
 705             invmap0(j)=invmap(j)\r
 706          enddo\r
 707       endif\r
 708 \r
 709 *****************************************************************\r
 710 *     initerative alignment, for different gap_open:\r
 711 *****************************************************************\r
 712       DO 44 i_gapp=2,n_gapp     !different gap panalties\r
 713          GAP_OPEN=gapp(i_gapp)  !gap panalty\r
 714          do 444 id=1,2          !maximum interation is 200\r
 715             call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
 716 *     Input: score(i,j), and gap_open\r
 717 *     Output: invmap(j)\r
 718             \r
 719             call get_score      !calculate TM-score, score(i,j)\r
 720 c     record the best alignment in whole search ---------->\r
 721             if(TM.gt.TMmax)then\r
 722                TMmax=TM\r
 723                do j=1,nseq2\r
 724                   invmap0(j)=invmap(j)\r
 725                enddo\r
 726             endif\r
 727  444     continue\r
 728  44   continue\r
 729 \r
 730 c^^^^^^^^^^^^^^^ best alignment invmap0(j) found ^^^^^^^^^^^^^^^^^^\r
 731       RETURN\r
 732       END\r
 733 \r
 734 **************************************************************\r
 735 *     get initial alignment invmap0(i) from gapless threading\r
 736 **************************************************************\r
 737       subroutine get_initial\r
 738       PARAMETER(nmax=5000)\r
 739       COMMON/BACKBONE/XA(3,nmax,0:1)\r
 740       common/length/nseq1,nseq2\r
 741       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
 742       common/alignrst/invmap0(nmax)\r
 743       common/zscore/zrms,n_al,rmsd_al\r
 744       common/TM/TM,TMmax\r
 745       common/init/invmap_i(nmax)\r
 746 \r
 747       aL=min(nseq1,nseq2)\r
 748       idel=aL/2.5               !minimum size of considered fragment\r
 749       if(idel.le.5)idel=5\r
 750       n1=-nseq2+idel\r
 751       n2=nseq1-idel\r
 752       GL_max=0\r
 753       do ishift=n1,n2\r
 754          L=0\r
 755          do j=1,nseq2\r
 756             i=j+ishift\r
 757             if(i.ge.1.and.i.le.nseq1)then\r
 758                L=L+1\r
 759                invmap(j)=i\r
 760             else\r
 761                invmap(j)=-1\r
 762             endif\r
 763          enddo\r
 764          if(L.ge.idel)then\r
 765             call get_GL(GL)\r
 766             if(GL.gt.GL_max)then\r
 767                GL_max=GL\r
 768                do i=1,nseq2\r
 769                   invmap_i(i)=invmap(i)\r
 770                enddo\r
 771             endif\r
 772          endif\r
 773       enddo\r
 774 \r
 775       return\r
 776       end\r
 777 \r
 778 **************************************************************\r
 779 *     get initial alignment invmap0(i) from secondary structure\r
 780 **************************************************************\r
 781       subroutine get_initial2\r
 782       PARAMETER(nmax=5000)\r
 783       COMMON/BACKBONE/XA(3,nmax,0:1)\r
 784       common/length/nseq1,nseq2\r
 785       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
 786       common/alignrst/invmap0(nmax)\r
 787       common/zscore/zrms,n_al,rmsd_al\r
 788       common/TM/TM,TMmax\r
 789       common/sec/isec(nmax),jsec(nmax)\r
 790       common/init/invmap_i(nmax)\r
 791 \r
 792 ********** assign secondary structures ***************\r
 793 c     1->coil, 2->helix, 3->turn, 4->strand\r
 794       do i=1,nseq1\r
 795          isec(i)=1\r
 796          j1=i-2\r
 797          j2=i-1\r
 798          j3=i\r
 799          j4=i+1\r
 800          j5=i+2\r
 801          if(j1.ge.1.and.j5.le.nseq1)then\r
 802             dis13=diszy(0,j1,j3)\r
 803             dis14=diszy(0,j1,j4)\r
 804             dis15=diszy(0,j1,j5)\r
 805             dis24=diszy(0,j2,j4)\r
 806             dis25=diszy(0,j2,j5)\r
 807             dis35=diszy(0,j3,j5)\r
 808             isec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35)\r
 809          endif\r
 810       enddo\r
 811       do i=1,nseq2\r
 812          jsec(i)=1\r
 813          j1=i-2\r
 814          j2=i-1\r
 815          j3=i\r
 816          j4=i+1\r
 817          j5=i+2\r
 818          if(j1.ge.1.and.j5.le.nseq2)then\r
 819             dis13=diszy(1,j1,j3)\r
 820             dis14=diszy(1,j1,j4)\r
 821             dis15=diszy(1,j1,j5)\r
 822             dis24=diszy(1,j2,j4)\r
 823             dis25=diszy(1,j2,j5)\r
 824             dis35=diszy(1,j3,j5)\r
 825             jsec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35)\r
 826          endif\r
 827       enddo\r
 828       call smooth               !smooth the assignment\r
 829 \r
 830 ********** score matrix **************************\r
 831       do i=1,nseq1\r
 832          do j=1,nseq2\r
 833             if(isec(i).eq.jsec(j))then\r
 834                score(i,j)=1\r
 835             else\r
 836                score(i,j)=0\r
 837             endif\r
 838          enddo\r
 839       enddo\r
 840 \r
 841 ********** find initial alignment: invmap(j) ************\r
 842       gap_open=-1.0             !should be -1\r
 843       call DP(NSEQ1,NSEQ2)      !produce alignment invmap(j)\r
 844       do i=1,nseq2\r
 845          invmap_i(i)=invmap(i)\r
 846       enddo\r
 847 \r
 848 *^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^\r
 849       return\r
 850       end\r
 851 \r
 852 **************************************************************\r
 853 *     get initial alignment invmap0(i) from secondary structure \r
 854 *     and previous alignments\r
 855 **************************************************************\r
 856       subroutine get_initial3\r
 857       PARAMETER(nmax=5000)\r
 858       COMMON/BACKBONE/XA(3,nmax,0:1)\r
 859       common/length/nseq1,nseq2\r
 860       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
 861       common/alignrst/invmap0(nmax)\r
 862       common/zscore/zrms,n_al,rmsd_al\r
 863       common/TM/TM,TMmax\r
 864       common/sec/isec(nmax),jsec(nmax)\r
 865       common/init/invmap_i(nmax)\r
 866 \r
 867 ********** score matrix **************************\r
 868       do i=1,nseq2\r
 869          invmap(i)=invmap0(i)\r
 870       enddo\r
 871       call get_score1           !get score(i,j) using RMSD martix\r
 872       do i=1,nseq1\r
 873          do j=1,nseq2\r
 874             if(isec(i).eq.jsec(j))then\r
 875                score(i,j)=0.5+score(i,j)\r
 876             else\r
 877                score(i,j)=score(i,j)\r
 878             endif\r
 879          enddo\r
 880       enddo\r
 881 \r
 882 ********** find initial alignment: invmap(j) ************\r
 883       gap_open=-1.0             !should be -1\r
 884       call DP(NSEQ1,NSEQ2)      !produce alignment invmap(j)\r
 885       do i=1,nseq2\r
 886          invmap_i(i)=invmap(i)\r
 887       enddo\r
 888 \r
 889 *^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^\r
 890       return\r
 891       end\r
 892 \r
 893 **************************************************************\r
 894 *     get initial alignment invmap0(i) from fragment gapless threading\r
 895 **************************************************************\r
 896       subroutine get_initial4\r
 897       PARAMETER(nmax=5000)\r
 898       COMMON/BACKBONE/XA(3,nmax,0:1)\r
 899       common/length/nseq1,nseq2\r
 900       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
 901       common/alignrst/invmap0(nmax)\r
 902       common/zscore/zrms,n_al,rmsd_al\r
 903       common/TM/TM,TMmax\r
 904       common/init/invmap_i(nmax)\r
 905       common/initial4/mm1(nmax),mm2(nmax)\r
 906       logical contin\r
 907 \r
 908       dimension ifr2(2,nmax,nmax),Lfr2(2,nmax),Lfr_max2(2),i_fr2(2)\r
 909       dimension ifr(nmax)\r
 910       dimension mm(2,nmax)\r
 911       \r
 912       fra_min=4                 !>=4,minimum fragment for search\r
 913       fra_min1=fra_min-1        !cutoff for shift, save time\r
 914       dcu0=3.85\r
 915       dcu_min=3.65\r
 916 \r
 917 ccc   Find the smallest continuous fragments -------->\r
 918       do i=1,nseq1\r
 919          mm(1,i)=mm1(i)\r
 920       enddo\r
 921       do i=1,nseq2\r
 922          mm(2,i)=mm2(i)\r
 923       enddo\r
 924       do k=1,2\r
 925          dcu=dcu0\r
 926          if(k.eq.1)then\r
 927             nseq0=nseq1\r
 928             r_min=nseq1/3.0     !minimum fragment, in case too small protein\r
 929          else\r
 930             nseq0=nseq2\r
 931             r_min=nseq2/3.0     !minimum fragment, in case too small protein\r
 932          endif\r
 933          if(r_min.gt.fra_min)r_min=fra_min\r
 934  20      nfr=1                  !number of fragments\r
 935          j=1                    !number of residue at nf-fragment\r
 936          ifr2(k,nfr,j)=1        !what residue\r
 937          Lfr2(k,nfr)=j          !length of the fragment\r
 938          do i=2,nseq0\r
 939             dis=diszy(k-1,i-1,i)\r
 940             contin=.false.\r
 941             if(dcu.gt.dcu0)then\r
 942                if(dis.lt.dcu)then\r
 943                   if(dis.gt.dcu_min)then\r
 944                      contin=.true.\r
 945                   endif\r
 946                endif\r
 947             elseif(mm(k,i).eq.(mm(k,i-1)+1))then\r
 948                if(dis.lt.dcu)then\r
 949                   if(dis.gt.dcu_min)then\r
 950                      contin=.true.\r
 951                   endif\r
 952                endif\r
 953             endif\r
 954             if(contin)then\r
 955                j=j+1\r
 956                ifr2(k,nfr,j)=i\r
 957                Lfr2(k,nfr)=j\r
 958             else\r
 959                nfr=nfr+1\r
 960                j=1\r
 961                ifr2(k,nfr,j)=i\r
 962                Lfr2(k,nfr)=j\r
 963             endif\r
 964          enddo\r
 965          Lfr_max=0\r
 966          i_fr2(k)=1             !ID of the maximum piece\r
 967          do i=1,nfr\r
 968             if(Lfr_max.lt.Lfr2(k,i))then\r
 969                Lfr_max=Lfr2(k,i)\r
 970                i_fr2(k)=i\r
 971             endif\r
 972          enddo\r
 973          if(Lfr_max.lt.r_min)then\r
 974             dcu=1.1*dcu\r
 975             goto 20\r
 976          endif\r
 977       enddo\r
 978 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
 979       \r
 980 ccc   select what piece will be used (this may araise ansysmetry, but\r
 981 ccc   only when L1=L2 and Lfr1=Lfr2 and L1 ne Lfr1\r
 982 ccc   if L1=Lfr1 and L2=Lfr2 (normal proteins), it will be the same as initial1\r
 983       mark=1\r
 984       if(Lfr2(1,i_fr2(1)).lt.Lfr2(2,i_fr2(2)))then\r
 985          mark=1\r
 986       elseif(Lfr2(1,i_fr2(1)).gt.Lfr2(2,i_fr2(2)))then\r
 987          mark=2\r
 988       else                      !Lfr1=Lfr2\r
 989          if(nseq1.le.nseq2)then\r
 990             mark=1\r
 991          else\r
 992             mark=2\r
 993          endif\r
 994       endif\r
 995 ccc   \r
 996       L_fr=Lfr2(mark,i_fr2(mark))\r
 997       do i=1,L_fr\r
 998          ifr(i)=ifr2(mark,i_fr2(mark),i)\r
 999       enddo\r
1000 ccc   \r
1001       if(mark.eq.1)then         !non-redundant to get_initial1\r
1002          nseq0=nseq1\r
1003       else\r
1004          nseq0=nseq2\r
1005       endif\r
1006       if(L_fr.eq.nseq0)then\r
1007          n1=int(nseq0*0.1)      !0\r
1008          n2=int(nseq0*0.89)     !2\r
1009          j=0\r
1010          do i=n1,n2\r
1011             j=j+1\r
1012             ifr(j)=ifr(n1+j)\r
1013          enddo\r
1014          L_fr=j\r
1015       endif\r
1016       \r
1017 ccc   get initial ------------->\r
1018       if(mark.eq.1)then    !nseq1 as the smallest one\r
1019          nseq1_=L_fr\r
1020          aL=min(nseq1_,nseq2)\r
1021          idel=aL/2.5            !minimum size of considered fragment\r
1022          if(idel.le.fra_min1)idel=fra_min1\r
1023          n1=-nseq2+idel         !shift1\r
1024          n2=nseq1_-idel         !shift2\r
1025          GL_max=0\r
1026          do ishift=n1,n2\r
1027             L=0\r
1028             do j=1,nseq2\r
1029                i=j+ishift\r
1030                if(i.ge.1.and.i.le.nseq1_)then\r
1031                   L=L+1\r
1032                   invmap(j)=ifr(i)\r
1033                else\r
1034                   invmap(j)=-1\r
1035                endif\r
1036             enddo\r
1037             if(L.ge.idel)then\r
1038                call get_GL(GL)\r
1039                if(GL.gt.GL_max)then\r
1040                   GL_max=GL\r
1041                   do i=1,nseq2\r
1042                      invmap_i(i)=invmap(i)\r
1043                   enddo\r
1044                endif\r
1045             endif\r
1046          enddo\r
1047       else                      !@@@@@@@@@@@@@@@@@@@@\r
1048          nseq2_=L_fr\r
1049          aL=min(nseq1,nseq2_)\r
1050          idel=aL/2.5            !minimum size of considered fragment\r
1051          if(idel.le.fra_min1)idel=fra_min1\r
1052          n1=-nseq2_+idel\r
1053          n2=nseq1-idel\r
1054          GL_max=0\r
1055          do ishift=n1,n2\r
1056             L=0\r
1057             do j=1,nseq2\r
1058                invmap(j)=-1\r
1059             enddo\r
1060             do j=1,nseq2_\r
1061                i=j+ishift\r
1062                if(i.ge.1.and.i.le.nseq1)then\r
1063                   L=L+1\r
1064                   invmap(ifr(j))=i\r
1065                endif\r
1066             enddo\r
1067             if(L.ge.idel)then\r
1068                call get_GL(GL)\r
1069                if(GL.gt.GL_max)then\r
1070                   GL_max=GL\r
1071                   do i=1,nseq2\r
1072                      invmap_i(i)=invmap(i)\r
1073                   enddo\r
1074                endif\r
1075             endif\r
1076          enddo\r
1077       endif\r
1078       \r
1079       return\r
1080       end\r
1081 \r
1082 **************************************************************\r
1083 *     smooth the secondary structure assignment\r
1084 **************************************************************\r
1085       subroutine smooth\r
1086       PARAMETER(nmax=5000)\r
1087       common/sec/isec(nmax),jsec(nmax)\r
1088       common/length/nseq1,nseq2\r
1089 \r
1090 ***   smooth single -------------->\r
1091 ***   --x-- => -----\r
1092       do i=1,nseq1\r
1093          if(isec(i).eq.2.or.isec(i).eq.4)then\r
1094             j=isec(i)\r
1095             if(isec(i-2).ne.j)then\r
1096                if(isec(i-1).ne.j)then\r
1097                   if(isec(i+1).ne.j)then\r
1098                      if(isec(i+1).ne.j)then\r
1099                         isec(i)=1\r
1100                      endif\r
1101                   endif\r
1102                endif\r
1103             endif\r
1104          endif\r
1105       enddo\r
1106       do i=1,nseq2\r
1107          if(jsec(i).eq.2.or.jsec(i).eq.4)then\r
1108             j=jsec(i)\r
1109             if(jsec(i-2).ne.j)then\r
1110                if(jsec(i-1).ne.j)then\r
1111                   if(jsec(i+1).ne.j)then\r
1112                      if(jsec(i+1).ne.j)then\r
1113                         jsec(i)=1\r
1114                      endif\r
1115                   endif\r
1116                endif\r
1117             endif\r
1118          endif\r
1119       enddo\r
1120 \r
1121 ***   smooth double -------------->\r
1122 ***   --xx-- => ------\r
1123       do i=1,nseq1\r
1124          if(isec(i).ne.2)then\r
1125          if(isec(i+1).ne.2)then\r
1126          if(isec(i+2).eq.2)then\r
1127          if(isec(i+3).eq.2)then\r
1128          if(isec(i+4).ne.2)then\r
1129          if(isec(i+5).ne.2)then\r
1130             isec(i+2)=1\r
1131             isec(i+3)=1\r
1132          endif\r
1133          endif\r
1134          endif\r
1135          endif\r
1136          endif\r
1137          endif\r
1138 \r
1139          if(isec(i).ne.4)then\r
1140          if(isec(i+1).ne.4)then\r
1141          if(isec(i+2).eq.4)then\r
1142          if(isec(i+3).eq.4)then\r
1143          if(isec(i+4).ne.4)then\r
1144          if(isec(i+5).ne.4)then\r
1145             isec(i+2)=1\r
1146             isec(i+3)=1\r
1147          endif\r
1148          endif\r
1149          endif\r
1150          endif\r
1151          endif\r
1152          endif\r
1153       enddo\r
1154       do i=1,nseq2\r
1155          if(jsec(i).ne.2)then\r
1156          if(jsec(i+1).ne.2)then\r
1157          if(jsec(i+2).eq.2)then\r
1158          if(jsec(i+3).eq.2)then\r
1159          if(jsec(i+4).ne.2)then\r
1160          if(jsec(i+5).ne.2)then\r
1161             jsec(i+2)=1\r
1162             jsec(i+3)=1\r
1163          endif\r
1164          endif\r
1165          endif\r
1166          endif\r
1167          endif\r
1168          endif\r
1169 \r
1170          if(jsec(i).ne.4)then\r
1171          if(jsec(i+1).ne.4)then\r
1172          if(jsec(i+2).eq.4)then\r
1173          if(jsec(i+3).eq.4)then\r
1174          if(jsec(i+4).ne.4)then\r
1175          if(jsec(i+5).ne.4)then\r
1176             jsec(i+2)=1\r
1177             jsec(i+3)=1\r
1178          endif\r
1179          endif\r
1180          endif\r
1181          endif\r
1182          endif\r
1183          endif\r
1184       enddo\r
1185 \r
1186 ***   connect -------------->\r
1187 ***   x-x => xxx\r
1188       do i=1,nseq1\r
1189          if(isec(i).eq.2)then\r
1190          if(isec(i+1).ne.2)then\r
1191          if(isec(i+2).eq.2)then\r
1192             isec(i+1)=2\r
1193          endif\r
1194          endif\r
1195          endif\r
1196 \r
1197          if(isec(i).eq.4)then\r
1198          if(isec(i+1).ne.4)then\r
1199          if(isec(i+2).eq.4)then\r
1200             isec(i+1)=4\r
1201          endif\r
1202          endif\r
1203          endif\r
1204       enddo\r
1205       do i=1,nseq2\r
1206          if(jsec(i).eq.2)then\r
1207          if(jsec(i+1).ne.2)then\r
1208          if(jsec(i+2).eq.2)then\r
1209             jsec(i+1)=2\r
1210          endif\r
1211          endif\r
1212          endif\r
1213 \r
1214          if(jsec(i).eq.4)then\r
1215          if(jsec(i+1).ne.4)then\r
1216          if(jsec(i+2).eq.4)then\r
1217             jsec(i+1)=4\r
1218          endif\r
1219          endif\r
1220          endif\r
1221       enddo\r
1222 \r
1223       return\r
1224       end\r
1225 \r
1226 *************************************************************\r
1227 *     assign secondary structure:\r
1228 *************************************************************\r
1229       function diszy(i,i1,i2)\r
1230       PARAMETER(nmax=5000)\r
1231       COMMON/BACKBONE/XA(3,nmax,0:1)\r
1232       diszy=sqrt((xa(1,i1,i)-xa(1,i2,i))**2\r
1233      &     +(xa(2,i1,i)-xa(2,i2,i))**2\r
1234      &     +(xa(3,i1,i)-xa(3,i2,i))**2)\r
1235       return\r
1236       end\r
1237 \r
1238 *************************************************************\r
1239 *     assign secondary structure:\r
1240 *************************************************************\r
1241       function make_sec(dis13,dis14,dis15,dis24,dis25,dis35)\r
1242       make_sec=1\r
1243       delta=2.1\r
1244       if(abs(dis15-6.37).lt.delta)then\r
1245          if(abs(dis14-5.18).lt.delta)then\r
1246             if(abs(dis25-5.18).lt.delta)then\r
1247                if(abs(dis13-5.45).lt.delta)then\r
1248                   if(abs(dis24-5.45).lt.delta)then\r
1249                      if(abs(dis35-5.45).lt.delta)then\r
1250                         make_sec=2 !helix\r
1251                         return\r
1252                      endif\r
1253                   endif\r
1254                endif\r
1255             endif\r
1256          endif\r
1257       endif\r
1258       delta=1.42\r
1259       if(abs(dis15-13).lt.delta)then\r
1260          if(abs(dis14-10.4).lt.delta)then\r
1261             if(abs(dis25-10.4).lt.delta)then\r
1262                if(abs(dis13-6.1).lt.delta)then\r
1263                   if(abs(dis24-6.1).lt.delta)then\r
1264                      if(abs(dis35-6.1).lt.delta)then\r
1265                         make_sec=4 !strand\r
1266                         return\r
1267                      endif\r
1268                   endif\r
1269                endif\r
1270             endif\r
1271          endif\r
1272       endif\r
1273       if(dis15.lt.8)then\r
1274          make_sec=3\r
1275       endif\r
1276 \r
1277       return\r
1278       end\r
1279 \r
1280 ****************************************************************\r
1281 *     quickly calculate TM-score with given invmap(i) in 3 iterations\r
1282 ****************************************************************\r
1283       subroutine get_GL(GL)\r
1284       PARAMETER(nmax=5000)\r
1285       common/length/nseq1,nseq2\r
1286       COMMON/BACKBONE/XA(3,nmax,0:1)\r
1287       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
1288       common/zscore/zrms,n_al,rmsd_al\r
1289       common/d0/d0,anseq\r
1290       dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
1291       dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
1292       common/TM/TM,TMmax\r
1293       common/n1n2/n1(nmax),n2(nmax)\r
1294       common/d00/d00,d002\r
1295 \r
1296       dimension xo1(nmax),yo1(nmax),zo1(nmax)\r
1297       dimension xo2(nmax),yo2(nmax),zo2(nmax)\r
1298       dimension dis2(nmax)\r
1299 \r
1300 ccc   RMSD:\r
1301       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
1302       double precision u(3,3),t(3),rms,drms !armsd is real\r
1303       data w /nmax*1.0/\r
1304 ccc   \r
1305 \r
1306 c     calculate RMSD between aligned structures and rotate the structures -->\r
1307       n_al=0\r
1308       do j=1,NSEQ2\r
1309          i=invmap(j)            !j aligned to i\r
1310          if(i.gt.0)then\r
1311             n_al=n_al+1\r
1312             r_1(1,n_al)=xa(1,i,0)\r
1313             r_1(2,n_al)=xa(2,i,0)\r
1314             r_1(3,n_al)=xa(3,i,0)\r
1315             r_2(1,n_al)=xa(1,j,1)\r
1316             r_2(2,n_al)=xa(2,j,1)\r
1317             r_2(3,n_al)=xa(3,j,1)\r
1318             xo1(n_al)=xa(1,i,0)\r
1319             yo1(n_al)=xa(2,i,0)\r
1320             zo1(n_al)=xa(3,i,0)\r
1321             xo2(n_al)=xa(1,j,1)\r
1322             yo2(n_al)=xa(2,j,1)\r
1323             zo2(n_al)=xa(3,j,1)\r
1324          endif\r
1325       enddo\r
1326       call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1327       GL=0\r
1328       do i=1,n_al\r
1329          xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)\r
1330          yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)\r
1331          zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)\r
1332          dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2\r
1333          GL=GL+1/(1+dis2(i)/(d0**2))\r
1334       enddo\r
1335 ccc   for next iteration------------->\r
1336       d002t=d002\r
1337  21   j=0\r
1338       do i=1,n_al\r
1339          if(dis2(i).le.d002t)then\r
1340             j=j+1\r
1341             r_1(1,j)=xo1(i)\r
1342             r_1(2,j)=yo1(i)\r
1343             r_1(3,j)=zo1(i)\r
1344             r_2(1,j)=xo2(i)\r
1345             r_2(2,j)=yo2(i)\r
1346             r_2(3,j)=zo2(i)\r
1347          endif\r
1348       enddo\r
1349       if(j.lt.3.and.n_al.gt.3)then\r
1350          d002t=d002t+.5\r
1351          goto 21\r
1352       endif\r
1353       L=j\r
1354       call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1355       G2=0\r
1356       do i=1,n_al\r
1357          xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)\r
1358          yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)\r
1359          zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)\r
1360          dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2\r
1361          G2=G2+1/(1+dis2(i)/(d0**2))\r
1362       enddo\r
1363 ccc   for next iteration------------->\r
1364       d002t=d002+1\r
1365  22   j=0\r
1366       do i=1,n_al\r
1367          if(dis2(i).le.d002t)then\r
1368             j=j+1\r
1369             r_1(1,j)=xo1(i)\r
1370             r_1(2,j)=yo1(i)\r
1371             r_1(3,j)=zo1(i)\r
1372             r_2(1,j)=xo2(i)\r
1373             r_2(2,j)=yo2(i)\r
1374             r_2(3,j)=zo2(i)\r
1375          endif\r
1376       enddo\r
1377       if(j.lt.3.and.n_al.gt.3)then\r
1378          d002t=d002t+.5\r
1379          goto 22\r
1380       endif\r
1381       L=j\r
1382       call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1383       G3=0\r
1384       do i=1,n_al\r
1385          xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)\r
1386          yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)\r
1387          zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)\r
1388          dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2\r
1389          G3=G3+1/(1+dis2(i)/(d0**2))\r
1390       enddo\r
1391       if(G2.gt.GL)GL=G2\r
1392       if(G3.gt.GL)GL=G3\r
1393 \r
1394 c^^^^^^^^^^^^^^^^ GL done ^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
1395       return\r
1396       end\r
1397 \r
1398 ****************************************************************\r
1399 *     with invmap(i) calculate TM-score and martix score(i,j) for rotation \r
1400 ****************************************************************\r
1401       subroutine get_score\r
1402       PARAMETER(nmax=5000)\r
1403       common/length/nseq1,nseq2\r
1404       COMMON/BACKBONE/XA(3,nmax,0:1)\r
1405       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
1406       common/zscore/zrms,n_al,rmsd_al\r
1407       common/d0/d0,anseq\r
1408       dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
1409       dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
1410       common/TM/TM,TMmax\r
1411       common/n1n2/n1(nmax),n2(nmax)\r
1412 \r
1413 ccc   RMSD:\r
1414       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
1415       double precision u(3,3),t(3),rms,drms !armsd is real\r
1416       data w /nmax*1.0/\r
1417 ccc   \r
1418 \r
1419 c     calculate RMSD between aligned structures and rotate the structures -->\r
1420       n_al=0\r
1421       do j=1,NSEQ2\r
1422          i=invmap(j)            !j aligned to i\r
1423          if(i.gt.0)then\r
1424             n_al=n_al+1\r
1425 ccc   for TM-score:\r
1426             xtm1(n_al)=xa(1,i,0) !for TM-score\r
1427             ytm1(n_al)=xa(2,i,0)\r
1428             ztm1(n_al)=xa(3,i,0)\r
1429             xtm2(n_al)=xa(1,j,1)\r
1430             ytm2(n_al)=xa(2,j,1)\r
1431             ztm2(n_al)=xa(3,j,1)\r
1432 ccc   for rotation matrix:\r
1433             r_1(1,n_al)=xa(1,i,0)\r
1434             r_1(2,n_al)=xa(2,i,0)\r
1435             r_1(3,n_al)=xa(3,i,0)\r
1436          endif\r
1437       enddo\r
1438 ***   calculate TM-score for the given alignment----------->\r
1439       d0_input=d0\r
1440       call TMscore8_search(d0_input,n_al,xtm1,ytm1,ztm1,n1,\r
1441      &     n_al,xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !simplified search engine\r
1442       TM=TM*n_al/anseq          !TM-score\r
1443 ***   calculate score matrix score(i,j)------------------>\r
1444       do i=1,n_al\r
1445          r_2(1,i)=xtm1(i)\r
1446          r_2(2,i)=ytm1(i)\r
1447          r_2(3,i)=ztm1(i)\r
1448       enddo\r
1449       call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1450       do i=1,nseq1\r
1451          xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)\r
1452          yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)\r
1453          zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)\r
1454          do j=1,nseq2\r
1455             dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2\r
1456             score(i,j)=1/(1+dd/d0**2)\r
1457          enddo\r
1458       enddo\r
1459       \r
1460 c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
1461       return\r
1462       end\r
1463 \r
1464 ****************************************************************\r
1465 *     with invmap(i) calculate score(i,j) using RMSD rotation \r
1466 ****************************************************************\r
1467       subroutine get_score1\r
1468       PARAMETER(nmax=5000)\r
1469       common/length/nseq1,nseq2\r
1470       COMMON/BACKBONE/XA(3,nmax,0:1)\r
1471       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
1472       common/zscore/zrms,n_al,rmsd_al\r
1473       common/d0/d0,anseq\r
1474       common/d0min/d0_min\r
1475       dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
1476       dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
1477       common/TM/TM,TMmax\r
1478       common/n1n2/n1(nmax),n2(nmax)\r
1479 \r
1480 ccc   RMSD:\r
1481       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
1482       double precision u(3,3),t(3),rms,drms !armsd is real\r
1483       data w /nmax*1.0/\r
1484 ccc   \r
1485 \r
1486 c     calculate RMSD between aligned structures and rotate the structures -->\r
1487       n_al=0\r
1488       do j=1,NSEQ2\r
1489          i=invmap(j)            !j aligned to i\r
1490          if(i.gt.0)then\r
1491             n_al=n_al+1\r
1492 ccc   for rotation matrix:\r
1493             r_1(1,n_al)=xa(1,i,0)\r
1494             r_1(2,n_al)=xa(2,i,0)\r
1495             r_1(3,n_al)=xa(3,i,0)\r
1496             r_2(1,n_al)=xa(1,j,1)\r
1497             r_2(2,n_al)=xa(2,j,1)\r
1498             r_2(3,n_al)=xa(3,j,1)\r
1499          endif\r
1500       enddo\r
1501 ***   calculate score matrix score(i,j)------------------>\r
1502       call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1503       d01=d0+1.5\r
1504       if(d01.lt.d0_min)d01=d0_min\r
1505       d02=d01*d01\r
1506       do i=1,nseq1\r
1507          xx=t(1)+u(1,1)*xa(1,i,0)+u(1,2)*xa(2,i,0)+u(1,3)*xa(3,i,0)\r
1508          yy=t(2)+u(2,1)*xa(1,i,0)+u(2,2)*xa(2,i,0)+u(2,3)*xa(3,i,0)\r
1509          zz=t(3)+u(3,1)*xa(1,i,0)+u(3,2)*xa(2,i,0)+u(3,3)*xa(3,i,0)\r
1510          do j=1,nseq2\r
1511             dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2\r
1512             score(i,j)=1/(1+dd/d02)\r
1513          enddo\r
1514       enddo\r
1515 \r
1516 c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
1517       return\r
1518       end\r
1519 \r
1520 \r
1521 *************************************************************************\r
1522 *************************************************************************\r
1523 *     This is a subroutine to compare two structures and find the \r
1524 *     superposition that has the maximum TM-score.\r
1525 *\r
1526 *     L1--Length of the first structure\r
1527 *     (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure\r
1528 *     n1(i)--Residue sequence number of i'th residue at the first structure\r
1529 *     L2--Length of the second structure\r
1530 *     (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure\r
1531 *     n2(i)--Residue sequence number of i'th residue at the second structure\r
1532 *     TM--TM-score of the comparison\r
1533 *     Rcomm--RMSD of two structures in the common aligned residues\r
1534 *     Lcomm--Length of the common aligned regions\r
1535 *\r
1536 *     Note: \r
1537 *     1, Always put native as the second structure, by which TM-score\r
1538 *        is normalized.\r
1539 *     2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after\r
1540 *        TM-score superposition.\r
1541 *************************************************************************\r
1542 *************************************************************************\r
1543 *** dis<8, simplified search engine\r
1544       subroutine TMscore8_search(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,\r
1545      &     TM,Rcomm,Lcomm)\r
1546       PARAMETER(nmax=5000)\r
1547       common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)\r
1548       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
1549       common/para/d,d0\r
1550       common/d0min/d0_min\r
1551       common/align/n_ali,iA(nmax),iB(nmax)\r
1552       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
1553       dimension k_ali(nmax),k_ali0(nmax)\r
1554       dimension L_ini(100),iq(nmax)\r
1555       common/scores/score\r
1556       double precision score,score_max\r
1557       dimension xa(nmax),ya(nmax),za(nmax)\r
1558       dimension iL0(nmax)\r
1559 \r
1560       dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)\r
1561       dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)\r
1562 \r
1563 ccc   RMSD:\r
1564       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
1565       double precision u(3,3),t(3),rms,drms !armsd is real\r
1566       data w /nmax*1.0/\r
1567 ccc   \r
1568 \r
1569 ********* convert input data ***************************\r
1570 *     because L1=L2 in this special case---------->\r
1571       nseqA=L1\r
1572       nseqB=L2\r
1573       do i=1,nseqA\r
1574          xa(i)=x1(i)\r
1575          ya(i)=y1(i)\r
1576          za(i)=z1(i)\r
1577          nresA(i)=n1(i)\r
1578          xb(i)=x2(i)\r
1579          yb(i)=y2(i)\r
1580          zb(i)=z2(i)\r
1581          nresB(i)=n2(i)\r
1582          iA(i)=i\r
1583          iB(i)=i\r
1584       enddo\r
1585       n_ali=L1                  !number of aligned residues\r
1586       Lcomm=L1\r
1587 \r
1588 ************/////\r
1589 *     parameters:\r
1590 *****************\r
1591 ***   d0------------->\r
1592       d0=dx\r
1593       if(d0.lt.d0_min)d0=d0_min\r
1594 ***   d0_search ----->\r
1595       d0_search=d0\r
1596       if(d0_search.gt.8)d0_search=8\r
1597       if(d0_search.lt.4.5)d0_search=4.5\r
1598 ***   iterative parameters ----->\r
1599       n_it=20                   !maximum number of iterations\r
1600       d_output=5                !for output alignment\r
1601       n_init_max=6              !maximum number of L_init\r
1602       n_init=0\r
1603       L_ini_min=4\r
1604       if(n_ali.lt.4)L_ini_min=n_ali\r
1605       do i=1,n_init_max-1\r
1606          n_init=n_init+1\r
1607          L_ini(n_init)=n_ali/2**(n_init-1)\r
1608          if(L_ini(n_init).le.L_ini_min)then\r
1609             L_ini(n_init)=L_ini_min\r
1610             goto 402\r
1611          endif\r
1612       enddo\r
1613       n_init=n_init+1\r
1614       L_ini(n_init)=L_ini_min\r
1615  402  continue\r
1616 \r
1617 ******************************************************************\r
1618 *     find the maximum score starting from local structures superposition\r
1619 ******************************************************************\r
1620       score_max=-1              !TM-score\r
1621       do 333 i_init=1,n_init\r
1622         L_init=L_ini(i_init)\r
1623         iL_max=n_ali-L_init+1\r
1624         k=0\r
1625         do i=1,iL_max,40        !this is the simplification!\r
1626           k=k+1\r
1627           iL0(k)=i\r
1628         enddo\r
1629         if(iL0(k).lt.iL_max)then\r
1630           k=k+1\r
1631           iL0(k)=iL_max\r
1632         endif\r
1633         n_shift=k\r
1634         do 300 i_shift=1,n_shift\r
1635            iL=iL0(i_shift)\r
1636            LL=0\r
1637            ka=0\r
1638            do i=1,L_init\r
1639               k=iL+i-1          ![1,n_ali] common aligned\r
1640               r_1(1,i)=xa(iA(k))\r
1641               r_1(2,i)=ya(iA(k))\r
1642               r_1(3,i)=za(iA(k))\r
1643               r_2(1,i)=xb(iB(k))\r
1644               r_2(2,i)=yb(iB(k))\r
1645               r_2(3,i)=zb(iB(k))\r
1646               LL=LL+1\r
1647               ka=ka+1\r
1648               k_ali(ka)=k\r
1649            enddo\r
1650            call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1651            if(i_init.eq.1)then  !global superposition\r
1652               armsd=dsqrt(rms/LL)\r
1653               Rcomm=armsd\r
1654            endif\r
1655            do j=1,nseqA\r
1656               xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
1657               yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
1658               zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
1659            enddo\r
1660            d=d0_search-1\r
1661            call score_fun8       !init, get scores, n_cut+i_ali(i) for iteration\r
1662            if(score_max.lt.score)then\r
1663               score_max=score\r
1664               ka0=ka\r
1665               do i=1,ka0\r
1666                  k_ali0(i)=k_ali(i)\r
1667               enddo\r
1668            endif\r
1669 ***   iteration for extending ---------------------------------->\r
1670            d=d0_search+1\r
1671            do 301 it=1,n_it\r
1672               LL=0\r
1673               ka=0\r
1674               do i=1,n_cut\r
1675                  m=i_ali(i)     ![1,n_ali]\r
1676                  r_1(1,i)=xa(iA(m))\r
1677                  r_1(2,i)=ya(iA(m))\r
1678                  r_1(3,i)=za(iA(m))\r
1679                  r_2(1,i)=xb(iB(m))\r
1680                  r_2(2,i)=yb(iB(m))\r
1681                  r_2(3,i)=zb(iB(m))\r
1682                  ka=ka+1\r
1683                  k_ali(ka)=m\r
1684                  LL=LL+1\r
1685               enddo\r
1686               call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1687               do j=1,nseqA\r
1688                  xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
1689                  yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
1690                  zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
1691               enddo\r
1692               call score_fun8    !get scores, n_cut+i_ali(i) for iteration\r
1693               if(score_max.lt.score)then\r
1694                  score_max=score\r
1695                  ka0=ka\r
1696                  do i=1,ka\r
1697                     k_ali0(i)=k_ali(i)\r
1698                  enddo\r
1699               endif\r
1700               if(it.eq.n_it)goto 302\r
1701               if(n_cut.eq.ka)then\r
1702                  neq=0\r
1703                  do i=1,n_cut\r
1704                     if(i_ali(i).eq.k_ali(i))neq=neq+1\r
1705                  enddo\r
1706                  if(n_cut.eq.neq)goto 302\r
1707               endif\r
1708  301       continue             !for iteration\r
1709  302       continue\r
1710  300    continue                !for shift\r
1711  333  continue                  !for initial length, L_ali/M\r
1712 \r
1713 ******** return the final rotation ****************\r
1714       LL=0\r
1715       do i=1,ka0\r
1716          m=k_ali0(i)            !record of the best alignment\r
1717          r_1(1,i)=xa(iA(m))\r
1718          r_1(2,i)=ya(iA(m))\r
1719          r_1(3,i)=za(iA(m))\r
1720          r_2(1,i)=xb(iB(m))\r
1721          r_2(2,i)=yb(iB(m))\r
1722          r_2(3,i)=zb(iB(m))\r
1723          LL=LL+1\r
1724       enddo\r
1725       call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1726       do j=1,nseqA\r
1727          x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
1728          y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
1729          z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
1730       enddo\r
1731       TM=score_max\r
1732 \r
1733 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
1734       return\r
1735       END\r
1736 \r
1737 *************************************************************************\r
1738 *************************************************************************\r
1739 *     This is a subroutine to compare two structures and find the \r
1740 *     superposition that has the maximum TM-score.\r
1741 *\r
1742 *     L1--Length of the first structure\r
1743 *     (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure\r
1744 *     n1(i)--Residue sequence number of i'th residue at the first structure\r
1745 *     L2--Length of the second structure\r
1746 *     (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure\r
1747 *     n2(i)--Residue sequence number of i'th residue at the second structure\r
1748 *     TM--TM-score of the comparison\r
1749 *     Rcomm--RMSD of two structures in the common aligned residues\r
1750 *     Lcomm--Length of the common aligned regions\r
1751 *\r
1752 *     Note: \r
1753 *     1, Always put native as the second structure, by which TM-score\r
1754 *        is normalized.\r
1755 *     2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after\r
1756 *        TM-score superposition.\r
1757 *************************************************************************\r
1758 *************************************************************************\r
1759 ***   dis<8, but same search engine\r
1760       subroutine TMscore8(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,\r
1761      &     TM,Rcomm,Lcomm)\r
1762       PARAMETER(nmax=5000)\r
1763       common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)\r
1764       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
1765       common/para/d,d0\r
1766       common/d0min/d0_min\r
1767       common/align/n_ali,iA(nmax),iB(nmax)\r
1768       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
1769       dimension k_ali(nmax),k_ali0(nmax)\r
1770       dimension L_ini(100),iq(nmax)\r
1771       common/scores/score\r
1772       double precision score,score_max\r
1773       dimension xa(nmax),ya(nmax),za(nmax)\r
1774 \r
1775       dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)\r
1776       dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)\r
1777 \r
1778 ccc   RMSD:\r
1779       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
1780       double precision u(3,3),t(3),rms,drms !armsd is real\r
1781       data w /nmax*1.0/\r
1782 ccc   \r
1783 \r
1784 ********* convert input data ***************************\r
1785 *     because L1=L2 in this special case---------->\r
1786       nseqA=L1\r
1787       nseqB=L2\r
1788       do i=1,nseqA\r
1789          xa(i)=x1(i)\r
1790          ya(i)=y1(i)\r
1791          za(i)=z1(i)\r
1792          nresA(i)=n1(i)\r
1793          xb(i)=x2(i)\r
1794          yb(i)=y2(i)\r
1795          zb(i)=z2(i)\r
1796          nresB(i)=n2(i)\r
1797          iA(i)=i\r
1798          iB(i)=i\r
1799       enddo\r
1800       n_ali=L1                  !number of aligned residues\r
1801       Lcomm=L1\r
1802 \r
1803 ************/////\r
1804 *     parameters:\r
1805 *****************\r
1806 ***   d0------------->\r
1807       d0=dx\r
1808       if(d0.lt.d0_min)d0=d0_min\r
1809 ***   d0_search ----->\r
1810       d0_search=d0\r
1811       if(d0_search.gt.8)d0_search=8\r
1812       if(d0_search.lt.4.5)d0_search=4.5\r
1813 ***   iterative parameters ----->\r
1814       n_it=20                   !maximum number of iterations\r
1815       d_output=5                !for output alignment\r
1816       n_init_max=6              !maximum number of L_init\r
1817       n_init=0\r
1818       L_ini_min=4\r
1819       if(n_ali.lt.4)L_ini_min=n_ali\r
1820       do i=1,n_init_max-1\r
1821          n_init=n_init+1\r
1822          L_ini(n_init)=n_ali/2**(n_init-1)\r
1823          if(L_ini(n_init).le.L_ini_min)then\r
1824             L_ini(n_init)=L_ini_min\r
1825             goto 402\r
1826          endif\r
1827       enddo\r
1828       n_init=n_init+1\r
1829       L_ini(n_init)=L_ini_min\r
1830  402  continue\r
1831 \r
1832 ******************************************************************\r
1833 *     find the maximum score starting from local structures superposition\r
1834 ******************************************************************\r
1835       score_max=-1              !TM-score\r
1836       do 333 i_init=1,n_init\r
1837         L_init=L_ini(i_init)\r
1838         iL_max=n_ali-L_init+1\r
1839         do 300 iL=1,iL_max    !on aligned residues, [1,nseqA]\r
1840            LL=0\r
1841            ka=0\r
1842            do i=1,L_init\r
1843               k=iL+i-1          ![1,n_ali] common aligned\r
1844               r_1(1,i)=xa(iA(k))\r
1845               r_1(2,i)=ya(iA(k))\r
1846               r_1(3,i)=za(iA(k))\r
1847               r_2(1,i)=xb(iB(k))\r
1848               r_2(2,i)=yb(iB(k))\r
1849               r_2(3,i)=zb(iB(k))\r
1850               LL=LL+1\r
1851               ka=ka+1\r
1852               k_ali(ka)=k\r
1853            enddo\r
1854            call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1855            if(i_init.eq.1)then  !global superposition\r
1856               armsd=dsqrt(rms/LL)\r
1857               Rcomm=armsd\r
1858            endif\r
1859            do j=1,nseqA\r
1860               xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
1861               yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
1862               zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
1863            enddo\r
1864            d=d0_search-1\r
1865            call score_fun8       !init, get scores, n_cut+i_ali(i) for iteration\r
1866            if(score_max.lt.score)then\r
1867               score_max=score\r
1868               ka0=ka\r
1869               do i=1,ka0\r
1870                  k_ali0(i)=k_ali(i)\r
1871               enddo\r
1872            endif\r
1873 ***   iteration for extending ---------------------------------->\r
1874            d=d0_search+1\r
1875            do 301 it=1,n_it\r
1876               LL=0\r
1877               ka=0\r
1878               do i=1,n_cut\r
1879                  m=i_ali(i)     ![1,n_ali]\r
1880                  r_1(1,i)=xa(iA(m))\r
1881                  r_1(2,i)=ya(iA(m))\r
1882                  r_1(3,i)=za(iA(m))\r
1883                  r_2(1,i)=xb(iB(m))\r
1884                  r_2(2,i)=yb(iB(m))\r
1885                  r_2(3,i)=zb(iB(m))\r
1886                  ka=ka+1\r
1887                  k_ali(ka)=m\r
1888                  LL=LL+1\r
1889               enddo\r
1890               call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1891               do j=1,nseqA\r
1892                  xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
1893                  yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
1894                  zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
1895               enddo\r
1896               call score_fun8    !get scores, n_cut+i_ali(i) for iteration\r
1897               if(score_max.lt.score)then\r
1898                  score_max=score\r
1899                  ka0=ka\r
1900                  do i=1,ka\r
1901                     k_ali0(i)=k_ali(i)\r
1902                  enddo\r
1903               endif\r
1904               if(it.eq.n_it)goto 302\r
1905               if(n_cut.eq.ka)then\r
1906                  neq=0\r
1907                  do i=1,n_cut\r
1908                     if(i_ali(i).eq.k_ali(i))neq=neq+1\r
1909                  enddo\r
1910                  if(n_cut.eq.neq)goto 302\r
1911               endif\r
1912  301       continue             !for iteration\r
1913  302       continue\r
1914  300    continue                !for shift\r
1915  333  continue                  !for initial length, L_ali/M\r
1916 \r
1917 ******** return the final rotation ****************\r
1918       LL=0\r
1919       do i=1,ka0\r
1920          m=k_ali0(i)            !record of the best alignment\r
1921          r_1(1,i)=xa(iA(m))\r
1922          r_1(2,i)=ya(iA(m))\r
1923          r_1(3,i)=za(iA(m))\r
1924          r_2(1,i)=xb(iB(m))\r
1925          r_2(2,i)=yb(iB(m))\r
1926          r_2(3,i)=zb(iB(m))\r
1927          LL=LL+1\r
1928       enddo\r
1929       call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
1930       do j=1,nseqA\r
1931          x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
1932          y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
1933          z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
1934       enddo\r
1935       TM=score_max\r
1936 \r
1937 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
1938       return\r
1939       END\r
1940 \r
1941 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
1942 c     1, collect those residues with dis<d;\r
1943 c     2, calculate score_GDT, score_maxsub, score_TM\r
1944 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
1945       subroutine score_fun8\r
1946       PARAMETER(nmax=5000)\r
1947 \r
1948       common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)\r
1949       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
1950       common/para/d,d0\r
1951       common/align/n_ali,iA(nmax),iB(nmax)\r
1952       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
1953       common/scores/score\r
1954       double precision score,score_max\r
1955       common/d8/d8\r
1956 \r
1957       d_tmp=d\r
1958  21   n_cut=0                   !number of residue-pairs dis<d, for iteration\r
1959       score_sum=0               !TMscore\r
1960       do k=1,n_ali\r
1961          i=iA(k)                ![1,nseqA] reoder number of structureA\r
1962          j=iB(k)                ![1,nseqB]\r
1963          dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)\r
1964          if(dis.lt.d_tmp)then\r
1965             n_cut=n_cut+1\r
1966             i_ali(n_cut)=k      ![1,n_ali], mark the residue-pairs in dis<d\r
1967          endif\r
1968          if(dis.le.d8)then\r
1969             score_sum=score_sum+1/(1+(dis/d0)**2)\r
1970          endif\r
1971       enddo\r
1972       if(n_cut.lt.3.and.n_ali.gt.3)then\r
1973          d_tmp=d_tmp+.5\r
1974          goto 21\r
1975       endif\r
1976       score=score_sum/float(nseqB) !TM-score\r
1977 \r
1978       return\r
1979       end\r
1980 \r
1981 *************************************************************************\r
1982 *************************************************************************\r
1983 *     This is a subroutine to compare two structures and find the \r
1984 *     superposition that has the maximum TM-score.\r
1985 *\r
1986 *     L1--Length of the first structure\r
1987 *     (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure\r
1988 *     n1(i)--Residue sequence number of i'th residue at the first structure\r
1989 *     L2--Length of the second structure\r
1990 *     (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure\r
1991 *     n2(i)--Residue sequence number of i'th residue at the second structure\r
1992 *     TM--TM-score of the comparison\r
1993 *     Rcomm--RMSD of two structures in the common aligned residues\r
1994 *     Lcomm--Length of the common aligned regions\r
1995 *\r
1996 *     Note: \r
1997 *     1, Always put native as the second structure, by which TM-score\r
1998 *        is normalized.\r
1999 *     2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after\r
2000 *        TM-score superposition.\r
2001 *************************************************************************\r
2002 *************************************************************************\r
2003 ***  normal TM-score:\r
2004       subroutine TMscore(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,\r
2005      &     TM,Rcomm,Lcomm)\r
2006       PARAMETER(nmax=5000)\r
2007       common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)\r
2008       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
2009       common/para/d,d0\r
2010       common/d0min/d0_min\r
2011       common/align/n_ali,iA(nmax),iB(nmax)\r
2012       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
2013       dimension k_ali(nmax),k_ali0(nmax)\r
2014       dimension L_ini(100),iq(nmax)\r
2015       common/scores/score\r
2016       double precision score,score_max\r
2017       dimension xa(nmax),ya(nmax),za(nmax)\r
2018 \r
2019       dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)\r
2020       dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)\r
2021 \r
2022 ccc   RMSD:\r
2023       double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
2024       double precision u(3,3),t(3),rms,drms !armsd is real\r
2025       data w /nmax*1.0/\r
2026 ccc   \r
2027 \r
2028 ********* convert input data ***************************\r
2029 *     because L1=L2 in this special case---------->\r
2030       nseqA=L1\r
2031       nseqB=L2\r
2032       do i=1,nseqA\r
2033          xa(i)=x1(i)\r
2034          ya(i)=y1(i)\r
2035          za(i)=z1(i)\r
2036          nresA(i)=n1(i)\r
2037          xb(i)=x2(i)\r
2038          yb(i)=y2(i)\r
2039          zb(i)=z2(i)\r
2040          nresB(i)=n2(i)\r
2041          iA(i)=i\r
2042          iB(i)=i\r
2043       enddo\r
2044       n_ali=L1                  !number of aligned residues\r
2045       Lcomm=L1\r
2046 \r
2047 ************/////\r
2048 *     parameters:\r
2049 *****************\r
2050 ***   d0------------->\r
2051 c      d0=1.24*(nseqB-15)**(1.0/3.0)-1.8\r
2052       d0=dx\r
2053       if(d0.lt.d0_min)d0=d0_min\r
2054 ***   d0_search ----->\r
2055       d0_search=d0\r
2056       if(d0_search.gt.8)d0_search=8\r
2057       if(d0_search.lt.4.5)d0_search=4.5\r
2058 ***   iterative parameters ----->\r
2059       n_it=20                   !maximum number of iterations\r
2060       d_output=5                !for output alignment\r
2061       n_init_max=6              !maximum number of L_init\r
2062       n_init=0\r
2063       L_ini_min=4\r
2064       if(n_ali.lt.4)L_ini_min=n_ali\r
2065       do i=1,n_init_max-1\r
2066          n_init=n_init+1\r
2067          L_ini(n_init)=n_ali/2**(n_init-1)\r
2068          if(L_ini(n_init).le.L_ini_min)then\r
2069             L_ini(n_init)=L_ini_min\r
2070             goto 402\r
2071          endif\r
2072       enddo\r
2073       n_init=n_init+1\r
2074       L_ini(n_init)=L_ini_min\r
2075  402  continue\r
2076 \r
2077 ******************************************************************\r
2078 *     find the maximum score starting from local structures superposition\r
2079 ******************************************************************\r
2080       score_max=-1              !TM-score\r
2081       do 333 i_init=1,n_init\r
2082         L_init=L_ini(i_init)\r
2083         iL_max=n_ali-L_init+1\r
2084         do 300 iL=1,iL_max      !on aligned residues, [1,nseqA]\r
2085            LL=0\r
2086            ka=0\r
2087            do i=1,L_init\r
2088               k=iL+i-1          ![1,n_ali] common aligned\r
2089               r_1(1,i)=xa(iA(k))\r
2090               r_1(2,i)=ya(iA(k))\r
2091               r_1(3,i)=za(iA(k))\r
2092               r_2(1,i)=xb(iB(k))\r
2093               r_2(2,i)=yb(iB(k))\r
2094               r_2(3,i)=zb(iB(k))\r
2095               LL=LL+1\r
2096               ka=ka+1\r
2097               k_ali(ka)=k\r
2098            enddo\r
2099            call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
2100            if(i_init.eq.1)then  !global superposition\r
2101               armsd=dsqrt(rms/LL)\r
2102               Rcomm=armsd\r
2103            endif\r
2104            do j=1,nseqA\r
2105               xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
2106               yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
2107               zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
2108            enddo\r
2109            d=d0_search-1\r
2110            call score_fun       !init, get scores, n_cut+i_ali(i) for iteration\r
2111            if(score_max.lt.score)then\r
2112               score_max=score\r
2113               ka0=ka\r
2114               do i=1,ka0\r
2115                  k_ali0(i)=k_ali(i)\r
2116               enddo\r
2117            endif\r
2118 ***   iteration for extending ---------------------------------->\r
2119            d=d0_search+1\r
2120            do 301 it=1,n_it\r
2121               LL=0\r
2122               ka=0\r
2123               do i=1,n_cut\r
2124                  m=i_ali(i)     ![1,n_ali]\r
2125                  r_1(1,i)=xa(iA(m))\r
2126                  r_1(2,i)=ya(iA(m))\r
2127                  r_1(3,i)=za(iA(m))\r
2128                  r_2(1,i)=xb(iB(m))\r
2129                  r_2(2,i)=yb(iB(m))\r
2130                  r_2(3,i)=zb(iB(m))\r
2131                  ka=ka+1\r
2132                  k_ali(ka)=m\r
2133                  LL=LL+1\r
2134               enddo\r
2135               call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
2136               do j=1,nseqA\r
2137                  xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
2138                  yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
2139                  zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
2140               enddo\r
2141               call score_fun    !get scores, n_cut+i_ali(i) for iteration\r
2142               if(score_max.lt.score)then\r
2143                  score_max=score\r
2144                  ka0=ka\r
2145                  do i=1,ka\r
2146                     k_ali0(i)=k_ali(i)\r
2147                  enddo\r
2148               endif\r
2149               if(it.eq.n_it)goto 302\r
2150               if(n_cut.eq.ka)then\r
2151                  neq=0\r
2152                  do i=1,n_cut\r
2153                     if(i_ali(i).eq.k_ali(i))neq=neq+1\r
2154                  enddo\r
2155                  if(n_cut.eq.neq)goto 302\r
2156               endif\r
2157  301       continue             !for iteration\r
2158  302       continue\r
2159  300    continue                !for shift\r
2160  333  continue                  !for initial length, L_ali/M\r
2161 \r
2162 ******** return the final rotation ****************\r
2163       LL=0\r
2164       do i=1,ka0\r
2165          m=k_ali0(i)            !record of the best alignment\r
2166          r_1(1,i)=xa(iA(m))\r
2167          r_1(2,i)=ya(iA(m))\r
2168          r_1(3,i)=za(iA(m))\r
2169          r_2(1,i)=xb(iB(m))\r
2170          r_2(2,i)=yb(iB(m))\r
2171          r_2(3,i)=zb(iB(m))\r
2172          LL=LL+1\r
2173       enddo\r
2174       call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
2175       do j=1,nseqA\r
2176          x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
2177          y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
2178          z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
2179       enddo\r
2180       TM=score_max\r
2181 \r
2182 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
2183       return\r
2184       END\r
2185 \r
2186 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
2187 c     1, collect those residues with dis<d;\r
2188 c     2, calculate score_GDT, score_maxsub, score_TM\r
2189 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
2190       subroutine score_fun\r
2191       PARAMETER(nmax=5000)\r
2192 \r
2193       common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)\r
2194       common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
2195       common/para/d,d0\r
2196       common/align/n_ali,iA(nmax),iB(nmax)\r
2197       common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
2198       common/scores/score\r
2199       double precision score\r
2200 \r
2201       d_tmp=d\r
2202  21   n_cut=0                   !number of residue-pairs dis<d, for iteration\r
2203       score_sum=0               !TMscore\r
2204       do k=1,n_ali\r
2205          i=iA(k)                ![1,nseqA] reoder number of structureA\r
2206          j=iB(k)                ![1,nseqB]\r
2207          dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)\r
2208          if(dis.lt.d_tmp)then\r
2209             n_cut=n_cut+1\r
2210             i_ali(n_cut)=k      ![1,n_ali], mark the residue-pairs in dis<d\r
2211          endif\r
2212          score_sum=score_sum+1/(1+(dis/d0)**2)\r
2213       enddo\r
2214       if(n_cut.lt.3.and.n_ali.gt.3)then\r
2215          d_tmp=d_tmp+.5\r
2216          goto 21\r
2217       endif\r
2218       score=score_sum/float(nseqB) !TM-score\r
2219 \r
2220       return\r
2221       end\r
2222 \r
2223 ********************************************************************\r
2224 *     Dynamic programming for alignment.\r
2225 *     Input: score(i,j), and gap_open\r
2226 *     Output: invmap(j)\r
2227 *     \r
2228 *     Please note this subroutine is not a correct implementation of \r
2229 *     the N-W dynamic programming because the score tracks back only \r
2230 *     one layer of the matrix. This code was exploited in TM-align \r
2231 *     because it is about 1.5 times faster than a complete N-W code\r
2232 *     and does not influence much the final structure alignment result.\r
2233 ********************************************************************\r
2234       SUBROUTINE DP(NSEQ1,NSEQ2)\r
2235       PARAMETER(nmax=5000)\r
2236       LOGICAL*1 DIR\r
2237       common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
2238       dimension DIR(0:nmax,0:nmax),VAL(0:nmax,0:nmax)\r
2239       REAL H,V\r
2240       \r
2241 ***   initialize the matrix:\r
2242       val(0,0)=0\r
2243       do i=1,nseq1\r
2244         dir(i,0)=.false.\r
2245         val(i,0)=0\r
2246       enddo\r
2247       do j=1,nseq2\r
2248         dir(0,j)=.false.\r
2249         val(0,j)=0\r
2250         invmap(j)=-1\r
2251       enddo\r
2252 \r
2253 ***   decide matrix and path:\r
2254       DO j=1,NSEQ2\r
2255         DO i=1,NSEQ1\r
2256           D=VAL(i-1,j-1)+SCORE(i,j)\r
2257           H=VAL(i-1,j)\r
2258           if(DIR(i-1,j))H=H+GAP_OPEN\r
2259           V=VAL(i,j-1)\r
2260           if(DIR(i,j-1))V=V+GAP_OPEN\r
2261           \r
2262           IF((D.GE.H).AND.(D.GE.V)) THEN\r
2263             DIR(I,J)=.true.\r
2264             VAL(i,j)=D\r
2265           ELSE\r
2266             DIR(I,J)=.false.\r
2267             if(V.GE.H)then\r
2268               val(i,j)=v\r
2269             else\r
2270               val(i,j)=h\r
2271             end if\r
2272           ENDIF\r
2273         ENDDO\r
2274       ENDDO\r
2275       \r
2276 ***   extract the alignment:\r
2277       i=NSEQ1\r
2278       j=NSEQ2\r
2279       DO WHILE((i.GT.0).AND.(j.GT.0))\r
2280         IF(DIR(i,j))THEN\r
2281           invmap(j)=i\r
2282           i=i-1\r
2283           j=j-1\r
2284         ELSE\r
2285           H=VAL(i-1,j)\r
2286           if(DIR(i-1,j))H=H+GAP_OPEN\r
2287           V=VAL(i,j-1)\r
2288           if(DIR(i,j-1))V=V+GAP_OPEN\r
2289           IF(V.GE.H) THEN\r
2290             j=j-1\r
2291           ELSE\r
2292             i=i-1\r
2293           ENDIF\r
2294         ENDIF\r
2295       ENDDO\r
2296       \r
2297 c^^^^^^^^^^^^^^^Dynamical programming done ^^^^^^^^^^^^^^^^^^^\r
2298       return\r
2299       END\r
2300 \r
2301 cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc\r
2302 c  w    - w(m) is weight for atom pair  c m           (given)\r
2303 c  x    - x(i,m) are coordinates of atom c m in set x       (given)\r
2304 c  y    - y(i,m) are coordinates of atom c m in set y       (given)\r
2305 c  n    - n is number of atom pairs                         (given)\r
2306 c  mode  - 0:calculate rms only                             (given)\r
2307 c          1:calculate rms,u,t                              (takes longer)\r
2308 c  rms   - sum of w*(ux+t-y)**2 over all atom pairs         (result)\r
2309 c  u    - u(i,j) is   rotation  matrix for best superposition  (result)\r
2310 c  t    - t(i)   is translation vector for best superposition  (result)\r
2311 c  ier  - 0: a unique optimal superposition has been determined(result)\r
2312 c       -1: superposition is not unique but optimal\r
2313 c       -2: no result obtained because of negative weights w\r
2314 c           or all weights equal to zero.\r
2315 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
2316       subroutine u3b(w, x, y, n, mode, rms, u, t, ier)\r
2317       integer ip(9), ip2312(4), i, j, k, l, m1, m, ier, n, mode\r
2318       double precision w(1), x(3, 1), y(3, 1), u(3, 3), t(3), rms, sigma\r
2319       double precision r(3, 3), xc(3), yc(3), wc, a(3, 3), b(3, 3), e0, \r
2320      &e(3), e1, e2, e3, d, spur, det, cof, h, g, cth, sth, sqrth, p, tol\r
2321      &, rr(6), rr1, rr2, rr3, rr4, rr5, rr6, ss(6), ss1, ss2, ss3, ss4, \r
2322      &ss5, ss6, zero, one, two, three, sqrt3\r
2323       equivalence (rr(1), rr1), (rr(2), rr2), (rr(3), rr3), (rr(4), rr4)\r
2324      &, (rr(5), rr5), (rr(6), rr6), (ss(1), ss1), (ss(2), ss2), (ss(3), \r
2325      &ss3), (ss(4), ss4), (ss(5), ss5), (ss(6), ss6), (e(1), e1), (e(2)\r
2326      &, e2), (e(3), e3)\r
2327       data sqrt3 / 1.73205080756888d+00 /\r
2328       data tol / 1.0d-2 /\r
2329       data zero / 0.0d+00 /\r
2330       data one / 1.0d+00 /\r
2331       data two / 2.0d+00 /\r
2332       data three / 3.0d+00 /\r
2333       data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /\r
2334       data ip2312 / 2, 3, 1, 2 /\r
2335 c 156 "rms.for"\r
2336       wc = zero\r
2337       rms = 0.0\r
2338       e0 = zero\r
2339       do 1 i = 1, 3\r
2340       xc(i) = zero\r
2341       yc(i) = zero\r
2342       t(i) = 0.0\r
2343       do 1 j = 1, 3\r
2344       d = zero\r
2345       if (i .eq. j) d = one\r
2346       u(i,j) = d\r
2347       a(i,j) = d\r
2348     1 r(i,j) = zero\r
2349       ier = -1\r
2350 c**** DETERMINE CENTROIDS OF BOTH VECTOR SETS X AND Y\r
2351 c 170 "rms.for"\r
2352       if (n .lt. 1) return \r
2353 c 172 "rms.for"\r
2354       ier = -2\r
2355       do 2 m = 1, n\r
2356       if (w(m) .lt. 0.0) return \r
2357       wc = wc + w(m)\r
2358       do 2 i = 1, 3\r
2359       xc(i) = xc(i) + (w(m) * x(i,m))\r
2360     2 yc(i) = yc(i) + (w(m) * y(i,m))\r
2361       if (wc .le. zero) return \r
2362       do 3 i = 1, 3\r
2363       xc(i) = xc(i) / wc\r
2364 c**** DETERMINE CORRELATION MATRIX R BETWEEN VECTOR SETS Y AND X\r
2365 c 182 "rms.for"\r
2366     3 yc(i) = yc(i) / wc\r
2367 c 184 "rms.for"\r
2368       do 4 m = 1, n\r
2369       do 4 i = 1, 3\r
2370       e0 = e0 + (w(m) * (((x(i,m) - xc(i)) ** 2) + ((y(i,m) - yc(i)) ** \r
2371      &2)))\r
2372 c 187 "rms.for"\r
2373       d = w(m) * (y(i,m) - yc(i))\r
2374       do 4 j = 1, 3\r
2375 c**** CALCULATE DETERMINANT OF R(I,J)\r
2376 c 189 "rms.for"\r
2377     4 r(i,j) = r(i,j) + (d * (x(j,m) - xc(j)))\r
2378 c 191 "rms.for"\r
2379       det = ((r(1,1) * ((r(2,2) * r(3,3)) - (r(2,3) * r(3,2)))) - (r(1,2\r
2380      &) * ((r(2,1) * r(3,3)) - (r(2,3) * r(3,1))))) + (r(1,3) * ((r(2,1)\r
2381      & * r(3,2)) - (r(2,2) * r(3,1))))\r
2382 c**** FORM UPPER TRIANGLE OF TRANSPOSED(R)*R\r
2383 c 194 "rms.for"\r
2384       sigma = det\r
2385 c 196 "rms.for"\r
2386       m = 0\r
2387       do 5 j = 1, 3\r
2388       do 5 i = 1, j\r
2389       m = m + 1\r
2390 c***************** EIGENVALUES *****************************************\r
2391 c**** FORM CHARACTERISTIC CUBIC  X**3-3*SPUR*X**2+3*COF*X-DET=0\r
2392 c 200 "rms.for"\r
2393     5 rr(m) = ((r(1,i) * r(1,j)) + (r(2,i) * r(2,j))) + (r(3,i) * r(3,j)\r
2394      &)\r
2395 c 203 "rms.for"\r
2396       spur = ((rr1 + rr3) + rr6) / three\r
2397       cof = ((((((rr3 * rr6) - (rr5 * rr5)) + (rr1 * rr6)) - (rr4 * rr4)\r
2398      &) + (rr1 * rr3)) - (rr2 * rr2)) / three\r
2399 c 205 "rms.for"\r
2400       det = det * det\r
2401       do 6 i = 1, 3\r
2402     6 e(i) = spur\r
2403 c**** REDUCE CUBIC TO STANDARD FORM Y**3-3HY+2G=0 BY PUTTING X=Y+SPUR\r
2404 c 208 "rms.for"\r
2405       if (spur .le. zero) goto 40\r
2406 c 210 "rms.for"\r
2407       d = spur * spur\r
2408       h = d - cof\r
2409 c**** SOLVE CUBIC. ROOTS ARE E1,E2,E3 IN DECREASING ORDER\r
2410 c 212 "rms.for"\r
2411       g = (((spur * cof) - det) / two) - (spur * h)\r
2412 c 214 "rms.for"\r
2413       if (h .le. zero) goto 8\r
2414       sqrth = dsqrt(h)\r
2415       d = ((h * h) * h) - (g * g)\r
2416       if (d .lt. zero) d = zero\r
2417       d = datan2(dsqrt(d),- g) / three\r
2418       cth = sqrth * dcos(d)\r
2419       sth = (sqrth * sqrt3) * dsin(d)\r
2420       e1 = (spur + cth) + cth\r
2421       e2 = (spur - cth) + sth\r
2422       e3 = (spur - cth) - sth\r
2423 c.....HANDLE SPECIAL CASE OF 3 IDENTICAL ROOTS\r
2424 c 224 "rms.for"\r
2425       if (mode) 10, 50, 10\r
2426 c**************** EIGENVECTORS *****************************************\r
2427 c 226 "rms.for"\r
2428     8 if (mode) 30, 50, 30\r
2429 c 228 "rms.for"\r
2430    10 do 15 l = 1, 3, 2\r
2431       d = e(l)\r
2432       ss1 = ((d - rr3) * (d - rr6)) - (rr5 * rr5)\r
2433       ss2 = ((d - rr6) * rr2) + (rr4 * rr5)\r
2434       ss3 = ((d - rr1) * (d - rr6)) - (rr4 * rr4)\r
2435       ss4 = ((d - rr3) * rr4) + (rr2 * rr5)\r
2436       ss5 = ((d - rr1) * rr5) + (rr2 * rr4)\r
2437       ss6 = ((d - rr1) * (d - rr3)) - (rr2 * rr2)\r
2438       j = 1\r
2439       if (dabs(ss1) .ge. dabs(ss3)) goto 12\r
2440       j = 2\r
2441       if (dabs(ss3) .ge. dabs(ss6)) goto 13\r
2442    11 j = 3\r
2443       goto 13\r
2444    12 if (dabs(ss1) .lt. dabs(ss6)) goto 11\r
2445    13 d = zero\r
2446       j = 3 * (j - 1)\r
2447       do 14 i = 1, 3\r
2448       k = ip(i + j)\r
2449       a(i,l) = ss(k)\r
2450    14 d = d + (ss(k) * ss(k))\r
2451       if (d .gt. zero) d = one / dsqrt(d)\r
2452       do 15 i = 1, 3\r
2453    15 a(i,l) = a(i,l) * d\r
2454       d = ((a(1,1) * a(1,3)) + (a(2,1) * a(2,3))) + (a(3,1) * a(3,3))\r
2455       m1 = 3\r
2456       m = 1\r
2457       if ((e1 - e2) .gt. (e2 - e3)) goto 16\r
2458       m1 = 1\r
2459       m = 3\r
2460    16 p = zero\r
2461       do 17 i = 1, 3\r
2462       a(i,m1) = a(i,m1) - (d * a(i,m))\r
2463    17 p = p + (a(i,m1) ** 2)\r
2464       if (p .le. tol) goto 19\r
2465       p = one / dsqrt(p)\r
2466       do 18 i = 1, 3\r
2467    18 a(i,m1) = a(i,m1) * p\r
2468       goto 21\r
2469    19 p = one\r
2470       do 20 i = 1, 3\r
2471       if (p .lt. dabs(a(i,m))) goto 20\r
2472       p = dabs(a(i,m))\r
2473       j = i\r
2474    20 continue\r
2475       k = ip2312(j)\r
2476       l = ip2312(j + 1)\r
2477       p = dsqrt((a(k,m) ** 2) + (a(l,m) ** 2))\r
2478       if (p .le. tol) goto 40\r
2479       a(j,m1) = zero\r
2480       a(k,m1) = - (a(l,m) / p)\r
2481       a(l,m1) = a(k,m) / p\r
2482    21 a(1,2) = (a(2,3) * a(3,1)) - (a(2,1) * a(3,3))\r
2483       a(2,2) = (a(3,3) * a(1,1)) - (a(3,1) * a(1,3))\r
2484 c****************** ROTATION MATRIX ************************************\r
2485 c 282 "rms.for"\r
2486       a(3,2) = (a(1,3) * a(2,1)) - (a(1,1) * a(2,3))\r
2487 c 284 "rms.for"\r
2488    30 do 32 l = 1, 2\r
2489       d = zero\r
2490       do 31 i = 1, 3\r
2491       b(i,l) = ((r(i,1) * a(1,l)) + (r(i,2) * a(2,l))) + (r(i,3) * a(3,l\r
2492      &))\r
2493 c 288 "rms.for"\r
2494    31 d = d + (b(i,l) ** 2)\r
2495       if (d .gt. zero) d = one / dsqrt(d)\r
2496       do 32 i = 1, 3\r
2497    32 b(i,l) = b(i,l) * d\r
2498       d = ((b(1,1) * b(1,2)) + (b(2,1) * b(2,2))) + (b(3,1) * b(3,2))\r
2499       p = zero\r
2500       do 33 i = 1, 3\r
2501       b(i,2) = b(i,2) - (d * b(i,1))\r
2502    33 p = p + (b(i,2) ** 2)\r
2503       if (p .le. tol) goto 35\r
2504       p = one / dsqrt(p)\r
2505       do 34 i = 1, 3\r
2506    34 b(i,2) = b(i,2) * p\r
2507       goto 37\r
2508    35 p = one\r
2509       do 36 i = 1, 3\r
2510       if (p .lt. dabs(b(i,1))) goto 36\r
2511       p = dabs(b(i,1))\r
2512       j = i\r
2513    36 continue\r
2514       k = ip2312(j)\r
2515       l = ip2312(j + 1)\r
2516       p = dsqrt((b(k,1) ** 2) + (b(l,1) ** 2))\r
2517       if (p .le. tol) goto 40\r
2518       b(j,2) = zero\r
2519       b(k,2) = - (b(l,1) / p)\r
2520       b(l,2) = b(k,1) / p\r
2521    37 b(1,3) = (b(2,1) * b(3,2)) - (b(2,2) * b(3,1))\r
2522       b(2,3) = (b(3,1) * b(1,2)) - (b(3,2) * b(1,1))\r
2523       b(3,3) = (b(1,1) * b(2,2)) - (b(1,2) * b(2,1))\r
2524       do 39 i = 1, 3\r
2525       do 39 j = 1, 3\r
2526 c****************** TRANSLATION VECTOR *********************************\r
2527 c 320 "rms.for"\r
2528    39 u(i,j) = ((b(i,1) * a(j,1)) + (b(i,2) * a(j,2))) + (b(i,3) * a(j,3\r
2529      &))\r
2530    40 do 41 i = 1, 3\r
2531 c********************** RMS ERROR **************************************\r
2532 c 323 "rms.for"\r
2533    41 t(i) = ((yc(i) - (u(i,1) * xc(1))) - (u(i,2) * xc(2))) - (u(i,3)\r
2534      & * xc(3))\r
2535    50 do 51 i = 1, 3\r
2536       if (e(i) .lt. zero) e(i) = zero\r
2537    51 e(i) = dsqrt(e(i))\r
2538       ier = 0\r
2539       if (e2 .le. (e1 * 1.0d-05)) ier = -1\r
2540       d = e3\r
2541       if (sigma .ge. 0.0) goto 52\r
2542       d = - d\r
2543       if ((e2 - e3) .le. (e1 * 1.0d-05)) ier = -1\r
2544    52 d = (d + e2) + e1\r
2545       rms = (e0 - d) - d\r
2546       if (rms .lt. 0.0) rms = 0.0\r
2547       return \r
2548 c.....END U3B...........................................................\r
2549 c----------------------------------------------------------\r
2550 c                       THE END\r
2551 c----------------------------------------------------------\r
2552 c 338 "rms.for"\r
2553       end\r
2554 \r