diff --git a/modules/bindings/pymod/tmtools.py b/modules/bindings/pymod/tmtools.py
index 522a2dd559e31fb6656dca3ad20dc9a42d360722..68eacab99ea6c1dd03ab854a2570e3de20528a32 100644
--- a/modules/bindings/pymod/tmtools.py
+++ b/modules/bindings/pymod/tmtools.py
@@ -63,8 +63,8 @@ def _ParseTmAlign(lines):
tf2[4], tf3[2], tf3[3], tf3[4])
tf=geom.Mat4(rot)
tf.PasteTranslation(geom.Vec3(tf1[1], tf2[1], tf3[1]))
- seq1 = seq.CreateSequence("1",lines[20].strip())
- seq2 = seq.CreateSequence("2",lines[22].strip())
+ seq1 = seq.CreateSequence("1",lines[26].strip())
+ seq2 = seq.CreateSequence("2",lines[28].strip())
alignment = seq.CreateAlignment()
alignment.AddSequence(seq2)
alignment.AddSequence(seq1)
diff --git a/modules/bindings/src/tmalign.f b/modules/bindings/src/tmalign.f
index 73bbd905d13d4196da88392504a7475abc1fc3a4..152453538d37421c7a836f56142a9d980be4028b 100644
--- a/modules/bindings/src/tmalign.f
+++ b/modules/bindings/src/tmalign.f
@@ -1,39 +1,50 @@
-*************************************************************************
+**************************************************************************
* This program is to identify the best alignment of two protein
-* structures to give the best TM-score. By default, TM-score is
-* normalized by the second protein. The program can be freely
-* copied or modified or redistributed.
-* (For comments, please email to: zhng@umich.edu)
+* structures that gives the highest TM-score. Input structures must
+* be in the PDB format. By default, TM-score is normalized by the
+* second protein. Users can obtain a brief instruction by simply
+* running the program without arguments. For comments/suggestions,
+* please contact email: zhng@umich.edu.
*
-* Reference:
+* Reference to cite:
* Yang Zhang, Jeffrey Skolnick, Nucl. Acid Res. 2005 33: 2303-9
-*
-************************ updating history *******************************
+*
+* Permission to use, copy, modify, and distribute this program for
+* any purpose, with or without fee, is hereby granted, provided that
+* the notices on the head, the reference information, and this
+* copyright notice appear in all copies or substantial portions of
+* the Software. It is provided "as is" without express or implied
+* warranty.
+************************ updating history ********************************
* 2005/06/01: A small bug of two-point superposition was fixed.
* 2005/10/19: the program was reformed so that the alignment
* results are not dependent on the specific compilers.
* 2006/06/20: select 'A' if there is altLoc when reading PDB file.
-* 2007/02/27: rotation matrix from Chain-1 to Chain-2 is added.
-* 2007/04/18: add options with TM-score normalized by average
+* 2007/02/27: rotation matrix from Chain-1 to Chain-2 was added.
+* 2007/04/18: added options with TM-score normalized by average
* length, shorter length, or longer length of two
* structures.
-* 2007/05/23: add additional output file 'TM.sup_all' for showing
+* 2007/05/23: added additional output file 'TM.sup_all' for showing
* all atoms while 'TM.sup' is only for aligned atoms
-* 2007/09/19: add a new feature alignment to deal with the problem
+* 2007/09/19: added a new feature alignment to deal with the problem
* of aligning fractional structures (e.g. protein
* interfaces).
* 2007/10/16: A bug for irregular bond-length models was fixed.
-* 2009/03/14: A new initial alignment is added and previous initial
-* alignments are further enhanced. This change increases
+* 2009/03/14: A new initial alignment was added and previous initial
+* alignments are further enhanced. This change increased
* accuracy by 4% but increases CPU cost by 40%.
* 2009/08/20: A bug for asymmetry alignment result was fixed.
-*
-*************************************************************************
+* 2010/08/02: A new RMSD matrix was used to remove obsolete statements.
+* Staled subroutines were deleted.
+* 2011/01/03: The length of pdb file names were extended to 500
+* 2011/01/24: Fixed a bug on output file name created on 2011/01/03
+* 2011/01/30: An open source license is attached to the program
+**************************************************************************
- program compares
+ program TMalign
PARAMETER(nmax=5000)
PARAMETER(nmax2=10000)
-
+
COMMON/BACKBONE/XA(3,nmax,0:1)
common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
common/alignrst/invmap0(nmax)
@@ -41,28 +52,26 @@
common/d0/d0,anseq
common/d0min/d0_min
common/d00/d00,d002
-
- character*100 fnam,pdb(100),outname
+
+ character*500 fnam,pdb(100),outname
character*3 aa(-1:20),aanam,ss1(nmax),ss2(nmax)
- character*100 s,du
- character*200 outnameall_tmp,outnameall
+ character*500 s,du
+ character*504 outnameall_tmp,outnameall
character seq1(0:nmax),seq2(0:nmax)
character aseq1(nmax2),aseq2(nmax2),aseq3(nmax2)
-
- dimension xx(nmax),yy(nmax),zz(nmax)
+
dimension m1(nmax),m2(nmax)
dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
common/init/invmap_i(nmax)
common/TM/TM,TMmax
- common/n1n2/n1(nmax),n2(nmax)
common/d8/d8
common/initial4/mm1(nmax),mm2(nmax)
ccc RMSD:
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms !armsd is real
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms !armsd is real
data w /nmax*1.0/
ccc
@@ -77,7 +86,7 @@ ccc
& 'D','N','L','K','E',
& 'Q','R','H','F','Y',
& 'W','C'/
-
+
call getarg(1,fnam)
if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then
write(*,*)
@@ -93,14 +102,14 @@ ccc
write(*,*)'2. Run TM-align and output the superposition ',
& 'to ''TM.sup'' and ''TM.sup_all'':'
write(*,*)' >TMalign structure.pdb target.pdb -o TM.sup'
- write(*,*)' To view the superimposed structures of the',
+ write(*,*)' To view the superimposed structures of the',
& ' aligned regions by rasmol:'
- write(*,*)' >rasmol -script TM.sup)'
- write(*,*)' To view the superimposed structures of all',
+ write(*,*)' >rasmol -script TM.sup'
+ write(*,*)' To view the superimposed structures of all',
& ' regions by rasmol:'
- write(*,*)' >rasmol -script TM.sup_all)'
+ write(*,*)' >rasmol -script TM.sup_all'
write(*,*)
- write(*,*)'3. If you want TM-score normalized by ',
+ write(*,*)'3. If you want output TM-score normalized by ',
& 'an assigned length, e.g. 100 aa:'
write(*,*)' >TMalign structure.pdb target.pdb -L 100'
write(*,*)' If you want TM-score normalized by the ',
@@ -113,18 +122,12 @@ ccc
& 'longer length of two structures:'
write(*,*)' >TMalign structure.pdb target.pdb -c'
write(*,*)
-c write(*,*)'5. If you want to set a minimum cutoff for d0, ',
-c & 'e.g. d0>3.0'
-c write(*,*)' (By default d0>0.5, this option need ',
-c & 'be considered only when L<35 aa)'
-c write(*,*)' >TMalign structure.pdb target.pdb -dmin 3.0'
-c write(*,*)
write(*,*)'(All above options do not change the ',
- & 'final structure alignment result)'
+ & 'final structure alignment results)'
write(*,*)
goto 9999
endif
-
+
******* options ----------->
m_out=-1 !decided output
m_fix=-1 !fixed length-scale only for output
@@ -156,7 +159,7 @@ c write(*,*)
i=i+1
call getarg(i,fnam)
read(fnam,*)d0_fix
- elseif(fnam.eq.'-a')then ! this will change the superposed output but not the alignment
+ elseif(fnam.eq.'-a')then !change superposed output but not the alignment
m_ave=1
i=i+1
elseif(fnam.eq.'-b')then
@@ -236,8 +239,8 @@ ccccccccc read data from the second CA file:
close(10)
nseq2=i
c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-*!!! Scale of TM-score in search is based on the smaller protein --------->
+
+* Scale of TM-score in search is based on the smaller protein --------->
d0_min=0.5
if(m_d0_min.eq.1)then
d0_min=d0_min_input !for search
@@ -257,13 +260,9 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if(d00.lt.4.5)d00=4.5
d002=d00**2
nseq=max(nseq1,nseq2)
- do i=1,nseq
- n1(i)=i
- n2(i)=i
- enddo
***** do alignment **************************
- CALL super_align !to find invmap(j)
+ CALL TM_align !to find invmap(j)
************************************************************
*** resuperpose to find residues of dis<d8 ------------------------>
@@ -283,10 +282,10 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
endif
enddo
d0_input=d0
- call TMscore8(d0_input,n_al,xtm1,ytm1,ztm1,n1,n_al,
- & xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !TM-score with dis<d8 only
+ call TMscore8(d0_input,n_al,xtm1,ytm1,ztm1,n_al,
+ & xtm2,ytm2,ztm2,TM,Rcomm,Lcomm) !TM-score with dis<d8 only
-*!!! Output TM-score is based on the second protein------------------>
+* Output TM-score is based on the second protein------------------>
d0_min=0.5 !for output
anseq=nseq2 !length for defining final TMscore
if(m_ave.eq.1)anseq=(nseq1+nseq2)/2.0 !<L>
@@ -326,8 +325,8 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
seq_id=float(n_eq)/(n_al+0.00000001)
n8_al=j
d0_input=d0
- call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n1,n8_al,
- & xtm2,ytm2,ztm2,n2,TM8,Rcomm,Lcomm) !normal TMscore
+ call TMscore(d0_input,n8_al,xtm1,ytm1,ztm1,n8_al,
+ & xtm2,ytm2,ztm2,TM8,Rcomm,Lcomm) !normal TMscore
rmsd8_al=Rcomm
TM8=TM8*n8_al/anseq !TM-score after cutoff
@@ -343,8 +342,8 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
& 'M-score *'
write(*,*)'* Reference: Y. Zhang and J. Skolnick, Nucl. Acids Re',
& 's. 2005 33, 2302-9 *'
- write(*,*)'* Comments on the program, please email to: yzhang@ku',
- & '.edu *'
+ write(*,*)'* Comments on the program, please email to: zhng@umic',
+ & 'h.edu *'
write(*,*)'*****************************************************',
& '*********************'
write(*,*)
@@ -373,19 +372,19 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
enddo
if(L.gt.3)then
call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2
- armsd=dsqrt(rms/L)
- write(*,*)'-------- rotation matrix to rotate Chain-1 to ',
+ write(*,*)'-------- Rotation matrix to rotate Chain-1 to ',
& 'Chain-2 ------'
- write(*,*)'i t(i) u(i,1) u(i,2) ',
- & ' u(i,3)'
+ write(*,*)'m t(m) u(m,1) u(m,2) ',
+ & ' u(m,3)'
do i=1,3
write(*,204)i,t(i),u(i,1),u(i,2),u(i,3)
enddo
-c do i=1,nseq1
-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)
-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)
-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)
-c enddo
+ write(*,*)'Code for rotating Chain-1 from (x,y,z) to (X,Y,Z):'
+ write(*,*)' do i=1,L'
+ write(*,*)' X(i)=t(1)+u(1,1)*x(i)+u(1,2)*y(i)+u(1,3)*z(i)'
+ write(*,*)' Y(i)=t(2)+u(2,1)*x(i)+u(2,2)*y(i)+u(2,3)*z(i)'
+ write(*,*)' Z(i)=t(3)+u(3,1)*x(i)+u(3,2)*y(i)+u(3,3)*z(i)'
+ write(*,*)' enddo'
write(*,*)
endif
204 format(I2,f18.10,f15.10,f15.10,f15.10)
@@ -449,7 +448,7 @@ ccc
k=0
outnameall_tmp=outname//'_all'
outnameall=''
- do i=1,200
+ do i=1,504
if(outnameall_tmp(i:i).ne.' ')then
k=k+1
outnameall(k:k)=outnameall_tmp(i:i)
@@ -546,8 +545,8 @@ ccc
aseq3(ii)=' '
enddo
write(*,50)
- 50 format('(":" denotes the residue pairs of distance < 5.0 ',
- & 'Angstrom)')
+ 50 format('(":" denotes residue pairs of d < 5 Angstrom,',
+ & ' "." denotes other aligned residues)')
write(*,10)(aseq1(i),i=1,ii)
write(*,10)(aseq3(i),i=1,ii)
write(*,10)(aseq2(i),i=1,ii)
@@ -564,7 +563,7 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
***********************************************************************
***********************************************************************
***********************************************************************
- SUBROUTINE super_align
+ SUBROUTINE TM_align
PARAMETER(nmax=5000)
COMMON/BACKBONE/XA(3,nmax,0:1)
common/length/nseq1,nseq2
@@ -573,19 +572,13 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
common/zscore/zrms,n_al,rmsd_al
common/TM/TM,TMmax
common/init/invmap_i(nmax)
- common/init1/invmap_ii(nmax)
dimension gapp(100)
-
+
TMmax=0
n_gapp=2
gapp(1)=-0.6
gapp(2)=0
-
-c n_gapp=11
-c do i=1,n_gapp
-c gapp(i)=-(n_gapp-i)
-c enddo
-
+
*11111111111111111111111111111111111111111111111111111111
* get initial alignment from gapless threading
**********************************************************
@@ -594,65 +587,25 @@ c enddo
invmap(i)=invmap_i(i) !with highest zcore
enddo
call get_score !TM, matrix score(i,j)
-
+
if(TM.gt.TMmax)then
TMmax=TM
do j=1,nseq2
invmap0(j)=invmap(j)
enddo
endif
-
+
*****************************************************************
* initerative alignment, for different gap_open:
*****************************************************************
- DO 111 i_gapp=1,n_gapp !different gap panalties
+ DO 1 i_gapp=1,n_gapp !different gap panalties
GAP_OPEN=gapp(i_gapp) !gap panalty
- do 222 id=1,30 !maximum interation is 200
+ do 11 id=1,30 !maximum interation is 200
call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
* Input: score(i,j), and gap_open
* Output: invmap(j)
-
- call get_score !calculate TM-score, score(i,j)
-c record the best alignment in whole search ---------->
- if(TM.gt.TMmax)then
- TMmax=TM
- do j=1,nseq2
- invmap0(j)=invmap(j)
- enddo
- endif
- if(id.gt.1)then
- diff=abs(TM-TM_old)
- if(diff.lt.0.000001)goto 33
- endif
- TM_old=TM
- 222 continue
- 33 continue
- 111 continue
-
- do i=1,nseq2
- invmap(i)=invmap_ii(i) !with highest zcore
- enddo
- call get_score !TM, matrix score(i,j)
-
- if(TM.gt.TMmax)then
- TMmax=TM
- do j=1,nseq2
- invmap0(j)=invmap(j)
- enddo
- endif
- if(TM.le.TMmax*0.6)goto 51
-****************************************************************
-* initerative alignment, for different gap_open:
-*****************************************************************
- DO 112 i_gapp=1,n_gapp !different gap panalties
- GAP_OPEN=gapp(i_gapp) !gap panalty
- do 223 id=1,2 !maximum interation is 200
- call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
-* Input: score(i,j), and gap_open
-* Output: invmap(j)
-
+
call get_score !calculate TM-score, score(i,j)
- !print *,"In1a",TM,TMmax
c record the best alignment in whole search ---------->
if(TM.gt.TMmax)then
TMmax=TM
@@ -662,14 +615,13 @@ c record the best alignment in whole search ---------->
endif
if(id.gt.1)then
diff=abs(TM-TM_old)
- if(diff.lt.0.000001)goto 34
+ if(diff.lt.0.000001)goto 111
endif
TM_old=TM
- 223 continue
- 34 continue
- 112 continue
- 51 continue
-
+ 11 continue
+ 111 continue
+ 1 continue
+
*222222222222222222222222222222222222222222222222222222222
* get initial alignment from secondary structure alignment
**********************************************************
@@ -684,17 +636,17 @@ c record the best alignment in whole search ---------->
invmap0(j)=invmap(j)
enddo
endif
- if(TM.le.TMmax*0.2)goto 52
+ if(TM.le.TMmax*0.2)goto 2222
*****************************************************************
* initerative alignment, for different gap_open:
*****************************************************************
- DO 1111 i_gapp=1,n_gapp !different gap panalties
+ DO 2 i_gapp=1,n_gapp !different gap panalties
GAP_OPEN=gapp(i_gapp) !gap panalty
- do 2222 id=1,30 !maximum interation is 200
+ do 22 id=1,30 !maximum interation is 200
call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
* Input: score(i,j), and gap_open
* Output: invmap(j)
-
+
call get_score !calculate TM-score, score(i,j)
c write(*,21)gap_open,rmsd_al,n_al,TM
c record the best alignment in whole search ---------->
@@ -706,14 +658,14 @@ c record the best alignment in whole search ---------->
endif
if(id.gt.1)then
diff=abs(TM-TM_old)
- if(diff.lt.0.000001)goto 333
+ if(diff.lt.0.000001)goto 222
endif
TM_old=TM
- 2222 continue
- 333 continue
- 1111 continue
- 52 continue
-
+ 22 continue
+ 222 continue
+ 2 continue
+ 2222 continue
+
*555555555555555555555555555555555555555555555555555555555555555555
* get initial alignment of local structure superposition
*******************************************************************
@@ -722,43 +674,42 @@ c record the best alignment in whole search ---------->
invmap(i)=invmap_i(i) !with highest zcore
enddo
call get_score !TM, matrix score(i,j)
-
+
if(TM.gt.TMmax)then
TMmax=TM
do j=1,nseq2
invmap0(j)=invmap(j)
- enddo
+ enddo
endif
- if(TM.le.TMmax*0.8)goto 53
+ if(TM.le.TMmax*0.8)goto 5555
*****************************************************************
* initerative alignment, for different gap_open:
*****************************************************************
- DO 88 i_gapp=1,n_gapp !different gap panalties
- GAP_OPEN=gapp(i_gapp) !gap panalty
- do 888 id=1,2 !maximum interation is 200
+ DO 5 i_gapp=1,n_gapp !different gap panalties
+ GAP_OPEN=gapp(i_gapp) !gap panalty
+ do 55 id=1,2 !maximum interation is 200
call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
* Input: score(i,j), and gap_open
* Output: invmap(j)
call get_score !calculate TM-score, score(i,j)
-c print *,"Ini5a",TM,TMmax
c record the best alignment in whole search ---------->
if(TM.gt.TMmax)then
TMmax=TM
do j=1,nseq2
invmap0(j)=invmap(j)
- enddo
+ enddo
endif
if(id.gt.1)then
diff=abs(TM-TM_old)
- if(diff.lt.0.000001)goto 880
+ if(diff.lt.0.000001)goto 555
endif
- TM_old=TM
- 888 continue
- 880 continue
- 88 continue
- 53 continue
-
+ TM_old=TM
+ 55 continue
+ 555 continue
+ 5 continue
+ 5555 continue
+
*333333333333333333333333333333333333333333333333333333333333
* get initial alignment from invmap0+SS
*************************************************************
@@ -773,13 +724,13 @@ c record the best alignment in whole search ---------->
invmap0(j)=invmap(j)
enddo
endif
- if(TM.le.TMmax*0.2)goto 54
+ if(TM.le.TMmax*0.2)goto 3333
*****************************************************************
* initerative alignment, for different gap_open:
*****************************************************************
- DO 1110 i_gapp=1,n_gapp !different gap panalties
+ DO 3 i_gapp=1,n_gapp !different gap panalties
GAP_OPEN=gapp(i_gapp) !gap panalty
- do 2220 id=1,30 !maximum interation is 200
+ do 33 id=1,30 !maximum interation is 200
call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
* Input: score(i,j), and gap_open
* Output: invmap(j)
@@ -795,13 +746,13 @@ c record the best alignment in whole search ---------->
endif
if(id.gt.1)then
diff=abs(TM-TM_old)
- if(diff.lt.0.000001)goto 330
+ if(diff.lt.0.000001)goto 333
endif
TM_old=TM
- 2220 continue
- 330 continue
- 1110 continue
- 54 continue
+ 33 continue
+ 333 continue
+ 3 continue
+ 3333 continue
*444444444444444444444444444444444444444444444444444444444
* get initial alignment of pieces from gapless threading
@@ -810,12 +761,6 @@ c record the best alignment in whole search ---------->
do i=1,nseq2
invmap(i)=invmap_i(i) !with highest zcore
enddo
- do i=1,nseq2
- do j=1,nseq1
- score(i,j)=0
- enddo
- enddo
-
call get_score !TM, matrix score(i,j)
if(TM.gt.TMmax)then
TMmax=TM
@@ -823,19 +768,18 @@ c record the best alignment in whole search ---------->
invmap0(j)=invmap(j)
enddo
endif
- if(TM.le.TMmax*0.3)goto 55
+ if(TM.le.TMmax*0.3)goto 4444
*****************************************************************
* initerative alignment, for different gap_open:
*****************************************************************
- DO 44 i_gapp=2,n_gapp !different gap panalties
+ DO 4 i_gapp=2,n_gapp !different gap panalties
GAP_OPEN=gapp(i_gapp) !gap panalty
- do 444 id=1,2 !maximum interation is 200
+ do 44 id=1,2 !maximum interation is 200
call DP(NSEQ1,NSEQ2) !produce alignment invmap(j)
* Input: score(i,j), and gap_open
* Output: invmap(j)
call get_score !calculate TM-score, score(i,j)
- ! print *,"Ini4a",TM,TMmax
c record the best alignment in whole search ---------->
if(TM.gt.TMmax)then
TMmax=TM
@@ -843,9 +787,9 @@ c record the best alignment in whole search ---------->
invmap0(j)=invmap(j)
enddo
endif
- 444 continue
- 44 continue
- 55 continue
+ 44 continue
+ 4 continue
+ 4444 continue
c^^^^^^^^^^^^^^^ best alignment invmap0(j) found ^^^^^^^^^^^^^^^^^^
RETURN
@@ -863,45 +807,20 @@ c^^^^^^^^^^^^^^^ best alignment invmap0(j) found ^^^^^^^^^^^^^^^^^^
common/zscore/zrms,n_al,rmsd_al
common/TM/TM,TMmax
common/init/invmap_i(nmax)
- common/init1/invmap_ii(nmax)
- common/d0/d0,anseq
- common/d0min/d0_min
- common/speed/align_len1
- real n_cut
- integer n_all
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms
- data w /nmax*1.0/
- dimension xtmm1(nmax),ytmm1(nmax),ztmm1(nmax)
- dimension xtmm2(nmax),ytmm2(nmax),ztmm2(nmax)
-
- d01=d0+1.5
- if(d01.lt.d0_min)d01=d0_min
- d02=d01*d01
- n_jump=1
- n_cut=0.99
aL=min(nseq1,nseq2)
- idel=aL/2.0 !minimum size of considered fragment
+ idel=aL/2.0 !minimum size of considered fragment
if(idel.le.5)idel=5
n1=-nseq2+idel
n2=nseq1-idel
GL_max=0
- GL_maxA=0
- !GL_max_cut=0.99
- do ishift=n1,n2,n_jump
+ do ishift=n1,n2
L=0
do j=1,nseq2
i=j+ishift
if(i.ge.1.and.i.le.nseq1)then
L=L+1
invmap(j)=i
- xtmm1(L)=xa(1,i,0)
- ytmm1(L)=xa(1,i,0)
- ztmm1(L)=xa(1,i,0)
- xtmm2(L)=xa(1,j,1)
- ytmm2(L)=xa(1,j,1)
- ztmm2(L)=xa(1,j,1)
else
invmap(j)=-1
endif
@@ -909,57 +828,14 @@ c^^^^^^^^^^^^^^^ best alignment invmap0(j) found ^^^^^^^^^^^^^^^^^^
if(L.ge.idel)then
call get_GL(GL)
if(GL.gt.GL_max)then
-
GL_max=GL
do i=1,nseq2
invmap_i(i)=invmap(i)
enddo
endif
- if(GL.gt.GL_max)then
- !if(nseq1.le.nseq2)then
- do ii=1,L
- r_1(1,ii)=xtmm1(ii)
- r_1(2,ii)=ytmm1(ii)
- r_1(3,ii)=ztmm1(ii)
- r_2(1,ii)=xtmm2(ii)
- r_2(2,ii)=ytmm2(ii)
- r_2(3,ii)=ztmm2(ii)
- enddo
-
- call u3b(w,r_1,r_2,L,1,rms,u,t,ier) !u rotate r_1 to r_2
-
- do i1=1,nseq2
- xx=t(1)+u(1,1)*xa(1,i1,1)+u(1,2)*xa(2,i1,1)+u(1,3)*
- & xa(3,i1,1)
- yy=t(2)+u(2,1)*xa(1,i1,1)+u(2,2)*xa(2,i1,1)+u(2,3)*
- & xa(3,i1,1)
- zz=t(3)+u(3,1)*xa(1,i1,1)+u(3,2)*xa(2,i1,1)+u(3,3)*
- & xa(3,i1,1)
- do j1=1,nseq1
- dd=(xx-xa(1,j1,0))**2+(yy-xa(2,j1,0))**2+(zz-xa
- & (3,j1,0))**2
- score(i1,j1)=1/(1+dd/d0)
- enddo
- enddo
-
-*********extract alignement with score(i,j) *****************
- align_len1=0
- call DP(NSEQ1,NSEQ2)
- if(align_len1.gt.L*n_cut.and.align_len1.gt.idel)then
- call get_GL(GL)
- !print *,GL
- if(GL.gt.GL_max.and.GL.gt.GL_maxA)then
- GL_maxA=GL
- do i1=1,aLL
- invmap_ii(i1)=invmap(i1)
- !print *,i,invmap_ii(i)
- enddo
- endif
- endif
- endif
endif
enddo
-
+
return
end
@@ -1093,15 +969,14 @@ c 1->coil, 2->helix, 3->turn, 4->strand
common/initial4/mm1(nmax),mm2(nmax)
logical contin
- dimension ifr2(2,nmax,nmax),Lfr2(2,nmax),Lfr_max2(2),i_fr2(2)
+ dimension ifr2(2,nmax,nmax),Lfr2(2,nmax),i_fr2(2)
dimension ifr(nmax)
dimension mm(2,nmax)
fra_min=4 !>=4,minimum fragment for search
fra_min1=fra_min-1 !cutoff for shift, save time
dcu0=3.85
-c dcu_min=3.65
-
+
ccc Find the smallest continuous fragments -------->
do i=1,nseq1
mm(1,i)=mm1(i)
@@ -1128,15 +1003,11 @@ ccc Find the smallest continuous fragments -------->
contin=.false.
if(dcu.gt.dcu0)then
if(dis.lt.dcu)then
-c if(dis.gt.dcu_min)then
- contin=.true.
-c endif
+ contin=.true.
endif
elseif(mm(k,i).eq.(mm(k,i-1)+1))then
if(dis.lt.dcu)then
-c if(dis.gt.dcu_min)then
- contin=.true.
-c endif
+ contin=.true.
endif
endif
if(contin)then
@@ -1272,7 +1143,7 @@ ccc get initial ------------->
* position. *
**************************************************************
subroutine get_initial5
- PARAMETER (nmax=5000)
+ PARAMETER(nmax=5000)
COMMON/BACKBONE/XA(3,nmax,0:1)
common/length/nseq1,nseq2
common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
@@ -1283,138 +1154,79 @@ ccc get initial ------------->
common/d0min/d0_min
common/init/invmap_i(nmax)
common/initial4/mm1(nmax),mm2(nmax)
- common/tinit/invmap_ii(nmax)
common/sec/isec(nmax),jsec(nmax)
- logical contin
- integer n_all,secmatch,coilcnt
- real n_max_cut,n_all_max,n_cutt
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms
data w /nmax*1.0/
common/inv/invmap_a(nmax)
- common/speed/align_len1
- integer aL,aLL,m1,m2,n_frag,n_frag1,n_frag2,rem
+ integer aL,m1,m2,n_frag
+
***** setting parameters ************************************
- m1=0
- m2=0
- n_frag=20
- if(nseq1.gt.250.and.nseq2.gt.250)then
+ d01=d0+1.5
+ if(d01.lt.d0_min)d01=d0_min
+ d02=d01*d01
+ GLmaxA=0
+
+ aL=min(nseq1,nseq2)
+ if(aL.gt.250)then
n_frag=50
- goto 100
- elseif(nseq1.gt.200.and.nseq2.gt.200)then
+ elseif(aL.gt.200)then
n_frag=40
- goto 100
- elseif(nseq1.gt.150.and.nseq2.gt.150)then
+ elseif(aL.gt.150)then
n_frag=30
- goto 100
- !elseif(nseq1.lt.100.and.nseq2.lt.100)then
- ! n_frag=15
- ! goto 100
- !elseif(nseq1.lt.200.and.nseq2.gt.200)then
- ! n_frag=22
- ! goto 100
- !elseif(nseq1.gt.150.and.nseq2.gt.150)then
- ! n_frag=20
- ! goto 100
- !elseif(nseq1.gt.150.and.nseq2.lt.150)then
- ! n_frag=18
- ! goto 100
- !elseif(nseq1.lt.150.and.nseq2.gt.150)then
- ! n_frag=18
- ! goto 100
- !
- !elseif(nseq1.gt.100.and.nseq2.gt.100)then
- ! n_frag=16
- ! goto 100
- !elseif(nseq1.gt.100.and.nseq2.lt.100)then
- ! n_frag=14
- ! goto 100
- !elseif(nseq1.lt.100.and.nseq2.lt.100)then
- ! n_frag=14
- ! goto 100
- !else
- ! n_frag=12
+ else
+ n_frag=20 !length of fragment for superposition
endif
+ if(n_frag.gt.aL/2)n_frag=aL/2
+ ns=20 !tail length to discard
+ if(ns.gt.aL/2)ns=aL/2
- 100 continue
-
-
- n_jump=n_frag
- !n_all=0
- !n_all_max=0
- !n_max_cut=0.7
- !n_cutt=0
- d01=d0+1.5
- if(d01.lt.d0_min)d01=d0_min
- d02=d01*d01
- !GLmaxB=0
- GLmaxA=0
- GL=0
- !GLmaxA_cut=0.12
- !n_frag1=n_frag-1
- !n_frag2=n_frag/2
+ m1=nseq1-n_frag-ns
+ m2=nseq2-n_frag-ns
+ do ii=ns,m1,n_frag
+ do jj=ns,m2,n_frag
+ do k=1,n_frag
+ iii=ii+k-1
+ jjj=jj+k-1
+ r_1(1,k)=xa(1,iii,0)
+ r_1(2,k)=xa(2,iii,0)
+ r_1(3,k)=xa(3,iii,0)
+ r_2(1,k)=xa(1,jjj,1)
+ r_2(2,k)=xa(2,jjj,1)
+ r_2(3,k)=xa(3,jjj,1)
+ enddo
-
- m1=nseq1-n_frag-20
- m2=nseq2-n_frag-20
-
- do ii=20,m1,n_jump
- do jj=20,m2,n_jump
-
- do k=1,n_frag
- iii=ii+k-1
- jjj=jj+k-1
- r_1(1,k)=xa(1,iii,0)
- r_1(2,k)=xa(2,iii,0)
- r_1(3,k)=xa(3,iii,0)
- r_2(1,k)=xa(1,jjj,1)
- r_2(2,k)=xa(2,jjj,1)
- r_2(3,k)=xa(3,jjj,1)
- enddo
-
*********superpose the two structures and rotate it *****************
- call u3b(w,r_1,r_2,n_frag,1,rms,u,t,ier) !u rotate r_1 to r_2
-
-
- do i1=1,nseq1
- xx=t(1)+u(1,1)*xa(1,i1,0)+u(1,2)*xa(2,i1,0)+u(1,3)
- & *xa(3,i1,0)
- yy=t(2)+u(2,1)*xa(1,i1,0)+u(2,2)*xa(2,i1,0)+u(2,3)
- & *xa(3,i1,0)
- zz=t(3)+u(3,1)*xa(1,i1,0)+u(3,2)*xa(2,i1,0)+u(3,3)
- & *xa(3,i1,0)
- do j1=1,nseq2
- dd=(xx-xa(1,j1,1))**2+(yy-xa(2,j1,1))**2+
- & (zz-xa(3,j1,1))**2
- score(i1,j1)=1/(1+dd/d02) ! changing
- enddo
+ call u3b(w,r_1,r_2,n_frag,1,rms,u,t,ier) !u rotate r_1 to r_2
+ do i1=1,nseq1
+ xx=t(1)+u(1,1)*xa(1,i1,0)+u(1,2)*xa(2,i1,0)+u(1,3)
+ & *xa(3,i1,0)
+ yy=t(2)+u(2,1)*xa(1,i1,0)+u(2,2)*xa(2,i1,0)+u(2,3)
+ & *xa(3,i1,0)
+ zz=t(3)+u(3,1)*xa(1,i1,0)+u(3,2)*xa(2,i1,0)+u(3,3)
+ & *xa(3,i1,0)
+ do j1=1,nseq2
+ dd=(xx-xa(1,j1,1))**2+(yy-xa(2,j1,1))**2+
+ & (zz-xa(3,j1,1))**2
+ score(i1,j1)=1/(1+dd/d02) ! changing
enddo
-
-
+ enddo
+
*********extract alignement with score(i,j) *****************
- !align_len1=0
- call DP(NSEQ1,NSEQ2)
-
-
- call get_GL(GL)
-
- if(GL.gt.GLmaxA)then
- GLmaxA=GL
- do j1=1,nseq2
- invmap_i(j1)=invmap(j1)
- enddo
- endif
-
+ call DP(NSEQ1,NSEQ2)
+ call get_GL(GL)
+ if(GL.gt.GLmaxA)then
+ GLmaxA=GL
+ do j1=1,nseq2
+ invmap_i(j1)=invmap(j1)
+ enddo
+ endif
enddo
-
enddo
-
+
return
end
-
-
-
**************************************************************
* smooth the secondary structure assignment
**************************************************************
@@ -1623,10 +1435,7 @@ ccc get initial ------------->
common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
common/zscore/zrms,n_al,rmsd_al
common/d0/d0,anseq
- dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
- dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
common/TM/TM,TMmax
- common/n1n2/n1(nmax),n2(nmax)
common/d00/d00,d002
dimension xo1(nmax),yo1(nmax),zo1(nmax)
@@ -1634,8 +1443,8 @@ ccc get initial ------------->
dimension dis2(nmax)
ccc RMSD:
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms !armsd is real
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms !armsd is real
data w /nmax*1.0/
ccc
@@ -1744,11 +1553,10 @@ c^^^^^^^^^^^^^^^^ GL done ^^^^^^^^^^^^^^^^^^^^^^^^^^^
dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
common/TM/TM,TMmax
- common/n1n2/n1(nmax),n2(nmax)
-
+
ccc RMSD:
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms !armsd is real
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms !armsd is real
data w /nmax*1.0/
ccc
@@ -1773,8 +1581,8 @@ ccc for rotation matrix:
enddo
*** calculate TM-score for the given alignment----------->
d0_input=d0
- call TMscore8_search(d0_input,n_al,xtm1,ytm1,ztm1,n1,
- & n_al,xtm2,ytm2,ztm2,n2,TM,Rcomm,Lcomm) !simplified search engine
+ call TMscore8_search(d0_input,n_al,xtm1,ytm1,ztm1,
+ & n_al,xtm2,ytm2,ztm2,TM,Rcomm,Lcomm) !simplified search engine
TM=TM*n_al/anseq !TM-score
*** calculate score matrix score(i,j)------------------>
do i=1,n_al
@@ -1808,14 +1616,11 @@ c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^
common/zscore/zrms,n_al,rmsd_al
common/d0/d0,anseq
common/d0min/d0_min
- dimension xtm1(nmax),ytm1(nmax),ztm1(nmax)
- dimension xtm2(nmax),ytm2(nmax),ztm2(nmax)
common/TM/TM,TMmax
- common/n1n2/n1(nmax),n2(nmax)
ccc RMSD:
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms !armsd is real
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms !armsd is real
data w /nmax*1.0/
ccc
@@ -1861,10 +1666,8 @@ c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^
*
* L1--Length of the first structure
* (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
-* n1(i)--Residue sequence number of i'th residue at the first structure
* L2--Length of the second structure
* (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
-* n2(i)--Residue sequence number of i'th residue at the second structure
* TM--TM-score of the comparison
* Rcomm--RMSD of two structures in the common aligned residues
* Lcomm--Length of the common aligned regions
@@ -1877,28 +1680,28 @@ c^^^^^^^^^^^^^^^^ score(i,j) done ^^^^^^^^^^^^^^^^^^^^^^^^^^^
*************************************************************************
*************************************************************************
*** dis<8, simplified search engine
- subroutine TMscore8_search(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,
+ subroutine TMscore8_search(dx,L1,x1,y1,z1,L2,x2,y2,z2,
& TM,Rcomm,Lcomm)
PARAMETER(nmax=5000)
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
+ common/nres/nseqA,nseqB
common/para/d,d0
common/d0min/d0_min
common/align/n_ali,iA(nmax),iB(nmax)
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
dimension k_ali(nmax),k_ali0(nmax)
- dimension L_ini(100),iq(nmax)
+ dimension L_ini(100)
common/scores/score
double precision score,score_max
dimension xa(nmax),ya(nmax),za(nmax)
dimension iL0(nmax)
- dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
- dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
+ dimension x1(nmax),y1(nmax),z1(nmax)
+ dimension x2(nmax),y2(nmax),z2(nmax)
ccc RMSD:
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms !armsd is real
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms !armsd is real
data w /nmax*1.0/
ccc
@@ -1910,11 +1713,9 @@ ccc
xa(i)=x1(i)
ya(i)=y1(i)
za(i)=z1(i)
- nresA(i)=n1(i)
xb(i)=x2(i)
yb(i)=y2(i)
zb(i)=z2(i)
- nresB(i)=n2(i)
iA(i)=i
iB(i)=i
enddo
@@ -1934,7 +1735,6 @@ ccc
*** iterative parameters ----->
n_it=20 !maximum number of iterations
- d_output=5 !for output alignment
n_init_max=6 !maximum number of L_init
n_init=0
L_ini_min=4
@@ -2079,10 +1879,8 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
*
* L1--Length of the first structure
* (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
-* n1(i)--Residue sequence number of i'th residue at the first structure
* L2--Length of the second structure
* (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
-* n2(i)--Residue sequence number of i'th residue at the second structure
* TM--TM-score of the comparison
* Rcomm--RMSD of two structures in the common aligned residues
* Lcomm--Length of the common aligned regions
@@ -2095,27 +1893,27 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
*************************************************************************
*************************************************************************
*** dis<8, but same search engine
- subroutine TMscore8(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,
+ subroutine TMscore8(dx,L1,x1,y1,z1,L2,x2,y2,z2,
& TM,Rcomm,Lcomm)
PARAMETER(nmax=5000)
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
+ common/nres/nseqA,nseqB
common/para/d,d0
common/d0min/d0_min
common/align/n_ali,iA(nmax),iB(nmax)
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
dimension k_ali(nmax),k_ali0(nmax)
- dimension L_ini(100),iq(nmax)
+ dimension L_ini(100)
common/scores/score
double precision score,score_max
dimension xa(nmax),ya(nmax),za(nmax)
- dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
- dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
+ dimension x1(nmax),y1(nmax),z1(nmax)
+ dimension x2(nmax),y2(nmax),z2(nmax)
ccc RMSD:
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms !armsd is real
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms !armsd is real
data w /nmax*1.0/
ccc
@@ -2127,11 +1925,9 @@ ccc
xa(i)=x1(i)
ya(i)=y1(i)
za(i)=z1(i)
- nresA(i)=n1(i)
xb(i)=x2(i)
yb(i)=y2(i)
zb(i)=z2(i)
- nresB(i)=n2(i)
iA(i)=i
iB(i)=i
enddo
@@ -2150,7 +1946,6 @@ ccc
if(d0_search.lt.4.5)d0_search=4.5
*** iterative parameters ----->
n_it=20 !maximum number of iterations
- d_output=5 !for output alignment
n_init_max=6 !maximum number of L_init
n_init=0
L_ini_min=4
@@ -2284,12 +2079,12 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
PARAMETER(nmax=5000)
common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
+ common/nres/nseqA,nseqB
common/para/d,d0
common/align/n_ali,iA(nmax),iB(nmax)
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
common/scores/score
- double precision score,score_max
+ double precision score
common/d8/d8
d_tmp=d
@@ -2323,10 +2118,8 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*
* L1--Length of the first structure
* (x1(i),y1(i),z1(i))--coordinates of i'th residue at the first structure
-* n1(i)--Residue sequence number of i'th residue at the first structure
* L2--Length of the second structure
* (x2(i),y2(i),z2(i))--coordinates of i'th residue at the second structure
-* n2(i)--Residue sequence number of i'th residue at the second structure
* TM--TM-score of the comparison
* Rcomm--RMSD of two structures in the common aligned residues
* Lcomm--Length of the common aligned regions
@@ -2339,27 +2132,27 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
*************************************************************************
*************************************************************************
*** normal TM-score:
- subroutine TMscore(dx,L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,
+ subroutine TMscore(dx,L1,x1,y1,z1,L2,x2,y2,z2,
& TM,Rcomm,Lcomm)
PARAMETER(nmax=5000)
common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
+ common/nres/nseqA,nseqB
common/para/d,d0
common/d0min/d0_min
common/align/n_ali,iA(nmax),iB(nmax)
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
dimension k_ali(nmax),k_ali0(nmax)
- dimension L_ini(100),iq(nmax)
+ dimension L_ini(100)
common/scores/score
double precision score,score_max
dimension xa(nmax),ya(nmax),za(nmax)
- dimension x1(nmax),y1(nmax),z1(nmax),n1(nmax)
- dimension x2(nmax),y2(nmax),z2(nmax),n2(nmax)
+ dimension x1(nmax),y1(nmax),z1(nmax)
+ dimension x2(nmax),y2(nmax),z2(nmax)
ccc RMSD:
- double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
- double precision u(3,3),t(3),rms,drms !armsd is real
+ double precision r_1(3,nmax),r_2(3,nmax),w(nmax)
+ double precision u(3,3),t(3),rms !armsd is real
data w /nmax*1.0/
ccc
@@ -2371,11 +2164,9 @@ ccc
xa(i)=x1(i)
ya(i)=y1(i)
za(i)=z1(i)
- nresA(i)=n1(i)
xb(i)=x2(i)
yb(i)=y2(i)
zb(i)=z2(i)
- nresB(i)=n2(i)
iA(i)=i
iB(i)=i
enddo
@@ -2395,7 +2186,6 @@ c d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
if(d0_search.lt.4.5)d0_search=4.5
*** iterative parameters ----->
n_it=20 !maximum number of iterations
- d_output=5 !for output alignment
n_init_max=6 !maximum number of L_init
n_init=0
L_ini_min=4
@@ -2529,7 +2319,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
PARAMETER(nmax=5000)
common/stru/xa(nmax),ya(nmax),za(nmax),xb(nmax),yb(nmax),zb(nmax)
- common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
+ common/nres/nseqA,nseqB
common/para/d,d0
common/align/n_ali,iA(nmax),iB(nmax)
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
@@ -2574,9 +2364,7 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
LOGICAL*1 DIR
common/dpc/score(nmax,nmax),gap_open,invmap(nmax)
dimension DIR(0:nmax,0:nmax),VAL(0:nmax,0:nmax)
- common/speed/align_len1
REAL H,V
- align_len1=0
*** initialize the matrix:
val(0,0)=0
do i=1,nseq1
@@ -2618,7 +2406,6 @@ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
DO WHILE((i.GT.0).AND.(j.GT.0))
IF(DIR(i,j))THEN
invmap(j)=i
- align_len1=align_len1+1
i=i-1
j=j-1
ELSE
@@ -2639,7 +2426,7 @@ c^^^^^^^^^^^^^^^Dynamical programming done ^^^^^^^^^^^^^^^^^^^
END
cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
-c w - w(m) is weight for atom pair c m (given)
+c w - w(m) is weight for atom pair c m (given)
c x - x(i,m) are coordinates of atom c m in set x (given)
c y - y(i,m) are coordinates of atom c m in set y (given)
c n - n is number of atom pairs (given)
@@ -2654,243 +2441,260 @@ c -2: no result obtained because of negative weights w
c or all weights equal to zero.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
- integer ip(9), ip2312(4), i, j, k, l, m1, m, ier, n, mode
- double precision w(1), x(3, 1), y(3, 1), u(3, 3), t(3), rms, sigma
- double precision r(3, 3), xc(3), yc(3), wc, a(3, 3), b(3, 3), e0,
- &e(3), e1, e2, e3, d, spur, det, cof, h, g, cth, sth, sqrth, p, tol
- &, rr(6), rr1, rr2, rr3, rr4, rr5, rr6, ss(6), ss1, ss2, ss3, ss4,
- &ss5, ss6, zero, one, two, three, sqrt3
- equivalence (rr(1), rr1), (rr(2), rr2), (rr(3), rr3), (rr(4), rr4)
- &, (rr(5), rr5), (rr(6), rr6), (ss(1), ss1), (ss(2), ss2), (ss(3),
- &ss3), (ss(4), ss4), (ss(5), ss5), (ss(6), ss6), (e(1), e1), (e(2)
- &, e2), (e(3), e3)
+ double precision w(*), x(3,*), y(3,*)
+ integer n, mode
+
+ double precision rms, u(3,3), t(3)
+ integer ier
+
+ integer i, j, k, l, m1, m
+ integer ip(9), ip2312(4)
+ double precision r(3,3), xc(3), yc(3), wc
+ double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
+ double precision e0, d, spur, det, cof, h, g
+ double precision cth, sth, sqrth, p, sigma
+
+ double precision sqrt3, tol, zero
+
data sqrt3 / 1.73205080756888d+00 /
data tol / 1.0d-2 /
data zero / 0.0d+00 /
- data one / 1.0d+00 /
- data two / 2.0d+00 /
- data three / 3.0d+00 /
data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
data ip2312 / 2, 3, 1, 2 /
-
-c 156 "rms.for"
wc = zero
- rms = 0.0
+ rms = zero
e0 = zero
- do 1 i = 1, 3
- xc(i) = zero
- yc(i) = zero
- t(i) = 0.0
- do 1 j = 1, 3
- d = zero
- if (i .eq. j) d = one
- u(i,j) = d
- a(i,j) = d
- 1 r(i,j) = zero
+
+ do i=1, 3
+ xc(i) = zero
+ yc(i) = zero
+ t(i) = zero
+ do j=1, 3
+ r(i,j) = zero
+ u(i,j) = zero
+ a(i,j) = zero
+ if( i .eq. j ) then
+ u(i,j) = 1.0
+ a(i,j) = 1.0
+ end if
+ end do
+ end do
+
ier = -1
-c**** DETERMINE CENTROIDS OF BOTH VECTOR SETS X AND Y
-c 170 "rms.for"
- if (n .lt. 1) return
-c 172 "rms.for"
+ if( n .lt. 1 ) return
ier = -2
- do 2 m = 1, n
- if (w(m) .lt. 0.0) return
- wc = wc + w(m)
- do 2 i = 1, 3
- xc(i) = xc(i) + (w(m) * x(i,m))
- 2 yc(i) = yc(i) + (w(m) * y(i,m))
- if (wc .le. zero) return
- do 3 i = 1, 3
- xc(i) = xc(i) / wc
-c**** DETERMINE CORRELATION MATRIX R BETWEEN VECTOR SETS Y AND X
-c 182 "rms.for"
- 3 yc(i) = yc(i) / wc
-c 184 "rms.for"
- do 4 m = 1, n
- do 4 i = 1, 3
- e0 = e0 + (w(m) * (((x(i,m) - xc(i)) ** 2) + ((y(i,m) - yc(i)) **
- &2)))
-c 187 "rms.for"
- d = w(m) * (y(i,m) - yc(i))
- do 4 j = 1, 3
-c**** CALCULATE DETERMINANT OF R(I,J)
-c 189 "rms.for"
- 4 r(i,j) = r(i,j) + (d * (x(j,m) - xc(j)))
-c 191 "rms.for"
- det = ((r(1,1) * ((r(2,2) * r(3,3)) - (r(2,3) * r(3,2)))) - (r(1,2
- &) * ((r(2,1) * r(3,3)) - (r(2,3) * r(3,1))))) + (r(1,3) * ((r(2,1)
- & * r(3,2)) - (r(2,2) * r(3,1))))
-c**** FORM UPPER TRIANGLE OF TRANSPOSED(R)*R
-c 194 "rms.for"
+ do m=1, n
+ if( w(m) .lt. 0.0 ) return
+ wc = wc + w(m)
+ do i=1, 3
+ xc(i) = xc(i) + w(m)*x(i,m)
+ yc(i) = yc(i) + w(m)*y(i,m)
+ end do
+ end do
+ if( wc .le. zero ) return
+ do i=1, 3
+ xc(i) = xc(i) / wc
+ yc(i) = yc(i) / wc
+ end do
+
+ do m=1, n
+ do i=1, 3
+ e0=e0+w(m)*((x(i,m)-xc(i))**2+(y(i,m)-yc(i))**2)
+ d = w(m) * ( y(i,m) - yc(i) )
+ do j=1, 3
+ r(i,j) = r(i,j) + d*( x(j,m) - xc(j) )
+ end do
+ end do
+ end do
+
+ det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
+ & - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
+ & + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
+
sigma = det
-c 196 "rms.for"
+
m = 0
- do 5 j = 1, 3
- do 5 i = 1, j
- m = m + 1
-c***************** EIGENVALUES *****************************************
-c**** FORM CHARACTERISTIC CUBIC X**3-3*SPUR*X**2+3*COF*X-DET=0
-c 200 "rms.for"
- 5 rr(m) = ((r(1,i) * r(1,j)) + (r(2,i) * r(2,j))) + (r(3,i) * r(3,j)
- &)
-c 203 "rms.for"
- spur = ((rr1 + rr3) + rr6) / three
- cof = ((((((rr3 * rr6) - (rr5 * rr5)) + (rr1 * rr6)) - (rr4 * rr4)
- &) + (rr1 * rr3)) - (rr2 * rr2)) / three
-c 205 "rms.for"
- det = det * det
- do 6 i = 1, 3
- 6 e(i) = spur
-c**** REDUCE CUBIC TO STANDARD FORM Y**3-3HY+2G=0 BY PUTTING X=Y+SPUR
-c 208 "rms.for"
- if (spur .le. zero) goto 40
-c 210 "rms.for"
- d = spur * spur
+ do j=1, 3
+ do i=1, j
+ m = m+1
+ rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
+ end do
+ end do
+
+ spur = (rr(1)+rr(3)+rr(6)) / 3.0
+ cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
+ & - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
+ det = det*det
+
+ do i=1, 3
+ e(i) = spur
+ end do
+ if( spur .le. zero ) goto 40
+ d = spur*spur
h = d - cof
-c**** SOLVE CUBIC. ROOTS ARE E1,E2,E3 IN DECREASING ORDER
-c 212 "rms.for"
- g = (((spur * cof) - det) / two) - (spur * h)
-c 214 "rms.for"
- if (h .le. zero) goto 8
+ g = (spur*cof - det)/2.0 - spur*h
+ if( h .le. zero ) then
+ if( mode .eq. 0 ) then
+ goto 50
+ else
+ goto 30
+ end if
+ end if
sqrth = dsqrt(h)
- d = ((h * h) * h) - (g * g)
- if (d .lt. zero) d = zero
- d = datan2(dsqrt(d),- g) / three
+ d = h*h*h - g*g
+ if( d .lt. zero ) d = zero
+ d = datan2( dsqrt(d), -g ) / 3.0
cth = sqrth * dcos(d)
- sth = (sqrth * sqrt3) * dsin(d)
- e1 = (spur + cth) + cth
- e2 = (spur - cth) + sth
- e3 = (spur - cth) - sth
-c.....HANDLE SPECIAL CASE OF 3 IDENTICAL ROOTS
-c 224 "rms.for"
- if (mode) 10, 50, 10
-c**************** EIGENVECTORS *****************************************
-c 226 "rms.for"
- 8 if (mode) 30, 50, 30
-c 228 "rms.for"
- 10 do 15 l = 1, 3, 2
- d = e(l)
- ss1 = ((d - rr3) * (d - rr6)) - (rr5 * rr5)
- ss2 = ((d - rr6) * rr2) + (rr4 * rr5)
- ss3 = ((d - rr1) * (d - rr6)) - (rr4 * rr4)
- ss4 = ((d - rr3) * rr4) + (rr2 * rr5)
- ss5 = ((d - rr1) * rr5) + (rr2 * rr4)
- ss6 = ((d - rr1) * (d - rr3)) - (rr2 * rr2)
- j = 1
- if (dabs(ss1) .ge. dabs(ss3)) goto 12
- j = 2
- if (dabs(ss3) .ge. dabs(ss6)) goto 13
- 11 j = 3
- goto 13
- 12 if (dabs(ss1) .lt. dabs(ss6)) goto 11
- 13 d = zero
- j = 3 * (j - 1)
- do 14 i = 1, 3
- k = ip(i + j)
- a(i,l) = ss(k)
- 14 d = d + (ss(k) * ss(k))
- if (d .gt. zero) d = one / dsqrt(d)
- do 15 i = 1, 3
- 15 a(i,l) = a(i,l) * d
- d = ((a(1,1) * a(1,3)) + (a(2,1) * a(2,3))) + (a(3,1) * a(3,3))
- m1 = 3
- m = 1
- if ((e1 - e2) .gt. (e2 - e3)) goto 16
- m1 = 1
- m = 3
- 16 p = zero
- do 17 i = 1, 3
- a(i,m1) = a(i,m1) - (d * a(i,m))
- 17 p = p + (a(i,m1) ** 2)
- if (p .le. tol) goto 19
- p = one / dsqrt(p)
- do 18 i = 1, 3
- 18 a(i,m1) = a(i,m1) * p
- goto 21
- 19 p = one
- do 20 i = 1, 3
- if (p .lt. dabs(a(i,m))) goto 20
- p = dabs(a(i,m))
- j = i
- 20 continue
- k = ip2312(j)
- l = ip2312(j + 1)
- p = dsqrt((a(k,m) ** 2) + (a(l,m) ** 2))
- if (p .le. tol) goto 40
- a(j,m1) = zero
- a(k,m1) = - (a(l,m) / p)
- a(l,m1) = a(k,m) / p
- 21 a(1,2) = (a(2,3) * a(3,1)) - (a(2,1) * a(3,3))
- a(2,2) = (a(3,3) * a(1,1)) - (a(3,1) * a(1,3))
-c****************** ROTATION MATRIX ************************************
-c 282 "rms.for"
- a(3,2) = (a(1,3) * a(2,1)) - (a(1,1) * a(2,3))
-c 284 "rms.for"
- 30 do 32 l = 1, 2
- d = zero
- do 31 i = 1, 3
- b(i,l) = ((r(i,1) * a(1,l)) + (r(i,2) * a(2,l))) + (r(i,3) * a(3,l
- &))
-c 288 "rms.for"
- 31 d = d + (b(i,l) ** 2)
- if (d .gt. zero) d = one / dsqrt(d)
- do 32 i = 1, 3
- 32 b(i,l) = b(i,l) * d
- d = ((b(1,1) * b(1,2)) + (b(2,1) * b(2,2))) + (b(3,1) * b(3,2))
+ sth = sqrth*sqrt3*dsin(d)
+ e(1) = (spur + cth) + cth
+ e(2) = (spur - cth) + sth
+ e(3) = (spur - cth) - sth
+
+ if( mode .eq. 0 ) then
+ goto 50
+ end if
+
+ do l=1, 3, 2
+ d = e(l)
+ ss(1) = (d-rr(3)) * (d-rr(6)) - rr(5)*rr(5)
+ ss(2) = (d-rr(6)) * rr(2) + rr(4)*rr(5)
+ ss(3) = (d-rr(1)) * (d-rr(6)) - rr(4)*rr(4)
+ ss(4) = (d-rr(3)) * rr(4) + rr(2)*rr(5)
+ ss(5) = (d-rr(1)) * rr(5) + rr(2)*rr(4)
+ ss(6) = (d-rr(1)) * (d-rr(3)) - rr(2)*rr(2)
+
+ if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
+ j=1
+ if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
+ else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
+ j = 2
+ else
+ j = 3
+ end if
+
+ d = zero
+ j = 3 * (j - 1)
+
+ do i=1, 3
+ k = ip(i+j)
+ a(i,l) = ss(k)
+ d = d + ss(k)*ss(k)
+ end do
+ if( d .gt. zero ) d = 1.0 / dsqrt(d)
+ do i=1, 3
+ a(i,l) = a(i,l) * d
+ end do
+ end do
+
+ d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
+ if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
+ m1 = 3
+ m = 1
+ else
+ m1 = 1
+ m = 3
+ endif
+
p = zero
- do 33 i = 1, 3
- b(i,2) = b(i,2) - (d * b(i,1))
- 33 p = p + (b(i,2) ** 2)
- if (p .le. tol) goto 35
- p = one / dsqrt(p)
- do 34 i = 1, 3
- 34 b(i,2) = b(i,2) * p
- goto 37
- 35 p = one
- do 36 i = 1, 3
- if (p .lt. dabs(b(i,1))) goto 36
- p = dabs(b(i,1))
- j = i
- 36 continue
- k = ip2312(j)
- l = ip2312(j + 1)
- p = dsqrt((b(k,1) ** 2) + (b(l,1) ** 2))
- if (p .le. tol) goto 40
- b(j,2) = zero
- b(k,2) = - (b(l,1) / p)
- b(l,2) = b(k,1) / p
- 37 b(1,3) = (b(2,1) * b(3,2)) - (b(2,2) * b(3,1))
- b(2,3) = (b(3,1) * b(1,2)) - (b(3,2) * b(1,1))
- b(3,3) = (b(1,1) * b(2,2)) - (b(1,2) * b(2,1))
- do 39 i = 1, 3
- do 39 j = 1, 3
-c****************** TRANSLATION VECTOR *********************************
-c 320 "rms.for"
- 39 u(i,j) = ((b(i,1) * a(j,1)) + (b(i,2) * a(j,2))) + (b(i,3) * a(j,3
- &))
- 40 do 41 i = 1, 3
-c********************** RMS ERROR **************************************
-c 323 "rms.for"
- 41 t(i) = ((yc(i) - (u(i,1) * xc(1))) - (u(i,2) * xc(2))) - (u(i,3)
- & * xc(3))
- 50 do 51 i = 1, 3
- if (e(i) .lt. zero) e(i) = zero
- 51 e(i) = dsqrt(e(i))
+ do i=1, 3
+ a(i,m1) = a(i,m1) - d*a(i,m)
+ p = p + a(i,m1)**2
+ end do
+ if( p .le. tol ) then
+ p = 1.0
+ do 21 i=1, 3
+ if (p .lt. dabs(a(i,m))) goto 21
+ p = dabs( a(i,m) )
+ j = i
+ 21 continue
+ k = ip2312(j)
+ l = ip2312(j+1)
+ p = dsqrt( a(k,m)**2 + a(l,m)**2 )
+ if( p .le. tol ) goto 40
+ a(j,m1) = zero
+ a(k,m1) = -a(l,m)/p
+ a(l,m1) = a(k,m)/p
+ else
+ p = 1.0 / dsqrt(p)
+ do i=1, 3
+ a(i,m1) = a(i,m1)*p
+ end do
+ end if
+
+ a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
+ a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
+ a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
+
+ 30 do l=1, 2
+ d = zero
+ do i=1, 3
+ b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
+ d = d + b(i,l)**2
+ end do
+ if( d .gt. zero ) d = 1.0 / dsqrt(d)
+ do i=1, 3
+ b(i,l) = b(i,l)*d
+ end do
+ end do
+ d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
+ p = zero
+
+ do i=1, 3
+ b(i,2) = b(i,2) - d*b(i,1)
+ p = p + b(i,2)**2
+ end do
+ if( p .le. tol ) then
+ p = 1.0
+ do 22 i=1, 3
+ if(p.lt.dabs(b(i,1)))goto 22
+ p = dabs( b(i,1) )
+ j = i
+ 22 continue
+ k = ip2312(j)
+ l = ip2312(j+1)
+ p = dsqrt( b(k,1)**2 + b(l,1)**2 )
+ if( p .le. tol ) goto 40
+ b(j,2) = zero
+ b(k,2) = -b(l,1)/p
+ b(l,2) = b(k,1)/p
+ else
+ p = 1.0 / dsqrt(p)
+ do i=1, 3
+ b(i,2) = b(i,2)*p
+ end do
+ end if
+
+ b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
+ b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
+ b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
+
+ do i=1, 3
+ do j=1, 3
+ u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
+ end do
+ end do
+
+ 40 do i=1, 3
+ t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
+ end do
+ 50 do i=1, 3
+ if( e(i) .lt. zero ) e(i) = zero
+ e(i) = dsqrt( e(i) )
+ end do
+
ier = 0
- if (e2 .le. (e1 * 1.0d-05)) ier = -1
- d = e3
- if (sigma .ge. 0.0) goto 52
- d = - d
- if ((e2 - e3) .le. (e1 * 1.0d-05)) ier = -1
- 52 d = (d + e2) + e1
+ if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
+
+ d = e(3)
+ if( sigma .lt. 0.0 ) then
+ d = - d
+ if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
+ end if
+ d = (d + e(2)) + e(1)
+
rms = (e0 - d) - d
- if (rms .lt. 0.0) rms = 0.0
-
- return
-c.....END U3B...........................................................
-c----------------------------------------------------------
-c THE END
-c----------------------------------------------------------
-c 338 "rms.for"
+ if( rms .lt. 0.0 ) rms = 0.0
+
+ return
end
diff --git a/modules/bindings/src/tmscore.f b/modules/bindings/src/tmscore.f
index 02886d7ab9ed606bd5e586682660347a3bc16900..86340daf7416f0f4bf14449a267767694ead8f55 100644
--- a/modules/bindings/src/tmscore.f
+++ b/modules/bindings/src/tmscore.f
@@ -1,19 +1,30 @@
*************************************************************************
* This program is to compare two protein structures and identify the
-* best superposition that has the maximum TM-score. The program can
-* be freely copied or modified.
-* For comments, please email to: yzhang@ku.edu
-*
+* best superposition that has the highest TM-score. Input structures
+* must be in the PDB format. By default, TM-score is normalized by
+* the second protein. Users can obtain a brief instruction by simply
+* running the program without arguments. For comments/suggestions,
+* please contact email: zhng@umich.edu.
+*
* Reference:
* Yang Zhang, Jeffrey Skolnick, Proteins, 2004 57:702-10.
-*
+*
+* Permission to use, copy, modify, and distribute this program for
+* any purpose, with or without fee, is hereby granted, provided that
+* the notices on the head, the reference information, and this
+* copyright notice appear in all copies or substantial portions of
+* the Software. It is provided "as is" without express or implied
+* warranty.
******************* Updating history ************************************
-* 2005/10/19: the program was reformed so that the score values
+* 2005/10/19: the program was reformed so that the score values.
* are not dependent on the specific compilers.
-* 2006/06/20: select 'A' if there is altLoc when reading PDB file
-* 2007/02/05: fix a bug with length<15 in TMscore_32
-* 2007/02/27: rotation matrix from Chain-1 to Chain-2 is added
-* 2007/12/06: GDT-HA score is added, fixed a bug for reading PDB
+* 2006/06/20: selected 'A' if there is altLoc when reading PDB file.
+* 2007/02/05: fixed a bug with length<15 in TMscore_32.
+* 2007/02/27: rotation matrix from Chain-1 to Chain-2 was added.
+* 2007/12/06: GDT-HA score was added, fixed a bug for reading PDB.
+* 2010/08/02: A new RMSD matrix was used and obsolete statement removed.
+* 2011/01/03: The length of pdb file names were extended to 500.
+* 2011/01/30: An open source license is attached to the program.
*************************************************************************
program TMscore
@@ -26,9 +37,9 @@
common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
dimension k_ali(nmax),k_ali0(nmax)
- character*100 fnam,pdb(100),outname
+ character*500 fnam,pdb(100),outname
character*3 aa(-1:20),seqA(nmax),seqB(nmax)
- character*100 s,du
+ character*500 s,du
character seq1A(nmax),seq1B(nmax),ali(nmax)
character sequenceA(nmax),sequenceB(nmax),sequenceM(nmax)
@@ -62,19 +73,22 @@ ccc
if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then
write(*,*)
write(*,*)'Brief instruction for running TM-score program:'
+ write(*,*)'(For detail: Zhang & Skolnick, Proteins, 2004',
+ & ' 57:702-10)'
write(*,*)
- write(*,*)'1. Run TM-score to compare ''model'' and ''native'':'
- write(*,*)' TMscore model native'
+ write(*,*)'1. Run TM-score to compare ''model'' and ''native',
+ & ':'
+ write(*,*)' >TMscore model native'
write(*,*)
- write(*,*)'2. Run TM-score with an assigned d0, e.g. 5 Angstroms:'
- write(*,*)' TMscore model native -d 5'
+ write(*,*)'2. Run TM-score with an assigned d0, e.g. 5',
+ & ' Angstroms:'
+ write(*,*)' >TMscore model native -d 5'
write(*,*)
write(*,*)'3. Run TM-score with superposition output, e.g. ',
& '''TM.sup'':'
- write(*,*)' TMscore model native -o TM.sup'
- write(*,*)
- write(*,*)'4, To view the superimposed structures by rasmol:'
- write(*,*)' rasmol -script TM.sup'
+ write(*,*)' >TMscore model native -o TM.sup'
+ write(*,*)' To view the superimposed structures by rasmol:'
+ write(*,*)' >rasmol -script TM.sup'
write(*,*)
goto 9999
endif
@@ -314,7 +328,6 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
- 303 format(i5,i5,i5,f17.14,f17.14,i5,f7.3)
if(it.eq.n_it)goto 302
if(n_cut.eq.ka)then
neq=0
@@ -337,17 +350,17 @@ c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
& '************************'
write(*,*)'* TM-SCORE ',
& ' *'
- write(*,*)'* A scoring function to assess the quality of protein',
- & ' structure predictions *'
+ write(*,*)'* A scoring function to assess the similarity of prot',
+ & 'ein structures *'
write(*,*)'* Based on statistics: ',
& ' *'
- write(*,*)'* 0.0 < TM-score < 0.17, Random predictions ',
- & ' *'
- write(*,*)'* 0.4 < TM-score < 1.00, Meaningful predictions',
+ write(*,*)'* 0.0 < TM-score < 0.17, random structural simi',
+ & 'larity *'
+ write(*,*)'* 0.5 < TM-score < 1.00, in about the same fold',
& ' *'
write(*,*)'* Reference: Yang Zhang and Jeffrey Skolnick, ',
& 'Proteins 2004 57: 702-710 *'
- write(*,*)'* For comments, please email to: yzhang@ku.edu ',
+ write(*,*)'* For comments, please email to: zhng@umich.edu ',
& ' *'
write(*,*)'*****************************************************',
& '************************'
@@ -648,240 +661,260 @@ c -2: no result obtained because of negative weights w
c or all weights equal to zero.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
- integer ip(9), ip2312(4), i, j, k, l, m1, m, ier, n, mode
- double precision w(1), x(3, 1), y(3, 1), u(3, 3), t(3), rms, sigma
- double precision r(3, 3), xc(3), yc(3), wc, a(3, 3), b(3, 3), e0,
- &e(3), e1, e2, e3, d, spur, det, cof, h, g, cth, sth, sqrth, p, tol
- &, rr(6), rr1, rr2, rr3, rr4, rr5, rr6, ss(6), ss1, ss2, ss3, ss4,
- &ss5, ss6, zero, one, two, three, sqrt3
- equivalence (rr(1), rr1), (rr(2), rr2), (rr(3), rr3), (rr(4), rr4)
- &, (rr(5), rr5), (rr(6), rr6), (ss(1), ss1), (ss(2), ss2), (ss(3),
- &ss3), (ss(4), ss4), (ss(5), ss5), (ss(6), ss6), (e(1), e1), (e(2)
- &, e2), (e(3), e3)
+ double precision w(*), x(3,*), y(3,*)
+ integer n, mode
+
+ double precision rms, u(3,3), t(3)
+ integer ier
+
+ integer i, j, k, l, m1, m
+ integer ip(9), ip2312(4)
+ double precision r(3,3), xc(3), yc(3), wc
+ double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
+ double precision e0, d, spur, det, cof, h, g
+ double precision cth, sth, sqrth, p, sigma
+
+ double precision sqrt3, tol, zero
+
data sqrt3 / 1.73205080756888d+00 /
data tol / 1.0d-2 /
data zero / 0.0d+00 /
- data one / 1.0d+00 /
- data two / 2.0d+00 /
- data three / 3.0d+00 /
data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
data ip2312 / 2, 3, 1, 2 /
-c 156 "rms.for"
+
wc = zero
- rms = 0.0
+ rms = zero
e0 = zero
- do 1 i = 1, 3
- xc(i) = zero
- yc(i) = zero
- t(i) = 0.0
- do 1 j = 1, 3
- d = zero
- if (i .eq. j) d = one
- u(i,j) = d
- a(i,j) = d
- 1 r(i,j) = zero
+
+ do i=1, 3
+ xc(i) = zero
+ yc(i) = zero
+ t(i) = zero
+ do j=1, 3
+ r(i,j) = zero
+ u(i,j) = zero
+ a(i,j) = zero
+ if( i .eq. j ) then
+ u(i,j) = 1.0
+ a(i,j) = 1.0
+ end if
+ end do
+ end do
+
ier = -1
-c**** DETERMINE CENTROIDS OF BOTH VECTOR SETS X AND Y
-c 170 "rms.for"
- if (n .lt. 1) return
-c 172 "rms.for"
+ if( n .lt. 1 ) return
ier = -2
- do 2 m = 1, n
- if (w(m) .lt. 0.0) return
- wc = wc + w(m)
- do 2 i = 1, 3
- xc(i) = xc(i) + (w(m) * x(i,m))
- 2 yc(i) = yc(i) + (w(m) * y(i,m))
- if (wc .le. zero) return
- do 3 i = 1, 3
- xc(i) = xc(i) / wc
-c**** DETERMINE CORRELATION MATRIX R BETWEEN VECTOR SETS Y AND X
-c 182 "rms.for"
- 3 yc(i) = yc(i) / wc
-c 184 "rms.for"
- do 4 m = 1, n
- do 4 i = 1, 3
- e0 = e0 + (w(m) * (((x(i,m) - xc(i)) ** 2) + ((y(i,m) - yc(i)) **
- &2)))
-c 187 "rms.for"
- d = w(m) * (y(i,m) - yc(i))
- do 4 j = 1, 3
-c**** CALCULATE DETERMINANT OF R(I,J)
-c 189 "rms.for"
- 4 r(i,j) = r(i,j) + (d * (x(j,m) - xc(j)))
-c 191 "rms.for"
- det = ((r(1,1) * ((r(2,2) * r(3,3)) - (r(2,3) * r(3,2)))) - (r(1,2
- &) * ((r(2,1) * r(3,3)) - (r(2,3) * r(3,1))))) + (r(1,3) * ((r(2,1)
- & * r(3,2)) - (r(2,2) * r(3,1))))
-c**** FORM UPPER TRIANGLE OF TRANSPOSED(R)*R
-c 194 "rms.for"
+ do m=1, n
+ if( w(m) .lt. 0.0 ) return
+ wc = wc + w(m)
+ do i=1, 3
+ xc(i) = xc(i) + w(m)*x(i,m)
+ yc(i) = yc(i) + w(m)*y(i,m)
+ end do
+ end do
+ if( wc .le. zero ) return
+ do i=1, 3
+ xc(i) = xc(i) / wc
+ yc(i) = yc(i) / wc
+ end do
+
+ do m=1, n
+ do i=1, 3
+ e0=e0+w(m)*((x(i,m)-xc(i))**2+(y(i,m)-yc(i))**2)
+ d = w(m) * ( y(i,m) - yc(i) )
+ do j=1, 3
+ r(i,j) = r(i,j) + d*( x(j,m) - xc(j) )
+ end do
+ end do
+ end do
+
+ det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
+ & - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
+ & + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
+
sigma = det
-c 196 "rms.for"
+
m = 0
- do 5 j = 1, 3
- do 5 i = 1, j
- m = m + 1
-c***************** EIGENVALUES *****************************************
-c**** FORM CHARACTERISTIC CUBIC X**3-3*SPUR*X**2+3*COF*X-DET=0
-c 200 "rms.for"
- 5 rr(m) = ((r(1,i) * r(1,j)) + (r(2,i) * r(2,j))) + (r(3,i) * r(3,j)
- &)
-c 203 "rms.for"
- spur = ((rr1 + rr3) + rr6) / three
- cof = ((((((rr3 * rr6) - (rr5 * rr5)) + (rr1 * rr6)) - (rr4 * rr4)
- &) + (rr1 * rr3)) - (rr2 * rr2)) / three
-c 205 "rms.for"
- det = det * det
- do 6 i = 1, 3
- 6 e(i) = spur
-c**** REDUCE CUBIC TO STANDARD FORM Y**3-3HY+2G=0 BY PUTTING X=Y+SPUR
-c 208 "rms.for"
- if (spur .le. zero) goto 40
-c 210 "rms.for"
- d = spur * spur
+ do j=1, 3
+ do i=1, j
+ m = m+1
+ rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
+ end do
+ end do
+
+ spur = (rr(1)+rr(3)+rr(6)) / 3.0
+ cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
+ & - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
+ det = det*det
+
+ do i=1, 3
+ e(i) = spur
+ end do
+ if( spur .le. zero ) goto 40
+ d = spur*spur
h = d - cof
-c**** SOLVE CUBIC. ROOTS ARE E1,E2,E3 IN DECREASING ORDER
-c 212 "rms.for"
- g = (((spur * cof) - det) / two) - (spur * h)
-c 214 "rms.for"
- if (h .le. zero) goto 8
+ g = (spur*cof - det)/2.0 - spur*h
+ if( h .le. zero ) then
+ if( mode .eq. 0 ) then
+ goto 50
+ else
+ goto 30
+ end if
+ end if
sqrth = dsqrt(h)
- d = ((h * h) * h) - (g * g)
- if (d .lt. zero) d = zero
- d = datan2(dsqrt(d),- g) / three
+ d = h*h*h - g*g
+ if( d .lt. zero ) d = zero
+ d = datan2( dsqrt(d), -g ) / 3.0
cth = sqrth * dcos(d)
- sth = (sqrth * sqrt3) * dsin(d)
- e1 = (spur + cth) + cth
- e2 = (spur - cth) + sth
- e3 = (spur - cth) - sth
-c.....HANDLE SPECIAL CASE OF 3 IDENTICAL ROOTS
-c 224 "rms.for"
- if (mode) 10, 50, 10
-c**************** EIGENVECTORS *****************************************
-c 226 "rms.for"
- 8 if (mode) 30, 50, 30
-c 228 "rms.for"
- 10 do 15 l = 1, 3, 2
- d = e(l)
- ss1 = ((d - rr3) * (d - rr6)) - (rr5 * rr5)
- ss2 = ((d - rr6) * rr2) + (rr4 * rr5)
- ss3 = ((d - rr1) * (d - rr6)) - (rr4 * rr4)
- ss4 = ((d - rr3) * rr4) + (rr2 * rr5)
- ss5 = ((d - rr1) * rr5) + (rr2 * rr4)
- ss6 = ((d - rr1) * (d - rr3)) - (rr2 * rr2)
- j = 1
- if (dabs(ss1) .ge. dabs(ss3)) goto 12
- j = 2
- if (dabs(ss3) .ge. dabs(ss6)) goto 13
- 11 j = 3
- goto 13
- 12 if (dabs(ss1) .lt. dabs(ss6)) goto 11
- 13 d = zero
- j = 3 * (j - 1)
- do 14 i = 1, 3
- k = ip(i + j)
- a(i,l) = ss(k)
- 14 d = d + (ss(k) * ss(k))
- if (d .gt. zero) d = one / dsqrt(d)
- do 15 i = 1, 3
- 15 a(i,l) = a(i,l) * d
- d = ((a(1,1) * a(1,3)) + (a(2,1) * a(2,3))) + (a(3,1) * a(3,3))
- m1 = 3
- m = 1
- if ((e1 - e2) .gt. (e2 - e3)) goto 16
- m1 = 1
- m = 3
- 16 p = zero
- do 17 i = 1, 3
- a(i,m1) = a(i,m1) - (d * a(i,m))
- 17 p = p + (a(i,m1) ** 2)
- if (p .le. tol) goto 19
- p = one / dsqrt(p)
- do 18 i = 1, 3
- 18 a(i,m1) = a(i,m1) * p
- goto 21
- 19 p = one
- do 20 i = 1, 3
- if (p .lt. dabs(a(i,m))) goto 20
- p = dabs(a(i,m))
- j = i
- 20 continue
- k = ip2312(j)
- l = ip2312(j + 1)
- p = dsqrt((a(k,m) ** 2) + (a(l,m) ** 2))
- if (p .le. tol) goto 40
- a(j,m1) = zero
- a(k,m1) = - (a(l,m) / p)
- a(l,m1) = a(k,m) / p
- 21 a(1,2) = (a(2,3) * a(3,1)) - (a(2,1) * a(3,3))
- a(2,2) = (a(3,3) * a(1,1)) - (a(3,1) * a(1,3))
-c****************** ROTATION MATRIX ************************************
-c 282 "rms.for"
- a(3,2) = (a(1,3) * a(2,1)) - (a(1,1) * a(2,3))
-c 284 "rms.for"
- 30 do 32 l = 1, 2
- d = zero
- do 31 i = 1, 3
- b(i,l) = ((r(i,1) * a(1,l)) + (r(i,2) * a(2,l))) + (r(i,3) * a(3,l
- &))
-c 288 "rms.for"
- 31 d = d + (b(i,l) ** 2)
- if (d .gt. zero) d = one / dsqrt(d)
- do 32 i = 1, 3
- 32 b(i,l) = b(i,l) * d
- d = ((b(1,1) * b(1,2)) + (b(2,1) * b(2,2))) + (b(3,1) * b(3,2))
+ sth = sqrth*sqrt3*dsin(d)
+ e(1) = (spur + cth) + cth
+ e(2) = (spur - cth) + sth
+ e(3) = (spur - cth) - sth
+
+ if( mode .eq. 0 ) then
+ goto 50
+ end if
+
+ do l=1, 3, 2
+ d = e(l)
+ ss(1) = (d-rr(3)) * (d-rr(6)) - rr(5)*rr(5)
+ ss(2) = (d-rr(6)) * rr(2) + rr(4)*rr(5)
+ ss(3) = (d-rr(1)) * (d-rr(6)) - rr(4)*rr(4)
+ ss(4) = (d-rr(3)) * rr(4) + rr(2)*rr(5)
+ ss(5) = (d-rr(1)) * rr(5) + rr(2)*rr(4)
+ ss(6) = (d-rr(1)) * (d-rr(3)) - rr(2)*rr(2)
+
+ if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
+ j=1
+ if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
+ else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
+ j = 2
+ else
+ j = 3
+ end if
+
+ d = zero
+ j = 3 * (j - 1)
+
+ do i=1, 3
+ k = ip(i+j)
+ a(i,l) = ss(k)
+ d = d + ss(k)*ss(k)
+ end do
+ if( d .gt. zero ) d = 1.0 / dsqrt(d)
+ do i=1, 3
+ a(i,l) = a(i,l) * d
+ end do
+ end do
+
+ d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
+ if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
+ m1 = 3
+ m = 1
+ else
+ m1 = 1
+ m = 3
+ endif
+
p = zero
- do 33 i = 1, 3
- b(i,2) = b(i,2) - (d * b(i,1))
- 33 p = p + (b(i,2) ** 2)
- if (p .le. tol) goto 35
- p = one / dsqrt(p)
- do 34 i = 1, 3
- 34 b(i,2) = b(i,2) * p
- goto 37
- 35 p = one
- do 36 i = 1, 3
- if (p .lt. dabs(b(i,1))) goto 36
- p = dabs(b(i,1))
- j = i
- 36 continue
- k = ip2312(j)
- l = ip2312(j + 1)
- p = dsqrt((b(k,1) ** 2) + (b(l,1) ** 2))
- if (p .le. tol) goto 40
- b(j,2) = zero
- b(k,2) = - (b(l,1) / p)
- b(l,2) = b(k,1) / p
- 37 b(1,3) = (b(2,1) * b(3,2)) - (b(2,2) * b(3,1))
- b(2,3) = (b(3,1) * b(1,2)) - (b(3,2) * b(1,1))
- b(3,3) = (b(1,1) * b(2,2)) - (b(1,2) * b(2,1))
- do 39 i = 1, 3
- do 39 j = 1, 3
-c****************** TRANSLATION VECTOR *********************************
-c 320 "rms.for"
- 39 u(i,j) = ((b(i,1) * a(j,1)) + (b(i,2) * a(j,2))) + (b(i,3) * a(j,3
- &))
- 40 do 41 i = 1, 3
-c********************** RMS ERROR **************************************
-c 323 "rms.for"
- 41 t(i) = ((yc(i) - (u(i,1) * xc(1))) - (u(i,2) * xc(2))) - (u(i,3)
- & * xc(3))
- 50 do 51 i = 1, 3
- if (e(i) .lt. zero) e(i) = zero
- 51 e(i) = dsqrt(e(i))
+ do i=1, 3
+ a(i,m1) = a(i,m1) - d*a(i,m)
+ p = p + a(i,m1)**2
+ end do
+ if( p .le. tol ) then
+ p = 1.0
+ do 21 i=1, 3
+ if (p .lt. dabs(a(i,m))) goto 21
+ p = dabs( a(i,m) )
+ j = i
+ 21 continue
+ k = ip2312(j)
+ l = ip2312(j+1)
+ p = dsqrt( a(k,m)**2 + a(l,m)**2 )
+ if( p .le. tol ) goto 40
+ a(j,m1) = zero
+ a(k,m1) = -a(l,m)/p
+ a(l,m1) = a(k,m)/p
+ else
+ p = 1.0 / dsqrt(p)
+ do i=1, 3
+ a(i,m1) = a(i,m1)*p
+ end do
+ end if
+
+ a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
+ a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
+ a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
+
+ 30 do l=1, 2
+ d = zero
+ do i=1, 3
+ b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
+ d = d + b(i,l)**2
+ end do
+ if( d .gt. zero ) d = 1.0 / dsqrt(d)
+ do i=1, 3
+ b(i,l) = b(i,l)*d
+ end do
+ end do
+ d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
+ p = zero
+
+ do i=1, 3
+ b(i,2) = b(i,2) - d*b(i,1)
+ p = p + b(i,2)**2
+ end do
+ if( p .le. tol ) then
+ p = 1.0
+ do 22 i=1, 3
+ if(p.lt.dabs(b(i,1)))goto 22
+ p = dabs( b(i,1) )
+ j = i
+ 22 continue
+ k = ip2312(j)
+ l = ip2312(j+1)
+ p = dsqrt( b(k,1)**2 + b(l,1)**2 )
+ if( p .le. tol ) goto 40
+ b(j,2) = zero
+ b(k,2) = -b(l,1)/p
+ b(l,2) = b(k,1)/p
+ else
+ p = 1.0 / dsqrt(p)
+ do i=1, 3
+ b(i,2) = b(i,2)*p
+ end do
+ end if
+
+ b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
+ b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
+ b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
+
+ do i=1, 3
+ do j=1, 3
+ u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
+ end do
+ end do
+
+ 40 do i=1, 3
+ t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
+ end do
+ 50 do i=1, 3
+ if( e(i) .lt. zero ) e(i) = zero
+ e(i) = dsqrt( e(i) )
+ end do
+
ier = 0
- if (e2 .le. (e1 * 1.0d-05)) ier = -1
- d = e3
- if (sigma .ge. 0.0) goto 52
- d = - d
- if ((e2 - e3) .le. (e1 * 1.0d-05)) ier = -1
- 52 d = (d + e2) + e1
+ if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
+
+ d = e(3)
+ if( sigma .lt. 0.0 ) then
+ d = - d
+ if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
+ end if
+ d = (d + e(2)) + e(1)
+
rms = (e0 - d) - d
- if (rms .lt. 0.0) rms = 0.0
- return
-c.....END U3B...........................................................
-c----------------------------------------------------------
-c THE END
-c----------------------------------------------------------
-c 338 "rms.for"
- end
\ No newline at end of file
+ if( rms .lt. 0.0 ) rms = 0.0
+
+ return
+ end