+++ /dev/null
-*************************************************************************\r
-* This program is to identify the best alignment of two protein \r
-* structures to give the best TM-score. By default, TM-score is \r
-* normalized by the second protein. The program can be freely \r
-* copied or modified or redistributed.\r
-* (For comments, please email to: yzhang@ku.edu)\r
-* \r
-* Reference:\r
-* Yang Zhang, Jeffrey Skolnick, Nucl. Acid Res. 2005 33: 2303-9\r
-*\r
-************************ updating history *******************************\r
-* 2005/06/01: A small bug of two-point superposition was fixed.\r
-* 2005/10/19: the program was reformed so that the alignment\r
-* results are not dependent on the specific compilers.\r
-* 2006/06/20: select 'A' if there is altLoc when reading PDB file.\r
-* 2007/02/27: rotation matrix from Chain-1 to Chain-2 is added.\r
-* 2007/04/18: add options with TM-score normalized by average\r
-* length, shorter length, or longer length of two \r
-* structures.\r
-* 2007/05/23: add additional output file 'TM.sup_all' for showing\r
-* all atoms while 'TM.sup' is only for aligned atoms\r
-* 2007/09/19: add a new feature alignment to deal with the problem\r
-* of aligning fractional structures (e.g. protein\r
-* interfaces).\r
-*************************************************************************\r
- \r
- program compares\r
- PARAMETER(nmax=5000)\r
- PARAMETER(nmax2=10000)\r
-\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/alignrst/invmap0(nmax)\r
- common/length/nseq1,nseq2\r
- common/d0/d0,anseq\r
- common/d0min/d0_min\r
- common/d00/d00,d002\r
-\r
- character*100 fnam,pdb(100),outname\r
- character*3 aa(-1:20),aanam,ss1(nmax),ss2(nmax)\r
- character*100 s,du\r
- character*200 outnameall_tmp,outnameall\r
- character seq1(0:nmax),seq2(0:nmax)\r
- character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2)\r
-\r
- dimension xx(nmax),yy(nmax),zz(nmax)\r
- dimension m1(nmax),m2(nmax)\r
- dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
- dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
- common/init/invmap_i(nmax)\r
-\r
- common/TM/TM,TMmax\r
- common/n1n2/n1(nmax),n2(nmax)\r
- common/d8/d8\r
- common/initial4/mm1(nmax),mm2(nmax)\r
-\r
-ccc RMSD:\r
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
- double precision u(3,3),t(3),rms,drms !armsd is real\r
- data w /nmax*1.0/\r
-ccc \r
-\r
- data aa/ 'BCK','GLY','ALA','SER','CYS','VAL','THR','ILE',\r
- & 'PRO','MET','ASP','ASN','LEU',\r
- & 'LYS','GLU','GLN','ARG',\r
- & 'HIS','PHE','TYR','TRP','CYX'/\r
- character*1 slc(-1:20)\r
- data slc/'X','G','A','S','C','V','T','I',\r
- & 'P','M','D','N','L','K','E','Q','R',\r
- & 'H','F','Y','W','C'/\r
-\r
- call getarg(1,fnam)\r
- if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then\r
- write(*,*)\r
- write(*,*)'Brief instruction for running TM-align program:'\r
- write(*,*)'(For detail: Zhang & Skolnick, Nucl. Acid Res.',\r
- & '2005 33, 2303)'\r
- write(*,*)\r
- write(*,*)'1. Align ''structure.pdb'' to ''target.pdb'''\r
- write(*,*)' (By default, TM-score is normalized by the ',\r
- & 'length of ''target.pdb'')'\r
- write(*,*)' >TMalign structure.pdb target.pdb'\r
- write(*,*)\r
- write(*,*)'2. Run TM-align and output the superposition ',\r
- & 'to ''TM.sup'' and ''TM.sup_all'':'\r
- write(*,*)' >TMalign structure.pdb target.pdb -o TM.sup'\r
- write(*,*)' To view the superimposed structures of the',\r
- & ' aligned regions by rasmol:'\r
- write(*,*)' >rasmol -script TM.sup)'\r
- write(*,*)' To view the superimposed structures of all',\r
- & ' regions by rasmol:'\r
- write(*,*)' >rasmol -script TM.sup_all)'\r
- write(*,*)\r
- write(*,*)'3. If you want TM-score normalized by ',\r
- & 'an assigned length, e.g. 100 aa:'\r
- write(*,*)' >TMalign structure.pdb target.pdb -L 100'\r
- write(*,*)' If you want TM-score normalized by the ',\r
- & 'average length of two structures:'\r
- write(*,*)' >TMalign structure.pdb target.pdb -a'\r
- write(*,*)' If you want TM-score normalized by the ',\r
- & 'shorter length of two structures:'\r
- write(*,*)' >TMalign structure.pdb target.pdb -b'\r
- write(*,*)' If you want TM-score normalized by the ',\r
- & 'longer length of two structures:'\r
- write(*,*)' >TMalign structure.pdb target.pdb -c'\r
- write(*,*)\r
-c write(*,*)'5. If you want to set a minimum cutoff for d0, ',\r
-c & 'e.g. d0>3.0'\r
-c write(*,*)' (By default d0>0.5, this option need ',\r
-c & 'be considered only when L<35 aa)'\r
-c write(*,*)' >TMalign structure.pdb target.pdb -dmin 3.0'\r
-c write(*,*)\r
- write(*,*)'(All above options does not change the ',\r
- & 'final structure alignment result)'\r
- write(*,*)\r
- goto 9999\r
- endif\r
-\r
-******* options ----------->\r
- m_out=-1 !decided output\r
- m_fix=-1 !fixed length-scale only for output\r
- m_ave=-1 !using average length\r
- m_d0_min=-1 !diminum d0 for search\r
- m_d0=-1 !given d0 for both search and output\r
- narg=iargc()\r
- i=0\r
- j=0\r
- 115 continue\r
- i=i+1\r
- call getarg(i,fnam)\r
- if(fnam.eq.'-o')then\r
- m_out=1\r
- i=i+1\r
- call getarg(i,outname)\r
- elseif(fnam.eq.'-L')then !change both L_all and d0\r
- m_fix=1\r
- i=i+1\r
- call getarg(i,fnam)\r
- read(fnam,*)L_fix\r
- elseif(fnam.eq.'-dmin')then\r
- m_d0_min=1\r
- i=i+1\r
- call getarg(i,fnam)\r
- read(fnam,*)d0_min_input\r
- elseif(fnam.eq.'-d0')then\r
- m_d0=1\r
- i=i+1\r
- call getarg(i,fnam)\r
- read(fnam,*)d0_fix\r
- elseif(fnam.eq.'-a')then ! this will change the superposed output but not the alignment\r
- m_ave=1\r
- i=i+1\r
- elseif(fnam.eq.'-b')then\r
- m_ave=2\r
- i=i+1\r
- elseif(fnam.eq.'-c')then\r
- m_ave=3\r
- i=i+1\r
- else\r
- j=j+1\r
- pdb(j)=fnam\r
- endif\r
- if(i.lt.narg)goto 115\r
- \r
-ccccccccc read data from first CA file:\r
- open(unit=10,file=pdb(1),status='old')\r
- i=0\r
- do while (.true.)\r
- read(10,9001,end=1010) s\r
- if(i.gt.0.and.s(1:3).eq.'TER')goto 1010\r
- if(s(1:3).eq.'ATO')then\r
- if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or\r
- & .s(13:16).eq.' CA')then\r
- if(s(17:17).eq.' '.or.s(17:17).eq.'A')then\r
- i=i+1\r
- read(s,9000)du,aanam,du,mm1(i),du,\r
- $ xa(1,i,0),xa(2,i,0),xa(3,i,0)\r
- do j=-1,20\r
- if(aanam.eq.aa(j))seq1(i)=slc(j)\r
- enddo\r
- ss1(i)=aanam\r
- if(i.ge.nmax)goto 1010\r
- endif\r
- endif\r
- endif\r
- enddo\r
- 1010 continue\r
- 9000 format(A17,A3,A2,i4,A4,3F8.3)\r
- 9001 format(A100)\r
- close(10)\r
- nseq1=i\r
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
-\r
-ccccccccc read data from the second CA file:\r
- open(unit=10,file=pdb(2),status='old')\r
- i=0\r
- do while (.true.)\r
- read(10,9001,end=1011) s\r
- if(i.gt.0.and.s(1:3).eq.'TER')goto 1011\r
- if(s(1:3).eq.'ATO')then\r
- if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.\r
- & s(13:16).eq.' CA')then\r
- if(s(17:17).eq.' '.or.s(17:17).eq.'A')then\r
- i=i+1\r
- read(s,9000)du,aanam,du,mm2(i),du,\r
- $ xa(1,i,1),xa(2,i,1),xa(3,i,1)\r
- do j=-1,20\r
- if(aanam.eq.aa(j))seq2(i)=slc(j)\r
- enddo\r
- ss2(i)=aanam\r
- if(i.ge.nmax)goto 1011\r
- endif\r
- endif\r
- endif\r
- enddo\r
- 1011 continue\r
- close(10)\r
- nseq2=i\r
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
-\r
-*!!! Scale of TM-score in search is based on the smaller protein --------->\r
- d0_min=0.5\r
- if(m_d0_min.eq.1)then\r
- d0_min=d0_min_input !for search\r
- endif\r
- anseq_min=min(nseq1,nseq2)\r
- anseq=anseq_min !length for defining TMscore in search\r
- d8=1.5*anseq_min**0.3+3.5 !remove pairs with dis>d8 during search & final\r
- if(anseq.gt.15)then\r
- d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score\r
- else\r
- d0=d0_min\r
- endif\r
- if(d0.lt.d0_min)d0=d0_min\r
- if(m_d0.eq.1)d0=d0_fix\r
- d00=d0 !for quickly calculate TM-score in searching\r
- if(d00.gt.8)d00=8\r
- if(d00.lt.4.5)d00=4.5\r
- d002=d00**2\r
- nseq=max(nseq1,nseq2)\r
- do i=1,nseq\r
- n1(i)=i\r
- n2(i)=i\r
- enddo\r
- \r
-***** do alignment **************************\r
- CALL super_align !to find invmap(j)\r
- \r
-************************************************************\r
-*** resuperpose to find residues of dis<d8 ------------------------>\r
- n_al=0\r
- do j=1,nseq2\r
- if(invmap0(j).gt.0)then\r
- i=invmap0(j)\r
- n_al=n_al+1\r
- xtm1(n_al)=xa(1,i,0)\r
- ytm1(n_al)=xa(2,i,0)\r
- ztm1(n_al)=xa(3,i,0)\r
- xtm2(n_al)=xa(1,j,1)\r
- ytm2(n_al)=xa(2,j,1)\r
- ztm2(n_al)=xa(3,j,1)\r
- m1(n_al)=i !for recording residue order\r
- m2(n_al)=j\r
- endif\r
- enddo\r
- d0_input=d0\r
- call TMscore8(d0_input,n_al,xtm1,ytm1,ztm1,n1,n_al,\r
- & xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !TM-score with dis<d8 only\r
-\r
-*!!! Output TM-score is based on the second protein------------------>\r
- d0_min=0.5 !for output\r
- anseq=nseq2 !length for defining final TMscore\r
- if(m_ave.eq.1)anseq=(nseq1+nseq2)/2.0 !<L>\r
- if(m_ave.eq.2)anseq=min(nseq1,nseq2)\r
- if(m_ave.eq.3)anseq=max(nseq1,nseq2)\r
- if(anseq.lt.anseq_min)anseq=anseq_min\r
- if(m_fix.eq.1)anseq=L_fix !input length\r
- if(anseq.gt.15)then\r
- d0=1.24*(anseq-15)**(1.0/3.0)-1.8 !scale for defining TM-score\r
- else\r
- d0=d0_min\r
- endif\r
- if(d0.lt.d0_min)d0=d0_min\r
- if(m_d0.eq.1)d0=d0_fix\r
- \r
-*** remove dis>d8 in normal TM-score calculation for final report----->\r
- j=0\r
- n_eq=0\r
- do i=1,n_al\r
- dis2=sqrt((xtm1(i)-xtm2(i))**2+(ytm1(i)-ytm2(i))**2+\r
- & (ztm1(i)-ztm2(i))**2)\r
- if(dis2.le.d8)then\r
- j=j+1\r
- xtm1(j)=xtm1(i)\r
- ytm1(j)=ytm1(i)\r
- ztm1(j)=ztm1(i)\r
- xtm2(j)=xtm2(i)\r
- ytm2(j)=ytm2(i)\r
- ztm2(j)=ztm2(i)\r
- m1(j)=m1(i)\r
- m2(j)=m2(i)\r
- if(ss1(m1(i)).eq.ss2(m2(i)))then\r
- n_eq=n_eq+1\r
- endif\r
- endif\r
- enddo\r
- seq_id=float(n_eq)/(n_al+0.00000001)\r
- n8_al=j\r
- d0_input=d0\r
- call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n1,n8_al,\r
- & xtm2,ytm2,ztm2,n2,TM8,Rcomm,Lcomm) !normal TMscore\r
- rmsd8_al=Rcomm\r
- TM8=TM8*n8_al/anseq !TM-score after cutoff\r
- \r
-********* for output summary ******************************\r
- write(*,*)\r
- write(*,*)'*****************************************************',\r
- & '*********************'\r
- write(*,*)'* TM-align ',\r
- & ' *'\r
- write(*,*)'* A protein structural alignment algorithm based on T',\r
- & 'M-score *'\r
- write(*,*)'* Reference: Y. Zhang and J. Skolnick, Nucl. Acids Re',\r
- & 's. 2005 33, 2302-9 *'\r
- write(*,*)'* Comments on the program, please email to: yzhang@ku',\r
- & '.edu *'\r
- write(*,*)'*****************************************************',\r
- & '*********************'\r
- write(*,*)\r
- write(*,101)pdb(1),nseq1\r
- 101 format('Chain 1:',A10,' Size=',I4)\r
- write(*,102)pdb(2),nseq2,int(anseq)\r
- 102 format('Chain 2:',A10,' Size=',I4,\r
- & ' (TM-score is normalized by ',I4,')')\r
- write(*,*)\r
- write(*,103)n8_al,rmsd8_al,TM8,seq_id\r
- 103 format('Aligned length=',I4,', RMSD=',f6.2,\r
- & ', TM-score=',f7.5,', ID=',f5.3)\r
- write(*,*)\r
-\r
-********* extract rotation matrix ------------>\r
- L=0\r
- do i=1,n8_al\r
- k=m1(i)\r
- L=L+1\r
- r_1(1,L)=xa(1,k,0)\r
- r_1(2,L)=xa(2,k,0)\r
- r_1(3,L)=xa(3,k,0)\r
- r_2(1,L)=xtm1(i)\r
- r_2(2,L)=ytm1(i)\r
- r_2(3,L)=ztm1(i)\r
- enddo\r
- if(L.gt.3)then\r
- call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- armsd=dsqrt(rms/L)\r
- write(*,*)'-------- rotation matrix to rotate Chain-1 to ',\r
- & 'Chain-2 ------'\r
- write(*,*)'i t(i) u(i,1) u(i,2) ',\r
- & ' u(i,3)'\r
- do i=1,3\r
- write(*,204)i,t(i),u(i,1),u(i,2),u(i,3)\r
- enddo\r
-c do i=1,nseq1\r
-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
-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
-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
-c enddo\r
- write(*,*)\r
- endif\r
- 204 format(I2,f18.10,f15.10,f15.10,f15.10)\r
- \r
-********* for output superposition ******************************\r
- if(m_out.eq.1)then\r
- 1237 format('ATOM ',i5,' CA ',A3,I6,4X,3F8.3)\r
- 1238 format('TER')\r
- 1239 format('CONECT',I5,I5)\r
- 900 format(A)\r
- 901 format('select atomno=',I4)\r
- 104 format('REMARK Chain 1:',A10,' Size=',I4)\r
- 105 format('REMARK Chain 2:',A10,' Size=',I4,\r
- & ' (TM-score is normalized by ',I4,')')\r
- 106 format('REMARK Aligned length=',I4,', RMSD=',f6.2,\r
- & ', TM-score=',f7.5,', ID=',f5.3)\r
- OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln\r
-*** script:\r
- write(7,900)'load inline'\r
- write(7,900)'select atomno<2000'\r
- write(7,900)'wireframe .45'\r
- write(7,900)'select none'\r
- write(7,900)'select atomno>2000'\r
- write(7,900)'wireframe .20'\r
- write(7,900)'color white'\r
- do i=1,n8_al\r
- dis2=sqrt((xtm1(i)-xtm2(i))**2+\r
- & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)\r
- if(dis2.le.5)then\r
- write(7,901)m1(i)\r
- write(7,900)'color red'\r
- write(7,901)2000+m2(i)\r
- write(7,900)'color red'\r
- endif\r
- enddo\r
- write(7,900)'select all'\r
- write(7,900)'exit'\r
- write(7,104)pdb(1),nseq1\r
- write(7,105)pdb(2),nseq2,int(anseq)\r
- write(7,106)n8_al,rmsd8_al,TM8,seq_id\r
-*** chain1:\r
- do i=1,n8_al\r
- write(7,1237)m1(i),ss1(m1(i)),mm1(m1(i)),\r
- & xtm1(i),ytm1(i),ztm1(i)\r
- enddo\r
- write(7,1238) !TER\r
- do i=2,n8_al\r
- write(7,1239)m1(i-1),m1(i) !connect atoms\r
- enddo\r
-*** chain2:\r
- do i=1,n8_al\r
- write(7,1237)2000+m2(i),ss2(m2(i)),mm2(m2(i)),\r
- $ xtm2(i),ytm2(i),ztm2(i)\r
- enddo\r
- write(7,1238)\r
- do i=2,n8_al\r
- write(7,1239)2000+m2(i-1),2000+m2(i)\r
- enddo\r
- close(7)\r
-ccc \r
- k=0\r
- outnameall_tmp=outname//'_all'\r
- outnameall=''\r
- do i=1,200\r
- if(outnameall_tmp(i:i).ne.' ')then\r
- k=k+1\r
- outnameall(k:k)=outnameall_tmp(i:i)\r
- endif\r
- enddo\r
- OPEN(unit=8,file=outnameall,status='unknown') !pdb1.aln + pdb2.aln\r
-*** script:\r
- write(8,900)'load inline'\r
- write(8,900)'select atomno<2000'\r
- write(8,900)'wireframe .45'\r
- write(8,900)'select none'\r
- write(8,900)'select atomno>2000'\r
- write(8,900)'wireframe .20'\r
- write(8,900)'color white'\r
- do i=1,n8_al\r
- dis2=sqrt((xtm1(i)-xtm2(i))**2+\r
- & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)\r
- if(dis2.le.5)then\r
- write(8,901)m1(i)\r
- write(8,900)'color red'\r
- write(8,901)2000+m2(i)\r
- write(8,900)'color red'\r
- endif\r
- enddo\r
- write(8,900)'select all'\r
- write(8,900)'exit'\r
- write(8,104)pdb(1),nseq1\r
- write(8,105)pdb(2),nseq2,int(anseq)\r
- write(8,106)n8_al,rmsd8_al,TM8,seq_id\r
-*** chain1:\r
- do i=1,nseq1\r
- 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
- 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
- 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
- write(8,1237)i,ss1(i),mm1(i),ax,ay,az\r
- enddo\r
- write(8,1238) !TER\r
- do i=2,nseq1\r
- write(8,1239)i-1,i\r
- enddo\r
-*** chain2:\r
- do i=1,nseq2\r
- write(8,1237)2000+i,ss2(i),mm2(i),\r
- $ xa(1,i,1),xa(2,i,1),xa(3,i,1)\r
- enddo\r
- write(8,1238)\r
- do i=2,nseq2\r
- write(8,1239)2000+i-1,2000+i\r
- enddo\r
- close(8)\r
- endif\r
-*^^^^^^^^^^^^^^^^^^ output finished ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
-\r
-************ output aligned sequences **************************\r
- ii=0\r
- i1_old=1\r
- i2_old=1\r
- do i=1,n8_al\r
- do j=i1_old,m1(i)-1\r
- ii=ii+1\r
- aseq1(ii)=seq1(j)\r
- aseq2(ii)='-'\r
- aseq3(ii)=' '\r
- enddo\r
- do j=i2_old,m2(i)-1\r
- ii=ii+1\r
- aseq1(ii)='-'\r
- aseq2(ii)=seq2(j)\r
- aseq3(ii)=' '\r
- enddo\r
- ii=ii+1\r
- aseq1(ii)=seq1(m1(i))\r
- aseq2(ii)=seq2(m2(i))\r
- dis2=sqrt((xtm1(i)-xtm2(i))**2+\r
- & (ytm1(i)-ytm2(i))**2+(ztm1(i)-ztm2(i))**2)\r
- if(dis2.le.5)then\r
- aseq3(ii)=':'\r
- else\r
- aseq3(ii)='.'\r
- endif\r
- i1_old=m1(i)+1\r
- i2_old=m2(i)+1\r
- enddo\r
- do i=i1_old,nseq1\r
- ii=ii+1\r
- aseq1(ii)=seq1(i)\r
- aseq2(ii)='-'\r
- aseq3(ii)=' '\r
- enddo\r
- do i=i2_old,nseq2\r
- ii=ii+1\r
- aseq1(ii)='-'\r
- aseq2(ii)=seq2(i)\r
- aseq3(ii)=' '\r
- enddo\r
- write(*,50)\r
- 50 format('(":" denotes the residue pairs of distance < 5.0 ',\r
- & 'Angstrom)')\r
- write(*,10)(aseq1(i),i=1,ii)\r
- write(*,10)(aseq3(i),i=1,ii)\r
- write(*,10)(aseq2(i),i=1,ii)\r
- 10 format(10000A1)\r
- write(*,*)\r
-\r
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- 9999 END\r
-\r
-***********************************************************************\r
-***********************************************************************\r
-* Structure superposition\r
-***********************************************************************\r
-***********************************************************************\r
-***********************************************************************\r
- SUBROUTINE super_align\r
- PARAMETER(nmax=5000)\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/length/nseq1,nseq2\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/alignrst/invmap0(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/TM/TM,TMmax\r
- common/init/invmap_i(nmax)\r
- dimension gapp(100)\r
-\r
- TMmax=0\r
- n_gapp=2\r
- gapp(1)=-0.6\r
- gapp(2)=0\r
-\r
-c n_gapp=11\r
-c do i=1,n_gapp\r
-c gapp(i)=-(n_gapp-i)\r
-c enddo\r
-\r
-*11111111111111111111111111111111111111111111111111111111\r
-* get initial alignment from gapless threading\r
-**********************************************************\r
- call get_initial !gapless threading\r
- do i=1,nseq2\r
- invmap(i)=invmap_i(i) !with highest zcore\r
- enddo\r
- call get_score !TM, matrix score(i,j)\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
-\r
-*****************************************************************\r
-* initerative alignment, for different gap_open:\r
-*****************************************************************\r
- DO 111 i_gapp=1,n_gapp !different gap panalties\r
- GAP_OPEN=gapp(i_gapp) !gap panalty\r
- do 222 id=1,30 !maximum interation is 200\r
- call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
-* Input: score(i,j), and gap_open\r
-* Output: invmap(j)\r
-\r
- call get_score !calculate TM-score, score(i,j)\r
-c record the best alignment in whole search ---------->\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
- if(id.gt.1)then\r
- diff=abs(TM-TM_old)\r
- if(diff.lt.0.000001)goto 33\r
- endif\r
- TM_old=TM\r
- 222 continue\r
- 33 continue\r
- 111 continue\r
-\r
-*222222222222222222222222222222222222222222222222222222222\r
-* get initial alignment from secondary structure alignment\r
-**********************************************************\r
- call get_initial2 !DP for secondary structure\r
- do i=1,nseq2\r
- invmap(i)=invmap_i(i) !with highest zcore\r
- enddo\r
- call get_score !TM, score(i,j)\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
-\r
-*****************************************************************\r
-* initerative alignment, for different gap_open:\r
-*****************************************************************\r
- DO 1111 i_gapp=1,n_gapp !different gap panalties\r
- GAP_OPEN=gapp(i_gapp) !gap panalty\r
- do 2222 id=1,30 !maximum interation is 200\r
- call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
-* Input: score(i,j), and gap_open\r
-* Output: invmap(j)\r
-\r
- call get_score !calculate TM-score, score(i,j)\r
-c write(*,21)gap_open,rmsd_al,n_al,TM\r
-c record the best alignment in whole search ---------->\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
- if(id.gt.1)then\r
- diff=abs(TM-TM_old)\r
- if(diff.lt.0.000001)goto 333\r
- endif\r
- TM_old=TM\r
- 2222 continue\r
- 333 continue\r
- 1111 continue\r
- \r
-*333333333333333333333333333333333333333333333333333333333333\r
-* get initial alignment from invmap0+SS\r
-*************************************************************\r
- call get_initial3 !invmap0+SS\r
- do i=1,nseq2\r
- invmap(i)=invmap_i(i) !with highest zcore\r
- enddo\r
- call get_score !TM, score(i,j)\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
-\r
-*****************************************************************\r
-* initerative alignment, for different gap_open:\r
-*****************************************************************\r
- DO 1110 i_gapp=1,n_gapp !different gap panalties\r
- GAP_OPEN=gapp(i_gapp) !gap panalty\r
- do 2220 id=1,30 !maximum interation is 200\r
- call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
-* Input: score(i,j), and gap_open\r
-* Output: invmap(j)\r
-\r
- call get_score !calculate TM-score, score(i,j)\r
-c write(*,21)gap_open,rmsd_al,n_al,TM\r
-c record the best alignment in whole search ---------->\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
- if(id.gt.1)then\r
- diff=abs(TM-TM_old)\r
- if(diff.lt.0.000001)goto 330\r
- endif\r
- TM_old=TM\r
- 2220 continue\r
- 330 continue\r
- 1110 continue\r
-\r
-*444444444444444444444444444444444444444444444444444444444\r
-* get initial alignment of pieces from gapless threading\r
-**********************************************************\r
- call get_initial4 !gapless threading\r
- do i=1,nseq2\r
- invmap(i)=invmap_i(i) !with highest zcore\r
- enddo\r
- call get_score !TM, matrix score(i,j)\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
-\r
-*****************************************************************\r
-* initerative alignment, for different gap_open:\r
-*****************************************************************\r
- DO 44 i_gapp=2,n_gapp !different gap panalties\r
- GAP_OPEN=gapp(i_gapp) !gap panalty\r
- do 444 id=1,2 !maximum interation is 200\r
- call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
-* Input: score(i,j), and gap_open\r
-* Output: invmap(j)\r
- \r
- call get_score !calculate TM-score, score(i,j)\r
-c record the best alignment in whole search ---------->\r
- if(TM.gt.TMmax)then\r
- TMmax=TM\r
- do j=1,nseq2\r
- invmap0(j)=invmap(j)\r
- enddo\r
- endif\r
- 444 continue\r
- 44 continue\r
-\r
-c^^^^^^^^^^^^^^^ best alignment invmap0(j) found ^^^^^^^^^^^^^^^^^^\r
- RETURN\r
- END\r
-\r
-**************************************************************\r
-* get initial alignment invmap0(i) from gapless threading\r
-**************************************************************\r
- subroutine get_initial\r
- PARAMETER(nmax=5000)\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/length/nseq1,nseq2\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/alignrst/invmap0(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/TM/TM,TMmax\r
- common/init/invmap_i(nmax)\r
-\r
- aL=min(nseq1,nseq2)\r
- idel=aL/2.5 !minimum size of considered fragment\r
- if(idel.le.5)idel=5\r
- n1=-nseq2+idel\r
- n2=nseq1-idel\r
- GL_max=0\r
- do ishift=n1,n2\r
- L=0\r
- do j=1,nseq2\r
- i=j+ishift\r
- if(i.ge.1.and.i.le.nseq1)then\r
- L=L+1\r
- invmap(j)=i\r
- else\r
- invmap(j)=-1\r
- endif\r
- enddo\r
- if(L.ge.idel)then\r
- call get_GL(GL)\r
- if(GL.gt.GL_max)then\r
- GL_max=GL\r
- do i=1,nseq2\r
- invmap_i(i)=invmap(i)\r
- enddo\r
- endif\r
- endif\r
- enddo\r
-\r
- return\r
- end\r
-\r
-**************************************************************\r
-* get initial alignment invmap0(i) from secondary structure\r
-**************************************************************\r
- subroutine get_initial2\r
- PARAMETER(nmax=5000)\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/length/nseq1,nseq2\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/alignrst/invmap0(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/TM/TM,TMmax\r
- common/sec/isec(nmax),jsec(nmax)\r
- common/init/invmap_i(nmax)\r
-\r
-********** assign secondary structures ***************\r
-c 1->coil, 2->helix, 3->turn, 4->strand\r
- do i=1,nseq1\r
- isec(i)=1\r
- j1=i-2\r
- j2=i-1\r
- j3=i\r
- j4=i+1\r
- j5=i+2\r
- if(j1.ge.1.and.j5.le.nseq1)then\r
- dis13=diszy(0,j1,j3)\r
- dis14=diszy(0,j1,j4)\r
- dis15=diszy(0,j1,j5)\r
- dis24=diszy(0,j2,j4)\r
- dis25=diszy(0,j2,j5)\r
- dis35=diszy(0,j3,j5)\r
- isec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35)\r
- endif\r
- enddo\r
- do i=1,nseq2\r
- jsec(i)=1\r
- j1=i-2\r
- j2=i-1\r
- j3=i\r
- j4=i+1\r
- j5=i+2\r
- if(j1.ge.1.and.j5.le.nseq2)then\r
- dis13=diszy(1,j1,j3)\r
- dis14=diszy(1,j1,j4)\r
- dis15=diszy(1,j1,j5)\r
- dis24=diszy(1,j2,j4)\r
- dis25=diszy(1,j2,j5)\r
- dis35=diszy(1,j3,j5)\r
- jsec(i)=make_sec(dis13,dis14,dis15,dis24,dis25,dis35)\r
- endif\r
- enddo\r
- call smooth !smooth the assignment\r
-\r
-********** score matrix **************************\r
- do i=1,nseq1\r
- do j=1,nseq2\r
- if(isec(i).eq.jsec(j))then\r
- score(i,j)=1\r
- else\r
- score(i,j)=0\r
- endif\r
- enddo\r
- enddo\r
-\r
-********** find initial alignment: invmap(j) ************\r
- gap_open=-1.0 !should be -1\r
- call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
- do i=1,nseq2\r
- invmap_i(i)=invmap(i)\r
- enddo\r
-\r
-*^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- end\r
-\r
-**************************************************************\r
-* get initial alignment invmap0(i) from secondary structure \r
-* and previous alignments\r
-**************************************************************\r
- subroutine get_initial3\r
- PARAMETER(nmax=5000)\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/length/nseq1,nseq2\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/alignrst/invmap0(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/TM/TM,TMmax\r
- common/sec/isec(nmax),jsec(nmax)\r
- common/init/invmap_i(nmax)\r
-\r
-********** score matrix **************************\r
- do i=1,nseq2\r
- invmap(i)=invmap0(i)\r
- enddo\r
- call get_score1 !get score(i,j) using RMSD martix\r
- do i=1,nseq1\r
- do j=1,nseq2\r
- if(isec(i).eq.jsec(j))then\r
- score(i,j)=0.5+score(i,j)\r
- else\r
- score(i,j)=score(i,j)\r
- endif\r
- enddo\r
- enddo\r
-\r
-********** find initial alignment: invmap(j) ************\r
- gap_open=-1.0 !should be -1\r
- call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)\r
- do i=1,nseq2\r
- invmap_i(i)=invmap(i)\r
- enddo\r
-\r
-*^^^^^^^^^^^^ initial alignment done ^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- end\r
-\r
-**************************************************************\r
-* get initial alignment invmap0(i) from fragment gapless threading\r
-**************************************************************\r
- subroutine get_initial4\r
- PARAMETER(nmax=5000)\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/length/nseq1,nseq2\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/alignrst/invmap0(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/TM/TM,TMmax\r
- common/init/invmap_i(nmax)\r
- common/initial4/mm1(nmax),mm2(nmax)\r
- logical contin\r
-\r
- dimension ifr2(2,nmax,nmax),Lfr2(2,nmax),Lfr_max2(2),i_fr2(2)\r
- dimension ifr(nmax)\r
- dimension mm(2,nmax)\r
- \r
- fra_min=4 !>=4,minimum fragment for search\r
- fra_min1=fra_min-1 !cutoff for shift, save time\r
- dcu0=3.85\r
- dcu_min=3.65\r
-\r
-ccc Find the smallest continuous fragments -------->\r
- do i=1,nseq1\r
- mm(1,i)=mm1(i)\r
- enddo\r
- do i=1,nseq2\r
- mm(2,i)=mm2(i)\r
- enddo\r
- do k=1,2\r
- dcu=dcu0\r
- if(k.eq.1)then\r
- nseq0=nseq1\r
- r_min=nseq1/3.0 !minimum fragment, in case too small protein\r
- else\r
- nseq0=nseq2\r
- r_min=nseq2/3.0 !minimum fragment, in case too small protein\r
- endif\r
- if(r_min.gt.fra_min)r_min=fra_min\r
- 20 nfr=1 !number of fragments\r
- j=1 !number of residue at nf-fragment\r
- ifr2(k,nfr,j)=1 !what residue\r
- Lfr2(k,nfr)=j !length of the fragment\r
- do i=2,nseq0\r
- dis=diszy(k-1,i-1,i)\r
- contin=.false.\r
- if(dcu.gt.dcu0)then\r
- if(dis.lt.dcu)then\r
- if(dis.gt.dcu_min)then\r
- contin=.true.\r
- endif\r
- endif\r
- elseif(mm(k,i).eq.(mm(k,i-1)+1))then\r
- if(dis.lt.dcu)then\r
- if(dis.gt.dcu_min)then\r
- contin=.true.\r
- endif\r
- endif\r
- endif\r
- if(contin)then\r
- j=j+1\r
- ifr2(k,nfr,j)=i\r
- Lfr2(k,nfr)=j\r
- else\r
- nfr=nfr+1\r
- j=1\r
- ifr2(k,nfr,j)=i\r
- Lfr2(k,nfr)=j\r
- endif\r
- enddo\r
- Lfr_max=0\r
- i_fr2(k)=1 !ID of the maximum piece\r
- do i=1,nfr\r
- if(Lfr_max.lt.Lfr2(k,i))then\r
- Lfr_max=Lfr2(k,i)\r
- i_fr2(k)=i\r
- endif\r
- enddo\r
- if(Lfr_max.lt.r_min)then\r
- dcu=1.1*dcu\r
- goto 20\r
- endif\r
- enddo\r
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- \r
-ccc select what piece will be used (this may araise ansysmetry, but\r
-ccc only when L1=L2 and Lfr1=Lfr2 and L1 ne Lfr1\r
-ccc if L1=Lfr1 and L2=Lfr2 (normal proteins), it will be the same as initial1\r
- mark=1\r
- if(Lfr2(1,i_fr2(1)).lt.Lfr2(2,i_fr2(2)))then\r
- mark=1\r
- elseif(Lfr2(1,i_fr2(1)).gt.Lfr2(2,i_fr2(2)))then\r
- mark=2\r
- else !Lfr1=Lfr2\r
- if(nseq1.le.nseq2)then\r
- mark=1\r
- else\r
- mark=2\r
- endif\r
- endif\r
-ccc \r
- L_fr=Lfr2(mark,i_fr2(mark))\r
- do i=1,L_fr\r
- ifr(i)=ifr2(mark,i_fr2(mark),i)\r
- enddo\r
-ccc \r
- if(mark.eq.1)then !non-redundant to get_initial1\r
- nseq0=nseq1\r
- else\r
- nseq0=nseq2\r
- endif\r
- if(L_fr.eq.nseq0)then\r
- n1=int(nseq0*0.1) !0\r
- n2=int(nseq0*0.89) !2\r
- j=0\r
- do i=n1,n2\r
- j=j+1\r
- ifr(j)=ifr(n1+j)\r
- enddo\r
- L_fr=j\r
- endif\r
- \r
-ccc get initial ------------->\r
- if(mark.eq.1)then !nseq1 as the smallest one\r
- nseq1_=L_fr\r
- aL=min(nseq1_,nseq2)\r
- idel=aL/2.5 !minimum size of considered fragment\r
- if(idel.le.fra_min1)idel=fra_min1\r
- n1=-nseq2+idel !shift1\r
- n2=nseq1_-idel !shift2\r
- GL_max=0\r
- do ishift=n1,n2\r
- L=0\r
- do j=1,nseq2\r
- i=j+ishift\r
- if(i.ge.1.and.i.le.nseq1_)then\r
- L=L+1\r
- invmap(j)=ifr(i)\r
- else\r
- invmap(j)=-1\r
- endif\r
- enddo\r
- if(L.ge.idel)then\r
- call get_GL(GL)\r
- if(GL.gt.GL_max)then\r
- GL_max=GL\r
- do i=1,nseq2\r
- invmap_i(i)=invmap(i)\r
- enddo\r
- endif\r
- endif\r
- enddo\r
- else !@@@@@@@@@@@@@@@@@@@@\r
- nseq2_=L_fr\r
- aL=min(nseq1,nseq2_)\r
- idel=aL/2.5 !minimum size of considered fragment\r
- if(idel.le.fra_min1)idel=fra_min1\r
- n1=-nseq2_+idel\r
- n2=nseq1-idel\r
- GL_max=0\r
- do ishift=n1,n2\r
- L=0\r
- do j=1,nseq2\r
- invmap(j)=-1\r
- enddo\r
- do j=1,nseq2_\r
- i=j+ishift\r
- if(i.ge.1.and.i.le.nseq1)then\r
- L=L+1\r
- invmap(ifr(j))=i\r
- endif\r
- enddo\r
- if(L.ge.idel)then\r
- call get_GL(GL)\r
- if(GL.gt.GL_max)then\r
- GL_max=GL\r
- do i=1,nseq2\r
- invmap_i(i)=invmap(i)\r
- enddo\r
- endif\r
- endif\r
- enddo\r
- endif\r
- \r
- return\r
- end\r
-\r
-**************************************************************\r
-* smooth the secondary structure assignment\r
-**************************************************************\r
- subroutine smooth\r
- PARAMETER(nmax=5000)\r
- common/sec/isec(nmax),jsec(nmax)\r
- common/length/nseq1,nseq2\r
-\r
-*** smooth single -------------->\r
-*** --x-- => -----\r
- do i=1,nseq1\r
- if(isec(i).eq.2.or.isec(i).eq.4)then\r
- j=isec(i)\r
- if(isec(i-2).ne.j)then\r
- if(isec(i-1).ne.j)then\r
- if(isec(i+1).ne.j)then\r
- if(isec(i+1).ne.j)then\r
- isec(i)=1\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- enddo\r
- do i=1,nseq2\r
- if(jsec(i).eq.2.or.jsec(i).eq.4)then\r
- j=jsec(i)\r
- if(jsec(i-2).ne.j)then\r
- if(jsec(i-1).ne.j)then\r
- if(jsec(i+1).ne.j)then\r
- if(jsec(i+1).ne.j)then\r
- jsec(i)=1\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- enddo\r
-\r
-*** smooth double -------------->\r
-*** --xx-- => ------\r
- do i=1,nseq1\r
- if(isec(i).ne.2)then\r
- if(isec(i+1).ne.2)then\r
- if(isec(i+2).eq.2)then\r
- if(isec(i+3).eq.2)then\r
- if(isec(i+4).ne.2)then\r
- if(isec(i+5).ne.2)then\r
- isec(i+2)=1\r
- isec(i+3)=1\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
-\r
- if(isec(i).ne.4)then\r
- if(isec(i+1).ne.4)then\r
- if(isec(i+2).eq.4)then\r
- if(isec(i+3).eq.4)then\r
- if(isec(i+4).ne.4)then\r
- if(isec(i+5).ne.4)then\r
- isec(i+2)=1\r
- isec(i+3)=1\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- enddo\r
- do i=1,nseq2\r
- if(jsec(i).ne.2)then\r
- if(jsec(i+1).ne.2)then\r
- if(jsec(i+2).eq.2)then\r
- if(jsec(i+3).eq.2)then\r
- if(jsec(i+4).ne.2)then\r
- if(jsec(i+5).ne.2)then\r
- jsec(i+2)=1\r
- jsec(i+3)=1\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
-\r
- if(jsec(i).ne.4)then\r
- if(jsec(i+1).ne.4)then\r
- if(jsec(i+2).eq.4)then\r
- if(jsec(i+3).eq.4)then\r
- if(jsec(i+4).ne.4)then\r
- if(jsec(i+5).ne.4)then\r
- jsec(i+2)=1\r
- jsec(i+3)=1\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- enddo\r
-\r
-*** connect -------------->\r
-*** x-x => xxx\r
- do i=1,nseq1\r
- if(isec(i).eq.2)then\r
- if(isec(i+1).ne.2)then\r
- if(isec(i+2).eq.2)then\r
- isec(i+1)=2\r
- endif\r
- endif\r
- endif\r
-\r
- if(isec(i).eq.4)then\r
- if(isec(i+1).ne.4)then\r
- if(isec(i+2).eq.4)then\r
- isec(i+1)=4\r
- endif\r
- endif\r
- endif\r
- enddo\r
- do i=1,nseq2\r
- if(jsec(i).eq.2)then\r
- if(jsec(i+1).ne.2)then\r
- if(jsec(i+2).eq.2)then\r
- jsec(i+1)=2\r
- endif\r
- endif\r
- endif\r
-\r
- if(jsec(i).eq.4)then\r
- if(jsec(i+1).ne.4)then\r
- if(jsec(i+2).eq.4)then\r
- jsec(i+1)=4\r
- endif\r
- endif\r
- endif\r
- enddo\r
-\r
- return\r
- end\r
-\r
-*************************************************************\r
-* assign secondary structure:\r
-*************************************************************\r
- function diszy(i,i1,i2)\r
- PARAMETER(nmax=5000)\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- diszy=sqrt((xa(1,i1,i)-xa(1,i2,i))**2\r
- & +(xa(2,i1,i)-xa(2,i2,i))**2\r
- & +(xa(3,i1,i)-xa(3,i2,i))**2)\r
- return\r
- end\r
-\r
-*************************************************************\r
-* assign secondary structure:\r
-*************************************************************\r
- function make_sec(dis13,dis14,dis15,dis24,dis25,dis35)\r
- make_sec=1\r
- delta=2.1\r
- if(abs(dis15-6.37).lt.delta)then\r
- if(abs(dis14-5.18).lt.delta)then\r
- if(abs(dis25-5.18).lt.delta)then\r
- if(abs(dis13-5.45).lt.delta)then\r
- if(abs(dis24-5.45).lt.delta)then\r
- if(abs(dis35-5.45).lt.delta)then\r
- make_sec=2 !helix\r
- return\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- delta=1.42\r
- if(abs(dis15-13).lt.delta)then\r
- if(abs(dis14-10.4).lt.delta)then\r
- if(abs(dis25-10.4).lt.delta)then\r
- if(abs(dis13-6.1).lt.delta)then\r
- if(abs(dis24-6.1).lt.delta)then\r
- if(abs(dis35-6.1).lt.delta)then\r
- make_sec=4 !strand\r
- return\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- endif\r
- if(dis15.lt.8)then\r
- make_sec=3\r
- endif\r
-\r
- return\r
- end\r
-\r
-****************************************************************\r
-* quickly calculate TM-score with given invmap(i) in 3 iterations\r
-****************************************************************\r
- subroutine get_GL(GL)\r
- PARAMETER(nmax=5000)\r
- common/length/nseq1,nseq2\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/d0/d0,anseq\r
- dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
- dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
- common/TM/TM,TMmax\r
- common/n1n2/n1(nmax),n2(nmax)\r
- common/d00/d00,d002\r
-\r
- dimension xo1(nmax),yo1(nmax),zo1(nmax)\r
- dimension xo2(nmax),yo2(nmax),zo2(nmax)\r
- dimension dis2(nmax)\r
-\r
-ccc RMSD:\r
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
- double precision u(3,3),t(3),rms,drms !armsd is real\r
- data w /nmax*1.0/\r
-ccc \r
-\r
-c calculate RMSD between aligned structures and rotate the structures -->\r
- n_al=0\r
- do j=1,NSEQ2\r
- i=invmap(j) !j aligned to i\r
- if(i.gt.0)then\r
- n_al=n_al+1\r
- r_1(1,n_al)=xa(1,i,0)\r
- r_1(2,n_al)=xa(2,i,0)\r
- r_1(3,n_al)=xa(3,i,0)\r
- r_2(1,n_al)=xa(1,j,1)\r
- r_2(2,n_al)=xa(2,j,1)\r
- r_2(3,n_al)=xa(3,j,1)\r
- xo1(n_al)=xa(1,i,0)\r
- yo1(n_al)=xa(2,i,0)\r
- zo1(n_al)=xa(3,i,0)\r
- xo2(n_al)=xa(1,j,1)\r
- yo2(n_al)=xa(2,j,1)\r
- zo2(n_al)=xa(3,j,1)\r
- endif\r
- enddo\r
- call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- GL=0\r
- do i=1,n_al\r
- xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)\r
- yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)\r
- zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)\r
- dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2\r
- GL=GL+1/(1+dis2(i)/(d0**2))\r
- enddo\r
-ccc for next iteration------------->\r
- d002t=d002\r
- 21 j=0\r
- do i=1,n_al\r
- if(dis2(i).le.d002t)then\r
- j=j+1\r
- r_1(1,j)=xo1(i)\r
- r_1(2,j)=yo1(i)\r
- r_1(3,j)=zo1(i)\r
- r_2(1,j)=xo2(i)\r
- r_2(2,j)=yo2(i)\r
- r_2(3,j)=zo2(i)\r
- endif\r
- enddo\r
- if(j.lt.3.and.n_al.gt.3)then\r
- d002t=d002t+.5\r
- goto 21\r
- endif\r
- L=j\r
- call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- G2=0\r
- do i=1,n_al\r
- xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)\r
- yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)\r
- zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)\r
- dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2\r
- G2=G2+1/(1+dis2(i)/(d0**2))\r
- enddo\r
-ccc for next iteration------------->\r
- d002t=d002+1\r
- 22 j=0\r
- do i=1,n_al\r
- if(dis2(i).le.d002t)then\r
- j=j+1\r
- r_1(1,j)=xo1(i)\r
- r_1(2,j)=yo1(i)\r
- r_1(3,j)=zo1(i)\r
- r_2(1,j)=xo2(i)\r
- r_2(2,j)=yo2(i)\r
- r_2(3,j)=zo2(i)\r
- endif\r
- enddo\r
- if(j.lt.3.and.n_al.gt.3)then\r
- d002t=d002t+.5\r
- goto 22\r
- endif\r
- L=j\r
- call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- G3=0\r
- do i=1,n_al\r
- xx=t(1)+u(1,1)*xo1(i)+u(1,2)*yo1(i)+u(1,3)*zo1(i)\r
- yy=t(2)+u(2,1)*xo1(i)+u(2,2)*yo1(i)+u(2,3)*zo1(i)\r
- zz=t(3)+u(3,1)*xo1(i)+u(3,2)*yo1(i)+u(3,3)*zo1(i)\r
- dis2(i)=(xx-xo2(i))**2+(yy-yo2(i))**2+(zz-zo2(i))**2\r
- G3=G3+1/(1+dis2(i)/(d0**2))\r
- enddo\r
- if(G2.gt.GL)GL=G2\r
- if(G3.gt.GL)GL=G3\r
-\r
-c^^^^^^^^^^^^^^^^ GL done ^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- end\r
-\r
-****************************************************************\r
-* with invmap(i) calculate TM-score and martix score(i,j) for rotation \r
-****************************************************************\r
- subroutine get_score\r
- PARAMETER(nmax=5000)\r
- common/length/nseq1,nseq2\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/d0/d0,anseq\r
- dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
- dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
- common/TM/TM,TMmax\r
- common/n1n2/n1(nmax),n2(nmax)\r
-\r
-ccc RMSD:\r
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
- double precision u(3,3),t(3),rms,drms !armsd is real\r
- data w /nmax*1.0/\r
-ccc \r
-\r
-c calculate RMSD between aligned structures and rotate the structures -->\r
- n_al=0\r
- do j=1,NSEQ2\r
- i=invmap(j) !j aligned to i\r
- if(i.gt.0)then\r
- n_al=n_al+1\r
-ccc for TM-score:\r
- xtm1(n_al)=xa(1,i,0) !for TM-score\r
- ytm1(n_al)=xa(2,i,0)\r
- ztm1(n_al)=xa(3,i,0)\r
- xtm2(n_al)=xa(1,j,1)\r
- ytm2(n_al)=xa(2,j,1)\r
- ztm2(n_al)=xa(3,j,1)\r
-ccc for rotation matrix:\r
- r_1(1,n_al)=xa(1,i,0)\r
- r_1(2,n_al)=xa(2,i,0)\r
- r_1(3,n_al)=xa(3,i,0)\r
- endif\r
- enddo\r
-*** calculate TM-score for the given alignment----------->\r
- d0_input=d0\r
- call TMscore8_search(d0_input,n_al,xtm1,ytm1,ztm1,n1,\r
- & n_al,xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !simplified search engine\r
- TM=TM*n_al/anseq !TM-score\r
-*** calculate score matrix score(i,j)------------------>\r
- do i=1,n_al\r
- r_2(1,i)=xtm1(i)\r
- r_2(2,i)=ytm1(i)\r
- r_2(3,i)=ztm1(i)\r
- enddo\r
- call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- do i=1,nseq1\r
- 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
- 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
- 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
- do j=1,nseq2\r
- dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2\r
- score(i,j)=1/(1+dd/d0**2)\r
- enddo\r
- enddo\r
- \r
-c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- end\r
-\r
-****************************************************************\r
-* with invmap(i) calculate score(i,j) using RMSD rotation \r
-****************************************************************\r
- subroutine get_score1\r
- PARAMETER(nmax=5000)\r
- common/length/nseq1,nseq2\r
- COMMON/BACKBONE/XA(3,nmax,0:1)\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- common/zscore/zrms,n_al,rmsd_al\r
- common/d0/d0,anseq\r
- common/d0min/d0_min\r
- dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)\r
- dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)\r
- common/TM/TM,TMmax\r
- common/n1n2/n1(nmax),n2(nmax)\r
-\r
-ccc RMSD:\r
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
- double precision u(3,3),t(3),rms,drms !armsd is real\r
- data w /nmax*1.0/\r
-ccc \r
-\r
-c calculate RMSD between aligned structures and rotate the structures -->\r
- n_al=0\r
- do j=1,NSEQ2\r
- i=invmap(j) !j aligned to i\r
- if(i.gt.0)then\r
- n_al=n_al+1\r
-ccc for rotation matrix:\r
- r_1(1,n_al)=xa(1,i,0)\r
- r_1(2,n_al)=xa(2,i,0)\r
- r_1(3,n_al)=xa(3,i,0)\r
- r_2(1,n_al)=xa(1,j,1)\r
- r_2(2,n_al)=xa(2,j,1)\r
- r_2(3,n_al)=xa(3,j,1)\r
- endif\r
- enddo\r
-*** calculate score matrix score(i,j)------------------>\r
- call u3b(w,r_1,r_2,n_al,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- d01=d0+1.5\r
- if(d01.lt.d0_min)d01=d0_min\r
- d02=d01*d01\r
- do i=1,nseq1\r
- 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
- 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
- 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
- do j=1,nseq2\r
- dd=(xx-xa(1,j,1))**2+(yy-xa(2,j,1))**2+(zz-xa(3,j,1))**2\r
- score(i,j)=1/(1+dd/d02)\r
- enddo\r
- enddo\r
-\r
-c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- end\r
-\r
-\r
-*************************************************************************\r
-*************************************************************************\r
-* This is a subroutine to compare two structures and find the \r
-* superposition that has the maximum TM-score.\r
-*\r
-* L1--Length of the first structure\r
-* (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure\r
-* n1(i)--Residue sequence number of i'th residue at the first structure\r
-* L2--Length of the second structure\r
-* (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure\r
-* n2(i)--Residue sequence number of i'th residue at the second structure\r
-* TM--TM-score of the comparison\r
-* Rcomm--RMSD of two structures in the common aligned residues\r
-* Lcomm--Length of the common aligned regions\r
-*\r
-* Note: \r
-* 1, Always put native as the second structure, by which TM-score\r
-* is normalized.\r
-* 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after\r
-* TM-score superposition.\r
-*************************************************************************\r
-*************************************************************************\r
-*** dis<8, simplified search engine\r
- subroutine TMscore8_search(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,\r
- & TM,Rcomm,Lcomm)\r
- PARAMETER(nmax=5000)\r
- common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)\r
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
- common/para/d,d0\r
- common/d0min/d0_min\r
- common/align/n_ali,iA(nmax),iB(nmax)\r
- common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
- dimension k_ali(nmax),k_ali0(nmax)\r
- dimension L_ini(100),iq(nmax)\r
- common/scores/score\r
- double precision score,score_max\r
- dimension xa(nmax),ya(nmax),za(nmax)\r
- dimension iL0(nmax)\r
-\r
- dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)\r
- dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)\r
-\r
-ccc RMSD:\r
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
- double precision u(3,3),t(3),rms,drms !armsd is real\r
- data w /nmax*1.0/\r
-ccc \r
-\r
-********* convert input data ***************************\r
-* because L1=L2 in this special case---------->\r
- nseqA=L1\r
- nseqB=L2\r
- do i=1,nseqA\r
- xa(i)=x1(i)\r
- ya(i)=y1(i)\r
- za(i)=z1(i)\r
- nresA(i)=n1(i)\r
- xb(i)=x2(i)\r
- yb(i)=y2(i)\r
- zb(i)=z2(i)\r
- nresB(i)=n2(i)\r
- iA(i)=i\r
- iB(i)=i\r
- enddo\r
- n_ali=L1 !number of aligned residues\r
- Lcomm=L1\r
-\r
-************/////\r
-* parameters:\r
-*****************\r
-*** d0------------->\r
- d0=dx\r
- if(d0.lt.d0_min)d0=d0_min\r
-*** d0_search ----->\r
- d0_search=d0\r
- if(d0_search.gt.8)d0_search=8\r
- if(d0_search.lt.4.5)d0_search=4.5\r
-*** iterative parameters ----->\r
- n_it=20 !maximum number of iterations\r
- d_output=5 !for output alignment\r
- n_init_max=6 !maximum number of L_init\r
- n_init=0\r
- L_ini_min=4\r
- if(n_ali.lt.4)L_ini_min=n_ali\r
- do i=1,n_init_max-1\r
- n_init=n_init+1\r
- L_ini(n_init)=n_ali/2**(n_init-1)\r
- if(L_ini(n_init).le.L_ini_min)then\r
- L_ini(n_init)=L_ini_min\r
- goto 402\r
- endif\r
- enddo\r
- n_init=n_init+1\r
- L_ini(n_init)=L_ini_min\r
- 402 continue\r
-\r
-******************************************************************\r
-* find the maximum score starting from local structures superposition\r
-******************************************************************\r
- score_max=-1 !TM-score\r
- do 333 i_init=1,n_init\r
- L_init=L_ini(i_init)\r
- iL_max=n_ali-L_init+1\r
- k=0\r
- do i=1,iL_max,40 !this is the simplification!\r
- k=k+1\r
- iL0(k)=i\r
- enddo\r
- if(iL0(k).lt.iL_max)then\r
- k=k+1\r
- iL0(k)=iL_max\r
- endif\r
- n_shift=k\r
- do 300 i_shift=1,n_shift\r
- iL=iL0(i_shift)\r
- LL=0\r
- ka=0\r
- do i=1,L_init\r
- k=iL+i-1 ![1,n_ali] common aligned\r
- r_1(1,i)=xa(iA(k))\r
- r_1(2,i)=ya(iA(k))\r
- r_1(3,i)=za(iA(k))\r
- r_2(1,i)=xb(iB(k))\r
- r_2(2,i)=yb(iB(k))\r
- r_2(3,i)=zb(iB(k))\r
- LL=LL+1\r
- ka=ka+1\r
- k_ali(ka)=k\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- if(i_init.eq.1)then !global superposition\r
- armsd=dsqrt(rms/LL)\r
- Rcomm=armsd\r
- endif\r
- do j=1,nseqA\r
- xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- d=d0_search-1\r
- call score_fun8 !init, get scores, n_cut+i_ali(i) for iteration\r
- if(score_max.lt.score)then\r
- score_max=score\r
- ka0=ka\r
- do i=1,ka0\r
- k_ali0(i)=k_ali(i)\r
- enddo\r
- endif\r
-*** iteration for extending ---------------------------------->\r
- d=d0_search+1\r
- do 301 it=1,n_it\r
- LL=0\r
- ka=0\r
- do i=1,n_cut\r
- m=i_ali(i) ![1,n_ali]\r
- r_1(1,i)=xa(iA(m))\r
- r_1(2,i)=ya(iA(m))\r
- r_1(3,i)=za(iA(m))\r
- r_2(1,i)=xb(iB(m))\r
- r_2(2,i)=yb(iB(m))\r
- r_2(3,i)=zb(iB(m))\r
- ka=ka+1\r
- k_ali(ka)=m\r
- LL=LL+1\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- do j=1,nseqA\r
- xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- call score_fun8 !get scores, n_cut+i_ali(i) for iteration\r
- if(score_max.lt.score)then\r
- score_max=score\r
- ka0=ka\r
- do i=1,ka\r
- k_ali0(i)=k_ali(i)\r
- enddo\r
- endif\r
- if(it.eq.n_it)goto 302\r
- if(n_cut.eq.ka)then\r
- neq=0\r
- do i=1,n_cut\r
- if(i_ali(i).eq.k_ali(i))neq=neq+1\r
- enddo\r
- if(n_cut.eq.neq)goto 302\r
- endif\r
- 301 continue !for iteration\r
- 302 continue\r
- 300 continue !for shift\r
- 333 continue !for initial length, L_ali/M\r
-\r
-******** return the final rotation ****************\r
- LL=0\r
- do i=1,ka0\r
- m=k_ali0(i) !record of the best alignment\r
- r_1(1,i)=xa(iA(m))\r
- r_1(2,i)=ya(iA(m))\r
- r_1(3,i)=za(iA(m))\r
- r_2(1,i)=xb(iB(m))\r
- r_2(2,i)=yb(iB(m))\r
- r_2(3,i)=zb(iB(m))\r
- LL=LL+1\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- do j=1,nseqA\r
- x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- TM=score_max\r
-\r
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- END\r
-\r
-*************************************************************************\r
-*************************************************************************\r
-* This is a subroutine to compare two structures and find the \r
-* superposition that has the maximum TM-score.\r
-*\r
-* L1--Length of the first structure\r
-* (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure\r
-* n1(i)--Residue sequence number of i'th residue at the first structure\r
-* L2--Length of the second structure\r
-* (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure\r
-* n2(i)--Residue sequence number of i'th residue at the second structure\r
-* TM--TM-score of the comparison\r
-* Rcomm--RMSD of two structures in the common aligned residues\r
-* Lcomm--Length of the common aligned regions\r
-*\r
-* Note: \r
-* 1, Always put native as the second structure, by which TM-score\r
-* is normalized.\r
-* 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after\r
-* TM-score superposition.\r
-*************************************************************************\r
-*************************************************************************\r
-*** dis<8, but same search engine\r
- subroutine TMscore8(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,\r
- & TM,Rcomm,Lcomm)\r
- PARAMETER(nmax=5000)\r
- common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)\r
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
- common/para/d,d0\r
- common/d0min/d0_min\r
- common/align/n_ali,iA(nmax),iB(nmax)\r
- common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
- dimension k_ali(nmax),k_ali0(nmax)\r
- dimension L_ini(100),iq(nmax)\r
- common/scores/score\r
- double precision score,score_max\r
- dimension xa(nmax),ya(nmax),za(nmax)\r
-\r
- dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)\r
- dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)\r
-\r
-ccc RMSD:\r
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
- double precision u(3,3),t(3),rms,drms !armsd is real\r
- data w /nmax*1.0/\r
-ccc \r
-\r
-********* convert input data ***************************\r
-* because L1=L2 in this special case---------->\r
- nseqA=L1\r
- nseqB=L2\r
- do i=1,nseqA\r
- xa(i)=x1(i)\r
- ya(i)=y1(i)\r
- za(i)=z1(i)\r
- nresA(i)=n1(i)\r
- xb(i)=x2(i)\r
- yb(i)=y2(i)\r
- zb(i)=z2(i)\r
- nresB(i)=n2(i)\r
- iA(i)=i\r
- iB(i)=i\r
- enddo\r
- n_ali=L1 !number of aligned residues\r
- Lcomm=L1\r
-\r
-************/////\r
-* parameters:\r
-*****************\r
-*** d0------------->\r
- d0=dx\r
- if(d0.lt.d0_min)d0=d0_min\r
-*** d0_search ----->\r
- d0_search=d0\r
- if(d0_search.gt.8)d0_search=8\r
- if(d0_search.lt.4.5)d0_search=4.5\r
-*** iterative parameters ----->\r
- n_it=20 !maximum number of iterations\r
- d_output=5 !for output alignment\r
- n_init_max=6 !maximum number of L_init\r
- n_init=0\r
- L_ini_min=4\r
- if(n_ali.lt.4)L_ini_min=n_ali\r
- do i=1,n_init_max-1\r
- n_init=n_init+1\r
- L_ini(n_init)=n_ali/2**(n_init-1)\r
- if(L_ini(n_init).le.L_ini_min)then\r
- L_ini(n_init)=L_ini_min\r
- goto 402\r
- endif\r
- enddo\r
- n_init=n_init+1\r
- L_ini(n_init)=L_ini_min\r
- 402 continue\r
-\r
-******************************************************************\r
-* find the maximum score starting from local structures superposition\r
-******************************************************************\r
- score_max=-1 !TM-score\r
- do 333 i_init=1,n_init\r
- L_init=L_ini(i_init)\r
- iL_max=n_ali-L_init+1\r
- do 300 iL=1,iL_max !on aligned residues, [1,nseqA]\r
- LL=0\r
- ka=0\r
- do i=1,L_init\r
- k=iL+i-1 ![1,n_ali] common aligned\r
- r_1(1,i)=xa(iA(k))\r
- r_1(2,i)=ya(iA(k))\r
- r_1(3,i)=za(iA(k))\r
- r_2(1,i)=xb(iB(k))\r
- r_2(2,i)=yb(iB(k))\r
- r_2(3,i)=zb(iB(k))\r
- LL=LL+1\r
- ka=ka+1\r
- k_ali(ka)=k\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- if(i_init.eq.1)then !global superposition\r
- armsd=dsqrt(rms/LL)\r
- Rcomm=armsd\r
- endif\r
- do j=1,nseqA\r
- xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- d=d0_search-1\r
- call score_fun8 !init, get scores, n_cut+i_ali(i) for iteration\r
- if(score_max.lt.score)then\r
- score_max=score\r
- ka0=ka\r
- do i=1,ka0\r
- k_ali0(i)=k_ali(i)\r
- enddo\r
- endif\r
-*** iteration for extending ---------------------------------->\r
- d=d0_search+1\r
- do 301 it=1,n_it\r
- LL=0\r
- ka=0\r
- do i=1,n_cut\r
- m=i_ali(i) ![1,n_ali]\r
- r_1(1,i)=xa(iA(m))\r
- r_1(2,i)=ya(iA(m))\r
- r_1(3,i)=za(iA(m))\r
- r_2(1,i)=xb(iB(m))\r
- r_2(2,i)=yb(iB(m))\r
- r_2(3,i)=zb(iB(m))\r
- ka=ka+1\r
- k_ali(ka)=m\r
- LL=LL+1\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- do j=1,nseqA\r
- xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- call score_fun8 !get scores, n_cut+i_ali(i) for iteration\r
- if(score_max.lt.score)then\r
- score_max=score\r
- ka0=ka\r
- do i=1,ka\r
- k_ali0(i)=k_ali(i)\r
- enddo\r
- endif\r
- if(it.eq.n_it)goto 302\r
- if(n_cut.eq.ka)then\r
- neq=0\r
- do i=1,n_cut\r
- if(i_ali(i).eq.k_ali(i))neq=neq+1\r
- enddo\r
- if(n_cut.eq.neq)goto 302\r
- endif\r
- 301 continue !for iteration\r
- 302 continue\r
- 300 continue !for shift\r
- 333 continue !for initial length, L_ali/M\r
-\r
-******** return the final rotation ****************\r
- LL=0\r
- do i=1,ka0\r
- m=k_ali0(i) !record of the best alignment\r
- r_1(1,i)=xa(iA(m))\r
- r_1(2,i)=ya(iA(m))\r
- r_1(3,i)=za(iA(m))\r
- r_2(1,i)=xb(iB(m))\r
- r_2(2,i)=yb(iB(m))\r
- r_2(3,i)=zb(iB(m))\r
- LL=LL+1\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- do j=1,nseqA\r
- x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- TM=score_max\r
-\r
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- END\r
-\r
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
-c 1, collect those residues with dis<d;\r
-c 2, calculate score_GDT, score_maxsub, score_TM\r
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
- subroutine score_fun8\r
- PARAMETER(nmax=5000)\r
-\r
- common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)\r
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
- common/para/d,d0\r
- common/align/n_ali,iA(nmax),iB(nmax)\r
- common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
- common/scores/score\r
- double precision score,score_max\r
- common/d8/d8\r
-\r
- d_tmp=d\r
- 21 n_cut=0 !number of residue-pairs dis<d, for iteration\r
- score_sum=0 !TMscore\r
- do k=1,n_ali\r
- i=iA(k) ![1,nseqA] reoder number of structureA\r
- j=iB(k) ![1,nseqB]\r
- dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)\r
- if(dis.lt.d_tmp)then\r
- n_cut=n_cut+1\r
- i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d\r
- endif\r
- if(dis.le.d8)then\r
- score_sum=score_sum+1/(1+(dis/d0)**2)\r
- endif\r
- enddo\r
- if(n_cut.lt.3.and.n_ali.gt.3)then\r
- d_tmp=d_tmp+.5\r
- goto 21\r
- endif\r
- score=score_sum/float(nseqB) !TM-score\r
-\r
- return\r
- end\r
-\r
-*************************************************************************\r
-*************************************************************************\r
-* This is a subroutine to compare two structures and find the \r
-* superposition that has the maximum TM-score.\r
-*\r
-* L1--Length of the first structure\r
-* (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure\r
-* n1(i)--Residue sequence number of i'th residue at the first structure\r
-* L2--Length of the second structure\r
-* (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure\r
-* n2(i)--Residue sequence number of i'th residue at the second structure\r
-* TM--TM-score of the comparison\r
-* Rcomm--RMSD of two structures in the common aligned residues\r
-* Lcomm--Length of the common aligned regions\r
-*\r
-* Note: \r
-* 1, Always put native as the second structure, by which TM-score\r
-* is normalized.\r
-* 2, The returned (x1(i),y1(i),z1(i)) are the rotated structure after\r
-* TM-score superposition.\r
-*************************************************************************\r
-*************************************************************************\r
-*** normal TM-score:\r
- subroutine TMscore(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,\r
- & TM,Rcomm,Lcomm)\r
- PARAMETER(nmax=5000)\r
- common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)\r
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
- common/para/d,d0\r
- common/d0min/d0_min\r
- common/align/n_ali,iA(nmax),iB(nmax)\r
- common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
- dimension k_ali(nmax),k_ali0(nmax)\r
- dimension L_ini(100),iq(nmax)\r
- common/scores/score\r
- double precision score,score_max\r
- dimension xa(nmax),ya(nmax),za(nmax)\r
-\r
- dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)\r
- dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)\r
-\r
-ccc RMSD:\r
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)\r
- double precision u(3,3),t(3),rms,drms !armsd is real\r
- data w /nmax*1.0/\r
-ccc \r
-\r
-********* convert input data ***************************\r
-* because L1=L2 in this special case---------->\r
- nseqA=L1\r
- nseqB=L2\r
- do i=1,nseqA\r
- xa(i)=x1(i)\r
- ya(i)=y1(i)\r
- za(i)=z1(i)\r
- nresA(i)=n1(i)\r
- xb(i)=x2(i)\r
- yb(i)=y2(i)\r
- zb(i)=z2(i)\r
- nresB(i)=n2(i)\r
- iA(i)=i\r
- iB(i)=i\r
- enddo\r
- n_ali=L1 !number of aligned residues\r
- Lcomm=L1\r
-\r
-************/////\r
-* parameters:\r
-*****************\r
-*** d0------------->\r
-c d0=1.24*(nseqB-15)**(1.0/3.0)-1.8\r
- d0=dx\r
- if(d0.lt.d0_min)d0=d0_min\r
-*** d0_search ----->\r
- d0_search=d0\r
- if(d0_search.gt.8)d0_search=8\r
- if(d0_search.lt.4.5)d0_search=4.5\r
-*** iterative parameters ----->\r
- n_it=20 !maximum number of iterations\r
- d_output=5 !for output alignment\r
- n_init_max=6 !maximum number of L_init\r
- n_init=0\r
- L_ini_min=4\r
- if(n_ali.lt.4)L_ini_min=n_ali\r
- do i=1,n_init_max-1\r
- n_init=n_init+1\r
- L_ini(n_init)=n_ali/2**(n_init-1)\r
- if(L_ini(n_init).le.L_ini_min)then\r
- L_ini(n_init)=L_ini_min\r
- goto 402\r
- endif\r
- enddo\r
- n_init=n_init+1\r
- L_ini(n_init)=L_ini_min\r
- 402 continue\r
-\r
-******************************************************************\r
-* find the maximum score starting from local structures superposition\r
-******************************************************************\r
- score_max=-1 !TM-score\r
- do 333 i_init=1,n_init\r
- L_init=L_ini(i_init)\r
- iL_max=n_ali-L_init+1\r
- do 300 iL=1,iL_max !on aligned residues, [1,nseqA]\r
- LL=0\r
- ka=0\r
- do i=1,L_init\r
- k=iL+i-1 ![1,n_ali] common aligned\r
- r_1(1,i)=xa(iA(k))\r
- r_1(2,i)=ya(iA(k))\r
- r_1(3,i)=za(iA(k))\r
- r_2(1,i)=xb(iB(k))\r
- r_2(2,i)=yb(iB(k))\r
- r_2(3,i)=zb(iB(k))\r
- LL=LL+1\r
- ka=ka+1\r
- k_ali(ka)=k\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- if(i_init.eq.1)then !global superposition\r
- armsd=dsqrt(rms/LL)\r
- Rcomm=armsd\r
- endif\r
- do j=1,nseqA\r
- xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- d=d0_search-1\r
- call score_fun !init, get scores, n_cut+i_ali(i) for iteration\r
- if(score_max.lt.score)then\r
- score_max=score\r
- ka0=ka\r
- do i=1,ka0\r
- k_ali0(i)=k_ali(i)\r
- enddo\r
- endif\r
-*** iteration for extending ---------------------------------->\r
- d=d0_search+1\r
- do 301 it=1,n_it\r
- LL=0\r
- ka=0\r
- do i=1,n_cut\r
- m=i_ali(i) ![1,n_ali]\r
- r_1(1,i)=xa(iA(m))\r
- r_1(2,i)=ya(iA(m))\r
- r_1(3,i)=za(iA(m))\r
- r_2(1,i)=xb(iB(m))\r
- r_2(2,i)=yb(iB(m))\r
- r_2(3,i)=zb(iB(m))\r
- ka=ka+1\r
- k_ali(ka)=m\r
- LL=LL+1\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- do j=1,nseqA\r
- xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- call score_fun !get scores, n_cut+i_ali(i) for iteration\r
- if(score_max.lt.score)then\r
- score_max=score\r
- ka0=ka\r
- do i=1,ka\r
- k_ali0(i)=k_ali(i)\r
- enddo\r
- endif\r
- if(it.eq.n_it)goto 302\r
- if(n_cut.eq.ka)then\r
- neq=0\r
- do i=1,n_cut\r
- if(i_ali(i).eq.k_ali(i))neq=neq+1\r
- enddo\r
- if(n_cut.eq.neq)goto 302\r
- endif\r
- 301 continue !for iteration\r
- 302 continue\r
- 300 continue !for shift\r
- 333 continue !for initial length, L_ali/M\r
-\r
-******** return the final rotation ****************\r
- LL=0\r
- do i=1,ka0\r
- m=k_ali0(i) !record of the best alignment\r
- r_1(1,i)=xa(iA(m))\r
- r_1(2,i)=ya(iA(m))\r
- r_1(3,i)=za(iA(m))\r
- r_2(1,i)=xb(iB(m))\r
- r_2(2,i)=yb(iB(m))\r
- r_2(3,i)=zb(iB(m))\r
- LL=LL+1\r
- enddo\r
- call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2\r
- do j=1,nseqA\r
- x1(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)\r
- y1(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)\r
- z1(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)\r
- enddo\r
- TM=score_max\r
-\r
-c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\r
- return\r
- END\r
-\r
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
-c 1, collect those residues with dis<d;\r
-c 2, calculate score_GDT, score_maxsub, score_TM\r
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
- subroutine score_fun\r
- PARAMETER(nmax=5000)\r
-\r
- common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)\r
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB\r
- common/para/d,d0\r
- common/align/n_ali,iA(nmax),iB(nmax)\r
- common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score\r
- common/scores/score\r
- double precision score\r
-\r
- d_tmp=d\r
- 21 n_cut=0 !number of residue-pairs dis<d, for iteration\r
- score_sum=0 !TMscore\r
- do k=1,n_ali\r
- i=iA(k) ![1,nseqA] reoder number of structureA\r
- j=iB(k) ![1,nseqB]\r
- dis=sqrt((xa(i)-xb(j))**2+(ya(i)-yb(j))**2+(za(i)-zb(j))**2)\r
- if(dis.lt.d_tmp)then\r
- n_cut=n_cut+1\r
- i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d\r
- endif\r
- score_sum=score_sum+1/(1+(dis/d0)**2)\r
- enddo\r
- if(n_cut.lt.3.and.n_ali.gt.3)then\r
- d_tmp=d_tmp+.5\r
- goto 21\r
- endif\r
- score=score_sum/float(nseqB) !TM-score\r
-\r
- return\r
- end\r
-\r
-********************************************************************\r
-* Dynamic programming for alignment.\r
-* Input: score(i,j), and gap_open\r
-* Output: invmap(j)\r
-* \r
-* Please note this subroutine is not a correct implementation of \r
-* the N-W dynamic programming because the score tracks back only \r
-* one layer of the matrix. This code was exploited in TM-align \r
-* because it is about 1.5 times faster than a complete N-W code\r
-* and does not influence much the final structure alignment result.\r
-********************************************************************\r
- SUBROUTINE DP(NSEQ1,NSEQ2)\r
- PARAMETER(nmax=5000)\r
- LOGICAL*1 DIR\r
- common/dpc/score(nmax,nmax),gap_open,invmap(nmax)\r
- dimension DIR(0:nmax,0:nmax),VAL(0:nmax,0:nmax)\r
- REAL H,V\r
- \r
-*** initialize the matrix:\r
- val(0,0)=0\r
- do i=1,nseq1\r
- dir(i,0)=.false.\r
- val(i,0)=0\r
- enddo\r
- do j=1,nseq2\r
- dir(0,j)=.false.\r
- val(0,j)=0\r
- invmap(j)=-1\r
- enddo\r
-\r
-*** decide matrix and path:\r
- DO j=1,NSEQ2\r
- DO i=1,NSEQ1\r
- D=VAL(i-1,j-1)+SCORE(i,j)\r
- H=VAL(i-1,j)\r
- if(DIR(i-1,j))H=H+GAP_OPEN\r
- V=VAL(i,j-1)\r
- if(DIR(i,j-1))V=V+GAP_OPEN\r
- \r
- IF((D.GE.H).AND.(D.GE.V)) THEN\r
- DIR(I,J)=.true.\r
- VAL(i,j)=D\r
- ELSE\r
- DIR(I,J)=.false.\r
- if(V.GE.H)then\r
- val(i,j)=v\r
- else\r
- val(i,j)=h\r
- end if\r
- ENDIF\r
- ENDDO\r
- ENDDO\r
- \r
-*** extract the alignment:\r
- i=NSEQ1\r
- j=NSEQ2\r
- DO WHILE((i.GT.0).AND.(j.GT.0))\r
- IF(DIR(i,j))THEN\r
- invmap(j)=i\r
- i=i-1\r
- j=j-1\r
- ELSE\r
- H=VAL(i-1,j)\r
- if(DIR(i-1,j))H=H+GAP_OPEN\r
- V=VAL(i,j-1)\r
- if(DIR(i,j-1))V=V+GAP_OPEN\r
- IF(V.GE.H) THEN\r
- j=j-1\r
- ELSE\r
- i=i-1\r
- ENDIF\r
- ENDIF\r
- ENDDO\r
- \r
-c^^^^^^^^^^^^^^^Dynamical programming done ^^^^^^^^^^^^^^^^^^^\r
- return\r
- END\r
-\r
-cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc\r
-c w - w(m) is weight for atom pair c m (given)\r
-c x - x(i,m) are coordinates of atom c m in set x (given)\r
-c y - y(i,m) are coordinates of atom c m in set y (given)\r
-c n - n is number of atom pairs (given)\r
-c mode - 0:calculate rms only (given)\r
-c 1:calculate rms,u,t (takes longer)\r
-c rms - sum of w*(ux+t-y)**2 over all atom pairs (result)\r
-c u - u(i,j) is rotation matrix for best superposition (result)\r
-c t - t(i) is translation vector for best superposition (result)\r
-c ier - 0: a unique optimal superposition has been determined(result)\r
-c -1: superposition is not unique but optimal\r
-c -2: no result obtained because of negative weights w\r
-c or all weights equal to zero.\r
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc\r
- subroutine u3b(w, x, y, n, mode, rms, u, t, ier)\r
- integer ip(9), ip2312(4), i, j, k, l, m1, m, ier, n, mode\r
- double precision w(1), x(3, 1), y(3, 1), u(3, 3), t(3), rms, sigma\r
- double precision r(3, 3), xc(3), yc(3), wc, a(3, 3), b(3, 3), e0, \r
- &e(3), e1, e2, e3, d, spur, det, cof, h, g, cth, sth, sqrth, p, tol\r
- &, rr(6), rr1, rr2, rr3, rr4, rr5, rr6, ss(6), ss1, ss2, ss3, ss4, \r
- &ss5, ss6, zero, one, two, three, sqrt3\r
- equivalence (rr(1), rr1), (rr(2), rr2), (rr(3), rr3), (rr(4), rr4)\r
- &, (rr(5), rr5), (rr(6), rr6), (ss(1), ss1), (ss(2), ss2), (ss(3), \r
- &ss3), (ss(4), ss4), (ss(5), ss5), (ss(6), ss6), (e(1), e1), (e(2)\r
- &, e2), (e(3), e3)\r
- data sqrt3 / 1.73205080756888d+00 /\r
- data tol / 1.0d-2 /\r
- data zero / 0.0d+00 /\r
- data one / 1.0d+00 /\r
- data two / 2.0d+00 /\r
- data three / 3.0d+00 /\r
- data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /\r
- data ip2312 / 2, 3, 1, 2 /\r
-c 156 "rms.for"\r
- wc = zero\r
- rms = 0.0\r
- e0 = zero\r
- do 1 i = 1, 3\r
- xc(i) = zero\r
- yc(i) = zero\r
- t(i) = 0.0\r
- do 1 j = 1, 3\r
- d = zero\r
- if (i .eq. j) d = one\r
- u(i,j) = d\r
- a(i,j) = d\r
- 1 r(i,j) = zero\r
- ier = -1\r
-c**** DETERMINE CENTROIDS OF BOTH VECTOR SETS X AND Y\r
-c 170 "rms.for"\r
- if (n .lt. 1) return \r
-c 172 "rms.for"\r
- ier = -2\r
- do 2 m = 1, n\r
- if (w(m) .lt. 0.0) return \r
- wc = wc + w(m)\r
- do 2 i = 1, 3\r
- xc(i) = xc(i) + (w(m) * x(i,m))\r
- 2 yc(i) = yc(i) + (w(m) * y(i,m))\r
- if (wc .le. zero) return \r
- do 3 i = 1, 3\r
- xc(i) = xc(i) / wc\r
-c**** DETERMINE CORRELATION MATRIX R BETWEEN VECTOR SETS Y AND X\r
-c 182 "rms.for"\r
- 3 yc(i) = yc(i) / wc\r
-c 184 "rms.for"\r
- do 4 m = 1, n\r
- do 4 i = 1, 3\r
- e0 = e0 + (w(m) * (((x(i,m) - xc(i)) ** 2) + ((y(i,m) - yc(i)) ** \r
- &2)))\r
-c 187 "rms.for"\r
- d = w(m) * (y(i,m) - yc(i))\r
- do 4 j = 1, 3\r
-c**** CALCULATE DETERMINANT OF R(I,J)\r
-c 189 "rms.for"\r
- 4 r(i,j) = r(i,j) + (d * (x(j,m) - xc(j)))\r
-c 191 "rms.for"\r
- det = ((r(1,1) * ((r(2,2) * r(3,3)) - (r(2,3) * r(3,2)))) - (r(1,2\r
- &) * ((r(2,1) * r(3,3)) - (r(2,3) * r(3,1))))) + (r(1,3) * ((r(2,1)\r
- & * r(3,2)) - (r(2,2) * r(3,1))))\r
-c**** FORM UPPER TRIANGLE OF TRANSPOSED(R)*R\r
-c 194 "rms.for"\r
- sigma = det\r
-c 196 "rms.for"\r
- m = 0\r
- do 5 j = 1, 3\r
- do 5 i = 1, j\r
- m = m + 1\r
-c***************** EIGENVALUES *****************************************\r
-c**** FORM CHARACTERISTIC CUBIC X**3-3*SPUR*X**2+3*COF*X-DET=0\r
-c 200 "rms.for"\r
- 5 rr(m) = ((r(1,i) * r(1,j)) + (r(2,i) * r(2,j))) + (r(3,i) * r(3,j)\r
- &)\r
-c 203 "rms.for"\r
- spur = ((rr1 + rr3) + rr6) / three\r
- cof = ((((((rr3 * rr6) - (rr5 * rr5)) + (rr1 * rr6)) - (rr4 * rr4)\r
- &) + (rr1 * rr3)) - (rr2 * rr2)) / three\r
-c 205 "rms.for"\r
- det = det * det\r
- do 6 i = 1, 3\r
- 6 e(i) = spur\r
-c**** REDUCE CUBIC TO STANDARD FORM Y**3-3HY+2G=0 BY PUTTING X=Y+SPUR\r
-c 208 "rms.for"\r
- if (spur .le. zero) goto 40\r
-c 210 "rms.for"\r
- d = spur * spur\r
- h = d - cof\r
-c**** SOLVE CUBIC. ROOTS ARE E1,E2,E3 IN DECREASING ORDER\r
-c 212 "rms.for"\r
- g = (((spur * cof) - det) / two) - (spur * h)\r
-c 214 "rms.for"\r
- if (h .le. zero) goto 8\r
- sqrth = dsqrt(h)\r
- d = ((h * h) * h) - (g * g)\r
- if (d .lt. zero) d = zero\r
- d = datan2(dsqrt(d),- g) / three\r
- cth = sqrth * dcos(d)\r
- sth = (sqrth * sqrt3) * dsin(d)\r
- e1 = (spur + cth) + cth\r
- e2 = (spur - cth) + sth\r
- e3 = (spur - cth) - sth\r
-c.....HANDLE SPECIAL CASE OF 3 IDENTICAL ROOTS\r
-c 224 "rms.for"\r
- if (mode) 10, 50, 10\r
-c**************** EIGENVECTORS *****************************************\r
-c 226 "rms.for"\r
- 8 if (mode) 30, 50, 30\r
-c 228 "rms.for"\r
- 10 do 15 l = 1, 3, 2\r
- d = e(l)\r
- ss1 = ((d - rr3) * (d - rr6)) - (rr5 * rr5)\r
- ss2 = ((d - rr6) * rr2) + (rr4 * rr5)\r
- ss3 = ((d - rr1) * (d - rr6)) - (rr4 * rr4)\r
- ss4 = ((d - rr3) * rr4) + (rr2 * rr5)\r
- ss5 = ((d - rr1) * rr5) + (rr2 * rr4)\r
- ss6 = ((d - rr1) * (d - rr3)) - (rr2 * rr2)\r
- j = 1\r
- if (dabs(ss1) .ge. dabs(ss3)) goto 12\r
- j = 2\r
- if (dabs(ss3) .ge. dabs(ss6)) goto 13\r
- 11 j = 3\r
- goto 13\r
- 12 if (dabs(ss1) .lt. dabs(ss6)) goto 11\r
- 13 d = zero\r
- j = 3 * (j - 1)\r
- do 14 i = 1, 3\r
- k = ip(i + j)\r
- a(i,l) = ss(k)\r
- 14 d = d + (ss(k) * ss(k))\r
- if (d .gt. zero) d = one / dsqrt(d)\r
- do 15 i = 1, 3\r
- 15 a(i,l) = a(i,l) * d\r
- d = ((a(1,1) * a(1,3)) + (a(2,1) * a(2,3))) + (a(3,1) * a(3,3))\r
- m1 = 3\r
- m = 1\r
- if ((e1 - e2) .gt. (e2 - e3)) goto 16\r
- m1 = 1\r
- m = 3\r
- 16 p = zero\r
- do 17 i = 1, 3\r
- a(i,m1) = a(i,m1) - (d * a(i,m))\r
- 17 p = p + (a(i,m1) ** 2)\r
- if (p .le. tol) goto 19\r
- p = one / dsqrt(p)\r
- do 18 i = 1, 3\r
- 18 a(i,m1) = a(i,m1) * p\r
- goto 21\r
- 19 p = one\r
- do 20 i = 1, 3\r
- if (p .lt. dabs(a(i,m))) goto 20\r
- p = dabs(a(i,m))\r
- j = i\r
- 20 continue\r
- k = ip2312(j)\r
- l = ip2312(j + 1)\r
- p = dsqrt((a(k,m) ** 2) + (a(l,m) ** 2))\r
- if (p .le. tol) goto 40\r
- a(j,m1) = zero\r
- a(k,m1) = - (a(l,m) / p)\r
- a(l,m1) = a(k,m) / p\r
- 21 a(1,2) = (a(2,3) * a(3,1)) - (a(2,1) * a(3,3))\r
- a(2,2) = (a(3,3) * a(1,1)) - (a(3,1) * a(1,3))\r
-c****************** ROTATION MATRIX ************************************\r
-c 282 "rms.for"\r
- a(3,2) = (a(1,3) * a(2,1)) - (a(1,1) * a(2,3))\r
-c 284 "rms.for"\r
- 30 do 32 l = 1, 2\r
- d = zero\r
- do 31 i = 1, 3\r
- b(i,l) = ((r(i,1) * a(1,l)) + (r(i,2) * a(2,l))) + (r(i,3) * a(3,l\r
- &))\r
-c 288 "rms.for"\r
- 31 d = d + (b(i,l) ** 2)\r
- if (d .gt. zero) d = one / dsqrt(d)\r
- do 32 i = 1, 3\r
- 32 b(i,l) = b(i,l) * d\r
- d = ((b(1,1) * b(1,2)) + (b(2,1) * b(2,2))) + (b(3,1) * b(3,2))\r
- p = zero\r
- do 33 i = 1, 3\r
- b(i,2) = b(i,2) - (d * b(i,1))\r
- 33 p = p + (b(i,2) ** 2)\r
- if (p .le. tol) goto 35\r
- p = one / dsqrt(p)\r
- do 34 i = 1, 3\r
- 34 b(i,2) = b(i,2) * p\r
- goto 37\r
- 35 p = one\r
- do 36 i = 1, 3\r
- if (p .lt. dabs(b(i,1))) goto 36\r
- p = dabs(b(i,1))\r
- j = i\r
- 36 continue\r
- k = ip2312(j)\r
- l = ip2312(j + 1)\r
- p = dsqrt((b(k,1) ** 2) + (b(l,1) ** 2))\r
- if (p .le. tol) goto 40\r
- b(j,2) = zero\r
- b(k,2) = - (b(l,1) / p)\r
- b(l,2) = b(k,1) / p\r
- 37 b(1,3) = (b(2,1) * b(3,2)) - (b(2,2) * b(3,1))\r
- b(2,3) = (b(3,1) * b(1,2)) - (b(3,2) * b(1,1))\r
- b(3,3) = (b(1,1) * b(2,2)) - (b(1,2) * b(2,1))\r
- do 39 i = 1, 3\r
- do 39 j = 1, 3\r
-c****************** TRANSLATION VECTOR *********************************\r
-c 320 "rms.for"\r
- 39 u(i,j) = ((b(i,1) * a(j,1)) + (b(i,2) * a(j,2))) + (b(i,3) * a(j,3\r
- &))\r
- 40 do 41 i = 1, 3\r
-c********************** RMS ERROR **************************************\r
-c 323 "rms.for"\r
- 41 t(i) = ((yc(i) - (u(i,1) * xc(1))) - (u(i,2) * xc(2))) - (u(i,3)\r
- & * xc(3))\r
- 50 do 51 i = 1, 3\r
- if (e(i) .lt. zero) e(i) = zero\r
- 51 e(i) = dsqrt(e(i))\r
- ier = 0\r
- if (e2 .le. (e1 * 1.0d-05)) ier = -1\r
- d = e3\r
- if (sigma .ge. 0.0) goto 52\r
- d = - d\r
- if ((e2 - e3) .le. (e1 * 1.0d-05)) ier = -1\r
- 52 d = (d + e2) + e1\r
- rms = (e0 - d) - d\r
- if (rms .lt. 0.0) rms = 0.0\r
- return \r
-c.....END U3B...........................................................\r
-c----------------------------------------------------------\r
-c THE END\r
-c----------------------------------------------------------\r
-c 338 "rms.for"\r
- end\r
-\r
git://source.jalview.org / jabaws.git / blobdiff