1 *************************************************************************
2 * This program is to compare two protein structures and identify the
3 * best superposition that has the highest TM-score. Input structures
4 * must be in the PDB format. By default, TM-score is normalized by
5 * the second protein. Users can obtain a brief instruction by simply
6 * running the program without arguments. For comments/suggestions,
7 * please contact email: zhng@umich.edu.
10 * Yang Zhang, Jeffrey Skolnick, Proteins, 2004 57:702-10.
12 * Permission to use, copy, modify, and distribute this program for
13 * any purpose, with or without fee, is hereby granted, provided that
14 * the notices on the head, the reference information, and this
15 * copyright notice appear in all copies or substantial portions of
16 * the Software. It is provided "as is" without express or implied
18 ******************* Updating history ************************************
19 * 2005/10/19: the program was reformed so that the score values.
20 * are not dependent on the specific compilers.
21 * 2006/06/20: selected 'A' if there is altLoc when reading PDB file.
22 * 2007/02/05: fixed a bug with length<15 in TMscore_32.
23 * 2007/02/27: rotation matrix from Chain-1 to Chain-2 was added.
24 * 2007/12/06: GDT-HA score was added, fixed a bug for reading PDB.
25 * 2010/08/02: A new RMSD matrix was used and obsolete statement removed.
26 * 2011/01/03: The length of pdb file names were extended to 500.
27 * 2011/01/30: An open source license is attached to the program.
28 * 2012/05/07: Improved RMSD calculation subroutine which speeds up
29 * TM-score program by 30%.
30 * 2012/06/05: Added option '-l L' which calculates TM-score (and maxsub
31 * and GDT scores) normalized by a specific length 'L'.
32 *************************************************************************
35 subroutine TMscore_sub(rmsd,gdt_ts,gdt_ha,tmscore,cfname,lprint)
38 include 'COMMON.CONTROL'
39 include 'COMMON.IOUNITS'
40 include 'COMMON.CHAIN'
41 include 'COMMON.INTERACT'
43 real*8 rmsd,gdt_ts,gdt_ha,tmscore
44 common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
45 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
46 common/para/d,d0,d0_fix
47 common/align/n_ali,iA(nmax),iB(nmax)
48 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
49 dimension k_ali(nmax),k_ali0(nmax)
51 character*500 fnam,pdb(100)!,outname
53 character*3 aa(-1:20),seqA(nmax),seqB(nmax)
55 character seq1A(nmax),seq1B(nmax),ali(nmax)
56 character sequenceA(nmax),sequenceB(nmax),sequenceM(nmax)
58 dimension L_ini(100),iq(nmax)
59 common/scores/score,score_maxsub,score_fix,score10
60 common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
61 double precision score,score_max,score_fix,score_fix_max
62 double precision score_maxsub,score10
63 dimension xa(nmax),ya(nmax),za(nmax)
66 double precision r_1(3,nmax),r_2(3,nmax),r_3(3,nmax),w(nmax)
67 double precision u(3,3),tt(3),rms,drms !armsd is real
73 data aa/ 'BCK','GLY','ALA','SER','CYS',
74 & 'VAL','THR','ILE','PRO','MET',
75 & 'ASP','ASN','LEU','LYS','GLU',
76 & 'GLN','ARG','HIS','PHE','TYR',
78 character*1 slc(-1:20)
79 data slc/'X','G','A','S','C',
80 & 'V','T','I','P','M',
81 & 'D','N','L','K','E',
82 & 'Q','R','H','F','Y',
85 *****instructions ----------------->
87 c if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then
89 c write(*,*)'Brief instruction for running TM-score program:'
90 c write(*,*)'(For detail: Zhang & Skolnick, Proteins, 2004',
93 c write(*,*)'1. Run TM-score to compare ''model'' and ',
95 c write(*,*)' >TMscore model native'
97 c write(*,*)'2. TM-score normalized with an assigned scale d0',
99 c write(*,*)' >TMscore model native -d 5'
101 c write(*,*)'3. TM-score normalized by a specific length, ',
103 c write(*,*)' >TMscore model native -l 120'
105 c write(*,*)'4. TM-score with superposition output, e.g. ',
107 c write(*,*)' >TMscore model native -o TM.sup'
108 c write(*,*)' To view the superimposed structures by rasmol:'
109 c write(*,*)' >rasmol -script TM.sup'
116 ******* options ----------->
125 c call getarg(i,fnam)
126 c if(fnam.eq.'-o')then
129 c call getarg(i,outname)
130 c elseif(fnam.eq.'-d')then
133 c call getarg(i,fnam)
135 c elseif(fnam.eq.'-l')then
138 c call getarg(i,fnam)
144 c if(i.lt.narg)goto 115
146 ccccccccc read data from first CA file:
147 c open(unit=10,file=pdb(1),status='old')
149 c 101 read(10,104,end=102) s
150 c if(s(1:3).eq.'TER') goto 102
151 c if(s(1:4).eq.'ATOM')then
152 c if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
154 c if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
156 c read(s,103)du,seqA(i),du,nresA(i),du,xa(i),ya(i),za(i)
158 c if(seqA(i).eq.aa(j))then
170 c 103 format(A17,A3,A2,i4,A4,3F8.3)
174 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
176 ccccccccc read data from first CA file:
177 c open(unit=10,file=pdb(2),status='old')
179 c 201 read(10,204,end=202) s
180 c if(s(1:3).eq.'TER') goto 202
181 c if(s(1:4).eq.'ATOM')then
182 c if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
184 c if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
186 c read(s,203)du,seqB(i),du,nresB(i),du,xb(i),yb(i),zb(i)
188 c if(seqB(i).eq.aa(j))then
200 c 203 format(A17,A3,A2,i4,A4,3F8.3)
204 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
206 ******************************************************************
207 * pickup the aligned residues:
208 ******************************************************************
212 c if(nresA(i).eq.nresB(j))then
221 c n_ali=k !number of aligned residues
223 c write(*,*)'There is no common residues in the input structures'
231 if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
234 if (itype(i).ne.ntyp1) then
245 d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
250 d0=1.24*(l0_fix-15)**(1.0/3.0)-1.8
253 if(m_fix.eq.1)d0=d0_fix
256 if(d0_search.gt.8)d0_search=8
257 if(d0_search.lt.4.5)d0_search=4.5
258 *** iterative parameters ----->
259 n_it=20 !maximum number of iterations
260 d_output=5 !for output alignment
261 if(m_fix.eq.1)d_output=d0_fix
262 n_init_max=6 !maximum number of L_init
265 if(n_ali.lt.4)L_ini_min=n_ali
268 L_ini(n_init)=n_ali/2**(n_init-1)
269 if(L_ini(n_init).le.L_ini_min)then
270 L_ini(n_init)=L_ini_min
275 L_ini(n_init)=L_ini_min
278 ******************************************************************
279 * find the maximum score starting from local structures superposition
280 ******************************************************************
281 score_max=-1 !TM-score
282 score_maxsub_max=-1 !MaxSub-score
283 score10_max=-1 !TM-score10
284 n_GDT05_max=-1 !number of residues<0.5
285 n_GDT1_max=-1 !number of residues<1
286 n_GDT2_max=-1 !number of residues<2
287 n_GDT4_max=-1 !number of residues<4
288 n_GDT8_max=-1 !number of residues<8
291 write (iout,*) "cref and ccref"
294 if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
297 if (itype(i).ne.ntyp1) then
302 xb(nnsup)=cref_pdb(1,i,1)
303 yb(nnsup)=cref_pdb(2,i,1)
304 zb(nnsup)=cref_pdb(3,i,1)
307 c ccref(j,nnsup)=cref_pdb(j,i,1)
310 write (iout,'(i5,3f10.5,5x,3f10.5)') nnsup,
311 & xa(nnsup),ya(nnsup),za(nnsup),xb(nnsup),yb(nnsup),zb(nnsup)
316 do 333 i_init=1,n_init
318 iL_max=n_ali-L_init+1
319 do 300 iL=1,iL_max !on aligned residues, [1,nseqA]
323 k=iL+i-1 ![1,n_ali] common aligned
334 if(i_init.eq.1)then !global superposition
335 call u3b(w,r_1,r_2,LL,2,rms,u,tt,ier) !0:rmsd; 1:u,t; 2:rmsd,u,t
339 call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
342 xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
343 yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
344 zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
347 call score_fun !init, get scores, n_cut+i_ali(i) for iteration
348 if(score_max.lt.score)then
355 if(score10_max.lt.score10)score10_max=score10
356 if(score_maxsub_max.lt.score_maxsub)score_maxsub_max=
358 if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
359 if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
360 if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
361 if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
362 if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
363 *** iteration for extending ---------------------------------->
369 m=i_ali(i) ![1,n_ali]
380 call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
382 xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
383 yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
384 zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
386 call score_fun !get scores, n_cut+i_ali(i) for iteration
387 if(score_max.lt.score)then
394 if(score10_max.lt.score10)score10_max=score10
395 if(score_maxsub_max.lt.score_maxsub)score_maxsub_max
397 if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
398 if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
399 if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
400 if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
401 if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
402 if(it.eq.n_it)goto 302
406 if(i_ali(i).eq.k_ali(i))neq=neq+1
408 if(n_cut.eq.neq)goto 302
410 301 continue !for iteration
412 300 continue !for shift
413 333 continue !for initial length, L_ali/M
417 ratio=float(nseqB)/float(l0_fix)
420 score_max=score_max*float(nseqB)/float(l0_fix)
422 score_GDT=(n_GDT1_max+n_GDT2_max+n_GDT4_max+n_GDT8_max)
424 score_GDT_HA=(n_GDT05_max+n_GDT1_max+n_GDT2_max+n_GDT4_max)
427 gdt_ts=score_GDT*ratio
428 gdt_ha=score_GDT_HA*ratio
430 ******************************************************************
432 ******************************************************************
433 *** output TM-scale ---------------------------->
438 write(iout,*)'**************************************************',
439 & '***************************'
440 write(iout,*)'* TM-SCORE ',
442 write(iout,*)'* A scoring function to assess the similarity of p',
443 & 'rotein structures *'
444 write(iout,*)'* Based on statistics: ',
446 write(iout,*)'* 0.0 < TM-score < 0.17, random structural s',
448 write(iout,*)'* 0.5 < TM-score < 1.00, in about the same f',
450 write(iout,*)'* Reference: Yang Zhang and Jeffrey Skolnick, ',
451 & 'Proteins 2004 57: 702-710 *'
452 write(iout,*)'* For comments, please email to: zhng@umich.edu ',
454 write(iout,*)'**************************************************',
455 & '***************************'
457 write(iout,501)pdb(1),nseqA
458 501 format('Structure1: ',A10,' Length= ',I4)
460 write(iout,411)pdb(2),nseqB
461 write(iout,412)l0_fix
463 write(iout,502)pdb(2),nseqB
465 411 format('Structure2: ',A10,' Length= ',I4)
466 412 format('TM-score is notmalized by ',I4)
467 502 format('Structure2: ',A10,' Length= ',I4,
468 & ' (by which all scores are normalized)')
470 503 format('Number of residues in common= ',I4)
471 write(iout,513)rmsd_ali
472 513 format('RMSD of the common residues= ',F8.3)
474 write(iout,504)score_max,d0
475 504 format('TM-score = ',f6.4,' (d0=',f5.2,')')
476 write(iout,505)score_maxsub_max*ratio
477 505 format('MaxSub-score= ',f6.4,' (d0= 3.50)')
478 write(iout,506)score_GDT*ratio,n_GDT1_max/float(nseqB)*ratio,
479 & n_GDT2_max/float(nseqB)*ratio,n_GDT4_max/float(nseqB)*ratio,
480 & n_GDT8_max/float(nseqB)*ratio
481 506 format('GDT-TS-score= ',f6.4,' %(d<1)=',f6.4,' %(d<2)=',f6.4,
482 $ ' %(d<4)=',f6.4,' %(d<8)=',f6.4)
483 write(iout,507)score_GDT_HA*ratio,n_GDT05_max/float(nseqB)*ratio,
484 & n_GDT1_max/float(nseqB)*ratio,n_GDT2_max/float(nseqB)*ratio,
485 & n_GDT4_max/float(nseqB)*ratio
486 507 format('GDT-HA-score= ',f6.4,' %(d<0.5)=',f6.4,' %(d<1)=',f6.4,
487 $ ' %(d<2)=',f6.4,' %(d<4)=',f6.4)
494 c------------------------------------------------------------------------
495 *** recall and output the superposition of maxiumum TM-score:
498 c m=k_ali0(i) !record of the best alignment
507 c call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
509 c xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
510 c yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
511 c zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
514 c********* extract rotation matrix ------------>
515 c write(*,*)'-------- rotation matrix to rotate Chain-1 to ',
517 c write(*,*)'i t(i) u(i,1) u(i,2) ',
520 c write(*,304)i,t(i),u(i,1),u(i,2),u(i,3)
523 cc xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
524 cc yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
525 cc zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
526 cc write(*,*)j,xt(j),yt(j),zt(j)
529 c 304 format(I2,f18.10,f15.10,f15.10,f15.10)
531 c********* rmsd in superposed regions --------------->
532 c d=d_output !for output
533 c call score_fun() !give i_ali(i), score_max=score now
536 c m=i_ali(i) ![1,nseqA]
545 c call u3b(w,r_1,r_2,LL,0,rms,u,t,ier)
546 c armsd=dsqrt(rms/LL)
549 c*** output rotated chain1 + chain2----->
550 c if(m_out.ne.1)goto 999
551 c OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
553 c 901 format('select ',I4)
554 c write(7,900)'load inline'
555 c write(7,900)'select atomno<1000'
556 cc write(7,900)'color [255,20,147]'
557 c write(7,900)'wireframe .45'
558 c write(7,900)'select none'
559 c write(7,900)'select atomno>1000'
560 cc write(7,900)'color [100,149,237]'
561 c write(7,900)'wireframe .15'
562 c write(7,900)'color white'
564 c write(7,901)nresA(iA(i_ali(i)))
565 c write(7,900)'color red'
567 c write(7,900)'select all'
569 c write(7,514)rmsd_ali
570 c 514 format('REMARK RMSD of the common residues=',F8.3)
571 c write(7,515)score_max,d0
572 c 515 format('REMARK TM-score=',f6.4,' (d0=',f5.2,')')
574 c write(7,1237)nresA(i),seqA(i),nresA(i),xt(i),yt(i),zt(i)
578 c write(7,1239)nresA(i-1),nresA(i)
581 c write(7,1237)2000+nresB(i),seqB(i),nresB(i),xb(i),yb(i),zb(i)
585 c write(7,1239)2000+nresB(i-1),2000+nresB(i)
587 c 1237 format('ATOM ',i5,' CA ',A3,I6,4X,3F8.3)
589 c 1239 format('CONECT',I5,I5)
592 c*** record aligned residues by i=[1,nseqA], for sequenceM()------------>
597 c j=iA(i_ali(i)) ![1,nseqA]
598 c k=iB(i_ali(i)) ![1,nseqB]
599 c dis=sqrt((xt(j)-xb(k))**2+(yt(j)-yb(k))**2+(zt(j)-zb(k))**2)
600 c if(dis.lt.d_output)then
604 c*******************************************************************
605 c*** output aligned sequences
610 c if(i.gt.nseqA.and.j.gt.nseqB)goto 802
611 c if(i.gt.nseqA.and.j.le.nseqB)then
614 c sequenceB(k)=seq1B(j)
619 c if(i.le.nseqA.and.j.gt.nseqB)then
621 c sequenceA(k)=seq1A(i)
627 c if(nresA(i).eq.nresB(j))then
629 c sequenceA(k)=seq1A(i)
630 c sequenceB(k)=seq1B(j)
639 c elseif(nresA(i).lt.nresB(j))then
641 c sequenceA(k)=seq1A(i)
646 c elseif(nresB(j).lt.nresA(i))then
649 c sequenceB(k)=seq1B(j)
656 c write(*,600)d_output,n_cut,rmsd
657 c 600 format('Superposition in the TM-score: Length(d<',f3.1,
658 c $ ')=',i3,' RMSD=',f6.2)
659 c write(*,603)d_output
660 c 603 format('(":" denotes the residue pairs of distance < ',f3.1,
662 c write(*,601)(sequenceA(i),i=1,k)
663 c write(*,601)(sequenceM(i),i=1,k)
664 c write(*,601)(sequenceB(i),i=1,k)
665 c write(*,602)(mod(i,10),i=1,k)
670 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
673 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
674 c 1, collect those residues with dis<d;
675 c 2, calculate score_GDT, score_maxsub, score_TM
676 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
680 common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
681 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
682 common/para/d,d0,d0_fix
683 common/align/n_ali,iA(nmax),iB(nmax)
684 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
685 common/scores/score,score_maxsub,score_fix,score10
686 common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
687 double precision score,score_max,score_fix,score_fix_max
688 double precision score_maxsub,score10
691 21 n_cut=0 !number of residue-pairs dis<d, for iteration
692 n_GDT05=0 !for GDT-score, # of dis<0.5
693 n_GDT1=0 !for GDT-score, # of dis<1
694 n_GDT2=0 !for GDT-score, # of dis<2
695 n_GDT4=0 !for GDT-score, # of dis<4
696 n_GDT8=0 !for GDT-score, # of dis<8
697 score_maxsub_sum=0 !Maxsub-score
699 score_sum10=0 !TMscore10
701 i=iA(k) ![1,nseqA] reoder number of structureA
703 dis=sqrt((xt(i)-xb(j))**2+(yt(i)-yb(j))**2+(zt(i)-zb(j))**2)
707 i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d
725 *** for MAXsub-score:
727 score_maxsub_sum=score_maxsub_sum+1/(1+(dis/3.5)**2)
730 score_sum=score_sum+1/(1+(dis/d0)**2)
733 score_sum10=score_sum10+1/(1+(dis/d0)**2)
736 if(n_cut.lt.3.and.n_ali.gt.3)then
740 score_maxsub=score_maxsub_sum/float(nseqB) !MAXsub-score
741 score=score_sum/float(nseqB) !TM-score
742 score10=score_sum10/float(nseqB) !TM-score10
747 cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
748 c w - w(m) is weight for atom pair c m (given)
749 c x - x(i,m) are coordinates of atom c m in set x (given)
750 c y - y(i,m) are coordinates of atom c m in set y (given)
751 c n - n is number of atom pairs (given)
752 c mode - 0:calculate rms only (given,short)
753 c 1:calculate u,t only (given,medium)
754 c 2:calculate rms,u,t (given,longer)
755 c rms - sum of w*(ux+t-y)**2 over all atom pairs (result)
756 c u - u(i,j) is rotation matrix for best superposition (result)
757 c t - t(i) is translation vector for best superposition (result)
758 c ier - 0: a unique optimal superposition has been determined(result)
759 c -1: superposition is not unique but optimal
760 c -2: no result obtained because of negative weights w
761 c or all weights equal to zero.
762 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
763 subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
764 double precision w(*), x(3,*), y(3,*)
767 double precision rms, u(3,3), t(3)
770 integer i, j, k, l, m1, m
771 integer ip(9), ip2312(4)
772 double precision r(3,3), xc(3), yc(3), wc
773 double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
774 double precision e0, d, spur, det, cof, h, g
775 double precision cth, sth, sqrth, p, sigma
776 double precision c1x, c1y, c1z, c2x, c2y, c2z
777 double precision s1x, s1y, s1z, s2x, s2y, s2z
778 double precision sxx, sxy, sxz, syx, syy, syz, szx, szy, szz
780 double precision sqrt3, tol, zero
782 data sqrt3 / 1.73205080756888d+00 /
784 data zero / 0.0d+00 /
785 data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
786 data ip2312 / 2, 3, 1, 2 /
823 if( n .lt. 1 ) return
863 if(mode.eq.2.or.mode.eq.0) then ! need rmsd
866 e0 = e0+ (x(i, m)-xc(i))**2 + (y(i, m)-yc(i))**2
871 r(1, 1) = sxx-s1x*s2x/n;
872 r(2, 1) = sxy-s1x*s2y/n;
873 r(3, 1) = sxz-s1x*s2z/n;
874 r(1, 2) = syx-s1y*s2x/n;
875 r(2, 2) = syy-s1y*s2y/n;
876 r(3, 2) = syz-s1y*s2z/n;
877 r(1, 3) = szx-s1z*s2x/n;
878 r(2, 3) = szy-s1z*s2y/n;
879 r(3, 3) = szz-s1z*s2z/n;
881 det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
882 & - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
883 & + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
891 rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
895 spur = (rr(1)+rr(3)+rr(6)) / 3.0
896 cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
897 & - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
903 if( spur .le. zero ) goto 40
906 g = (spur*cof - det)/2.0 - spur*h
907 if( h .le. zero ) then
908 if( mode .eq. 0 ) then
916 if( d .lt. zero ) d = zero
917 d = datan2( dsqrt(d), -g ) / 3.0
918 cth = sqrth * dcos(d)
919 sth = sqrth*sqrt3*dsin(d)
920 e(1) = (spur + cth) + cth
921 e(2) = (spur - cth) + sth
922 e(3) = (spur - cth) - sth
924 if( mode .eq. 0 ) then
930 ss(1) = (d-rr(3)) * (d-rr(6)) - rr(5)*rr(5)
931 ss(2) = (d-rr(6)) * rr(2) + rr(4)*rr(5)
932 ss(3) = (d-rr(1)) * (d-rr(6)) - rr(4)*rr(4)
933 ss(4) = (d-rr(3)) * rr(4) + rr(2)*rr(5)
934 ss(5) = (d-rr(1)) * rr(5) + rr(2)*rr(4)
935 ss(6) = (d-rr(1)) * (d-rr(3)) - rr(2)*rr(2)
937 if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
939 if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
940 else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
954 if( d .gt. zero ) d = 1.0 / dsqrt(d)
960 d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
961 if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
971 a(i,m1) = a(i,m1) - d*a(i,m)
974 if( p .le. tol ) then
977 if (p .lt. dabs(a(i,m))) goto 21
983 p = dsqrt( a(k,m)**2 + a(l,m)**2 )
984 if( p .le. tol ) goto 40
995 a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
996 a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
997 a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
1002 b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
1005 if( d .gt. zero ) d = 1.0 / dsqrt(d)
1010 d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
1014 b(i,2) = b(i,2) - d*b(i,1)
1017 if( p .le. tol ) then
1020 if(p.lt.dabs(b(i,1)))goto 22
1026 p = dsqrt( b(k,1)**2 + b(l,1)**2 )
1027 if( p .le. tol ) goto 40
1038 b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
1039 b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
1040 b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
1044 u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
1049 t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
1052 if( e(i) .lt. zero ) e(i) = zero
1053 e(i) = dsqrt( e(i) )
1057 if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
1060 if( sigma .lt. 0.0 ) then
1062 if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
1064 d = (d + e(2)) + e(1)
1066 if(mode .eq. 2.or.mode.eq.0) then ! need rmsd
1068 if( rms .lt. 0.0 ) rms = 0.0