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
69 integer ii,ipermmin,iperm
74 data aa/ 'BCK','GLY','ALA','SER','CYS',
75 & 'VAL','THR','ILE','PRO','MET',
76 & 'ASP','ASN','LEU','LYS','GLU',
77 & 'GLN','ARG','HIS','PHE','TYR',
79 character*1 slc(-1:20)
80 data slc/'X','G','A','S','C',
81 & 'V','T','I','P','M',
82 & 'D','N','L','K','E',
83 & 'Q','R','H','F','Y',
86 *****instructions ----------------->
88 c if(fnam.eq.' '.or.fnam.eq.'?'.or.fnam.eq.'-h')then
90 c write(*,*)'Brief instruction for running TM-score program:'
91 c write(*,*)'(For detail: Zhang & Skolnick, Proteins, 2004',
94 c write(*,*)'1. Run TM-score to compare ''model'' and ',
96 c write(*,*)' >TMscore model native'
98 c write(*,*)'2. TM-score normalized with an assigned scale d0',
100 c write(*,*)' >TMscore model native -d 5'
102 c write(*,*)'3. TM-score normalized by a specific length, ',
104 c write(*,*)' >TMscore model native -l 120'
106 c write(*,*)'4. TM-score with superposition output, e.g. ',
108 c write(*,*)' >TMscore model native -o TM.sup'
109 c write(*,*)' To view the superimposed structures by rasmol:'
110 c write(*,*)' >rasmol -script TM.sup'
117 ******* options ----------->
126 c call getarg(i,fnam)
127 c if(fnam.eq.'-o')then
130 c call getarg(i,outname)
131 c elseif(fnam.eq.'-d')then
134 c call getarg(i,fnam)
136 c elseif(fnam.eq.'-l')then
139 c call getarg(i,fnam)
145 c if(i.lt.narg)goto 115
147 ccccccccc read data from first CA file:
148 c open(unit=10,file=pdb(1),status='old')
150 c 101 read(10,104,end=102) s
151 c if(s(1:3).eq.'TER') goto 102
152 c if(s(1:4).eq.'ATOM')then
153 c if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
155 c if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
157 c read(s,103)du,seqA(i),du,nresA(i),du,xa(i),ya(i),za(i)
159 c if(seqA(i).eq.aa(j))then
171 c 103 format(A17,A3,A2,i4,A4,3F8.3)
175 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
177 ccccccccc read data from first CA file:
178 c open(unit=10,file=pdb(2),status='old')
180 c 201 read(10,204,end=202) s
181 c if(s(1:3).eq.'TER') goto 202
182 c if(s(1:4).eq.'ATOM')then
183 c if(s(13:16).eq.'CA '.or.s(13:16).eq.' CA '.or.s(13:16).
185 c if(s(17:17).eq.' '.or.s(17:17).eq.'A')then
187 c read(s,203)du,seqB(i),du,nresB(i),du,xb(i),yb(i),zb(i)
189 c if(seqB(i).eq.aa(j))then
201 c 203 format(A17,A3,A2,i4,A4,3F8.3)
205 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
207 ******************************************************************
208 * pickup the aligned residues:
209 ******************************************************************
213 c if(nresA(i).eq.nresB(j))then
222 c n_ali=k !number of aligned residues
224 c write(*,*)'There is no common residues in the input structures'
235 if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
238 if (itype(i).ne.ntyp1) then
249 d0=1.24*(nseqB-15)**(1.0/3.0)-1.8
254 d0=1.24*(l0_fix-15)**(1.0/3.0)-1.8
257 if(m_fix.eq.1)d0=d0_fix
260 if(d0_search.gt.8)d0_search=8
261 if(d0_search.lt.4.5)d0_search=4.5
262 *** iterative parameters ----->
263 n_it=20 !maximum number of iterations
264 d_output=5 !for output alignment
265 if(m_fix.eq.1)d_output=d0_fix
266 n_init_max=6 !maximum number of L_init
269 if(n_ali.lt.4)L_ini_min=n_ali
272 L_ini(n_init)=n_ali/2**(n_init-1)
273 if(L_ini(n_init).le.L_ini_min)then
274 L_ini(n_init)=L_ini_min
279 L_ini(n_init)=L_ini_min
282 ******************************************************************
283 * find the maximum score starting from local structures superposition
284 ******************************************************************
285 score_max=-1 !TM-score
286 score_maxsub_max=-1 !MaxSub-score
287 score10_max=-1 !TM-score10
288 n_GDT05_max=-1 !number of residues<0.5
289 n_GDT1_max=-1 !number of residues<1
290 n_GDT2_max=-1 !number of residues<2
291 n_GDT4_max=-1 !number of residues<4
292 n_GDT8_max=-1 !number of residues<8
295 write (iout,*) "cref and ccref"
298 if (nres.gt.nsup+nnt-1) noverlap=nsup+nnt-1
301 if (itype(i).ne.ntyp1) then
303 xa(nnsup)=c(1,iperm(i,ii))
304 ya(nnsup)=c(2,iperm(i,ii))
305 za(nnsup)=c(3,iperm(i,ii))
306 xb(nnsup)=cref_pdb(1,i)
307 yb(nnsup)=cref_pdb(2,i)
308 zb(nnsup)=cref_pdb(3,i)
311 c ccref(j,nnsup)=cref_pdb(j,i,1)
314 write (iout,'(i5,3f10.5,5x,3f10.5)') nnsup,
315 & xa(nnsup),ya(nnsup),za(nnsup),xb(nnsup),yb(nnsup),zb(nnsup)
320 do 333 i_init=1,n_init
322 iL_max=n_ali-L_init+1
323 do 300 iL=1,iL_max !on aligned residues, [1,nseqA]
327 k=iL+i-1 ![1,n_ali] common aligned
338 if(i_init.eq.1)then !global superposition
339 call u3b(w,r_1,r_2,LL,2,rms,u,tt,ier) !0:rmsd; 1:u,t; 2:rmsd,u,t
343 call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
346 xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
347 yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
348 zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
351 call score_fun !init, get scores, n_cut+i_ali(i) for iteration
352 if(score_max.lt.score)then
359 if(score10_max.lt.score10)score10_max=score10
360 if(score_maxsub_max.lt.score_maxsub)score_maxsub_max=
362 if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
363 if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
364 if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
365 if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
366 if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
367 *** iteration for extending ---------------------------------->
373 m=i_ali(i) ![1,n_ali]
384 call u3b(w,r_1,r_2,LL,1,rms,u,tt,ier) !u rotate r_1 to r_2
386 xt(j)=tt(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
387 yt(j)=tt(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
388 zt(j)=tt(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
390 call score_fun !get scores, n_cut+i_ali(i) for iteration
391 if(score_max.lt.score)then
398 if(score10_max.lt.score10)score10_max=score10
399 if(score_maxsub_max.lt.score_maxsub)score_maxsub_max
401 if(n_GDT05_max.lt.n_GDT05)n_GDT05_max=n_GDT05
402 if(n_GDT1_max.lt.n_GDT1)n_GDT1_max=n_GDT1
403 if(n_GDT2_max.lt.n_GDT2)n_GDT2_max=n_GDT2
404 if(n_GDT4_max.lt.n_GDT4)n_GDT4_max=n_GDT4
405 if(n_GDT8_max.lt.n_GDT8)n_GDT8_max=n_GDT8
406 if(it.eq.n_it)goto 302
410 if(i_ali(i).eq.k_ali(i))neq=neq+1
412 if(n_cut.eq.neq)goto 302
414 301 continue !for iteration
416 300 continue !for shift
417 333 continue !for initial length, L_ali/M
421 ratio=float(nseqB)/float(l0_fix)
424 score_max=score_max*float(nseqB)/float(l0_fix)
426 score_GDT=(n_GDT1_max+n_GDT2_max+n_GDT4_max+n_GDT8_max)
428 score_GDT_HA=(n_GDT05_max+n_GDT1_max+n_GDT2_max+n_GDT4_max)
431 gdt_ts=score_GDT*ratio
432 gdt_ha=score_GDT_HA*ratio
435 if (ii.eq.1 .or. rmsd.lt.rmsd_min) then
450 ******************************************************************
452 ******************************************************************
453 *** output TM-scale ---------------------------->
458 write(iout,*)'**************************************************',
459 & '***************************'
460 write(iout,*)'* TM-SCORE ',
462 write(iout,*)'* A scoring function to assess the similarity of p',
463 & 'rotein structures *'
464 write(iout,*)'* Based on statistics: ',
466 write(iout,*)'* 0.0 < TM-score < 0.17, random structural s',
468 write(iout,*)'* 0.5 < TM-score < 1.00, in about the same f',
470 write(iout,*)'* Reference: Yang Zhang and Jeffrey Skolnick, ',
471 & 'Proteins 2004 57: 702-710 *'
472 write(iout,*)'* For comments, please email to: zhng@umich.edu ',
474 write(iout,*)'**************************************************',
475 & '***************************'
477 write(iout,501)pdb(1),nseqA
478 501 format('Structure1: ',A10,' Length= ',I4)
480 write(iout,411)pdb(2),nseqB
481 write(iout,412)l0_fix
483 write(iout,502)pdb(2),nseqB
485 411 format('Structure2: ',A10,' Length= ',I4)
486 412 format('TM-score is notmalized by ',I4)
487 502 format('Structure2: ',A10,' Length= ',I4,
488 & ' (by which all scores are normalized)')
490 503 format('Number of residues in common= ',I4)
491 write(iout,513)rmsd_ali
492 513 format('RMSD of the common residues= ',F8.3)
494 write(iout,504)score_max,d0
495 504 format('TM-score = ',f6.4,' (d0=',f5.2,')')
496 write(iout,505)score_maxsub_max*ratio
497 505 format('MaxSub-score= ',f6.4,' (d0= 3.50)')
498 write(iout,506)score_GDT*ratio,n_GDT1_max/float(nseqB)*ratio,
499 & n_GDT2_max/float(nseqB)*ratio,n_GDT4_max/float(nseqB)*ratio,
500 & n_GDT8_max/float(nseqB)*ratio
501 506 format('GDT-TS-score= ',f6.4,' %(d<1)=',f6.4,' %(d<2)=',f6.4,
502 $ ' %(d<4)=',f6.4,' %(d<8)=',f6.4)
503 write(iout,507)score_GDT_HA*ratio,n_GDT05_max/float(nseqB)*ratio,
504 & n_GDT1_max/float(nseqB)*ratio,n_GDT2_max/float(nseqB)*ratio,
505 & n_GDT4_max/float(nseqB)*ratio
506 507 format('GDT-HA-score= ',f6.4,' %(d<0.5)=',f6.4,' %(d<1)=',f6.4,
507 $ ' %(d<2)=',f6.4,' %(d<4)=',f6.4)
508 write (iout,*) "Permutation",ipermmin
515 c------------------------------------------------------------------------
516 *** recall and output the superposition of maxiumum TM-score:
519 c m=k_ali0(i) !record of the best alignment
528 c call u3b(w,r_1,r_2,LL,1,rms,u,t,ier) !u rotate r_1 to r_2
530 c xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
531 c yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
532 c zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
535 c********* extract rotation matrix ------------>
536 c write(*,*)'-------- rotation matrix to rotate Chain-1 to ',
538 c write(*,*)'i t(i) u(i,1) u(i,2) ',
541 c write(*,304)i,t(i),u(i,1),u(i,2),u(i,3)
544 cc xt(j)=t(1)+u(1,1)*xa(j)+u(1,2)*ya(j)+u(1,3)*za(j)
545 cc yt(j)=t(2)+u(2,1)*xa(j)+u(2,2)*ya(j)+u(2,3)*za(j)
546 cc zt(j)=t(3)+u(3,1)*xa(j)+u(3,2)*ya(j)+u(3,3)*za(j)
547 cc write(*,*)j,xt(j),yt(j),zt(j)
550 c 304 format(I2,f18.10,f15.10,f15.10,f15.10)
552 c********* rmsd in superposed regions --------------->
553 c d=d_output !for output
554 c call score_fun() !give i_ali(i), score_max=score now
557 c m=i_ali(i) ![1,nseqA]
566 c call u3b(w,r_1,r_2,LL,0,rms,u,t,ier)
567 c armsd=dsqrt(rms/LL)
570 c*** output rotated chain1 + chain2----->
571 c if(m_out.ne.1)goto 999
572 c OPEN(unit=7,file=outname,status='unknown') !pdb1.aln + pdb2.aln
574 c 901 format('select ',I4)
575 c write(7,900)'load inline'
576 c write(7,900)'select atomno<1000'
577 cc write(7,900)'color [255,20,147]'
578 c write(7,900)'wireframe .45'
579 c write(7,900)'select none'
580 c write(7,900)'select atomno>1000'
581 cc write(7,900)'color [100,149,237]'
582 c write(7,900)'wireframe .15'
583 c write(7,900)'color white'
585 c write(7,901)nresA(iA(i_ali(i)))
586 c write(7,900)'color red'
588 c write(7,900)'select all'
590 c write(7,514)rmsd_ali
591 c 514 format('REMARK RMSD of the common residues=',F8.3)
592 c write(7,515)score_max,d0
593 c 515 format('REMARK TM-score=',f6.4,' (d0=',f5.2,')')
595 c write(7,1237)nresA(i),seqA(i),nresA(i),xt(i),yt(i),zt(i)
599 c write(7,1239)nresA(i-1),nresA(i)
602 c write(7,1237)2000+nresB(i),seqB(i),nresB(i),xb(i),yb(i),zb(i)
606 c write(7,1239)2000+nresB(i-1),2000+nresB(i)
608 c 1237 format('ATOM ',i5,' CA ',A3,I6,4X,3F8.3)
610 c 1239 format('CONECT',I5,I5)
613 c*** record aligned residues by i=[1,nseqA], for sequenceM()------------>
618 c j=iA(i_ali(i)) ![1,nseqA]
619 c k=iB(i_ali(i)) ![1,nseqB]
620 c dis=sqrt((xt(j)-xb(k))**2+(yt(j)-yb(k))**2+(zt(j)-zb(k))**2)
621 c if(dis.lt.d_output)then
625 c*******************************************************************
626 c*** output aligned sequences
631 c if(i.gt.nseqA.and.j.gt.nseqB)goto 802
632 c if(i.gt.nseqA.and.j.le.nseqB)then
635 c sequenceB(k)=seq1B(j)
640 c if(i.le.nseqA.and.j.gt.nseqB)then
642 c sequenceA(k)=seq1A(i)
648 c if(nresA(i).eq.nresB(j))then
650 c sequenceA(k)=seq1A(i)
651 c sequenceB(k)=seq1B(j)
660 c elseif(nresA(i).lt.nresB(j))then
662 c sequenceA(k)=seq1A(i)
667 c elseif(nresB(j).lt.nresA(i))then
670 c sequenceB(k)=seq1B(j)
677 c write(*,600)d_output,n_cut,rmsd
678 c 600 format('Superposition in the TM-score: Length(d<',f3.1,
679 c $ ')=',i3,' RMSD=',f6.2)
680 c write(*,603)d_output
681 c 603 format('(":" denotes the residue pairs of distance < ',f3.1,
683 c write(*,601)(sequenceA(i),i=1,k)
684 c write(*,601)(sequenceM(i),i=1,k)
685 c write(*,601)(sequenceB(i),i=1,k)
686 c write(*,602)(mod(i,10),i=1,k)
691 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
694 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
695 c 1, collect those residues with dis<d;
696 c 2, calculate score_GDT, score_maxsub, score_TM
697 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
701 common/stru/xt(nmax),yt(nmax),zt(nmax),xb(nmax),yb(nmax),zb(nmax)
702 common/nres/nresA(nmax),nresB(nmax),nseqA,nseqB
703 common/para/d,d0,d0_fix
704 common/align/n_ali,iA(nmax),iB(nmax)
705 common/nscore/i_ali(nmax),n_cut ![1,n_ali],align residues for the score
706 common/scores/score,score_maxsub,score_fix,score10
707 common/GDT/n_GDT05,n_GDT1,n_GDT2,n_GDT4,n_GDT8
708 double precision score,score_max,score_fix,score_fix_max
709 double precision score_maxsub,score10
712 21 n_cut=0 !number of residue-pairs dis<d, for iteration
713 n_GDT05=0 !for GDT-score, # of dis<0.5
714 n_GDT1=0 !for GDT-score, # of dis<1
715 n_GDT2=0 !for GDT-score, # of dis<2
716 n_GDT4=0 !for GDT-score, # of dis<4
717 n_GDT8=0 !for GDT-score, # of dis<8
718 score_maxsub_sum=0 !Maxsub-score
720 score_sum10=0 !TMscore10
722 i=iA(k) ![1,nseqA] reoder number of structureA
724 dis=sqrt((xt(i)-xb(j))**2+(yt(i)-yb(j))**2+(zt(i)-zb(j))**2)
728 i_ali(n_cut)=k ![1,n_ali], mark the residue-pairs in dis<d
746 *** for MAXsub-score:
748 score_maxsub_sum=score_maxsub_sum+1/(1+(dis/3.5)**2)
751 score_sum=score_sum+1/(1+(dis/d0)**2)
754 score_sum10=score_sum10+1/(1+(dis/d0)**2)
757 if(n_cut.lt.3.and.n_ali.gt.3)then
761 score_maxsub=score_maxsub_sum/float(nseqB) !MAXsub-score
762 score=score_sum/float(nseqB) !TM-score
763 score10=score_sum10/float(nseqB) !TM-score10
768 cccccccccccccccc Calculate sum of (r_d-r_m)^2 cccccccccccccccccccccccccc
769 c w - w(m) is weight for atom pair c m (given)
770 c x - x(i,m) are coordinates of atom c m in set x (given)
771 c y - y(i,m) are coordinates of atom c m in set y (given)
772 c n - n is number of atom pairs (given)
773 c mode - 0:calculate rms only (given,short)
774 c 1:calculate u,t only (given,medium)
775 c 2:calculate rms,u,t (given,longer)
776 c rms - sum of w*(ux+t-y)**2 over all atom pairs (result)
777 c u - u(i,j) is rotation matrix for best superposition (result)
778 c t - t(i) is translation vector for best superposition (result)
779 c ier - 0: a unique optimal superposition has been determined(result)
780 c -1: superposition is not unique but optimal
781 c -2: no result obtained because of negative weights w
782 c or all weights equal to zero.
783 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
784 subroutine u3b(w, x, y, n, mode, rms, u, t, ier)
785 double precision w(*), x(3,*), y(3,*)
788 double precision rms, u(3,3), t(3)
791 integer i, j, k, l, m1, m
792 integer ip(9), ip2312(4)
793 double precision r(3,3), xc(3), yc(3), wc
794 double precision a(3,3), b(3,3), e(3), rr(6), ss(6)
795 double precision e0, d, spur, det, cof, h, g
796 double precision cth, sth, sqrth, p, sigma
797 double precision c1x, c1y, c1z, c2x, c2y, c2z
798 double precision s1x, s1y, s1z, s2x, s2y, s2z
799 double precision sxx, sxy, sxz, syx, syy, syz, szx, szy, szz
801 double precision sqrt3, tol, zero
803 data sqrt3 / 1.73205080756888d+00 /
805 data zero / 0.0d+00 /
806 data ip / 1, 2, 4, 2, 3, 5, 4, 5, 6 /
807 data ip2312 / 2, 3, 1, 2 /
844 if( n .lt. 1 ) return
884 if(mode.eq.2.or.mode.eq.0) then ! need rmsd
887 e0 = e0+ (x(i, m)-xc(i))**2 + (y(i, m)-yc(i))**2
892 r(1, 1) = sxx-s1x*s2x/n;
893 r(2, 1) = sxy-s1x*s2y/n;
894 r(3, 1) = sxz-s1x*s2z/n;
895 r(1, 2) = syx-s1y*s2x/n;
896 r(2, 2) = syy-s1y*s2y/n;
897 r(3, 2) = syz-s1y*s2z/n;
898 r(1, 3) = szx-s1z*s2x/n;
899 r(2, 3) = szy-s1z*s2y/n;
900 r(3, 3) = szz-s1z*s2z/n;
902 det = r(1,1) * ( (r(2,2)*r(3,3)) - (r(2,3)*r(3,2)) )
903 & - r(1,2) * ( (r(2,1)*r(3,3)) - (r(2,3)*r(3,1)) )
904 & + r(1,3) * ( (r(2,1)*r(3,2)) - (r(2,2)*r(3,1)) )
912 rr(m) = r(1,i)*r(1,j) + r(2,i)*r(2,j) + r(3,i)*r(3,j)
916 spur = (rr(1)+rr(3)+rr(6)) / 3.0
917 cof = (((((rr(3)*rr(6) - rr(5)*rr(5)) + rr(1)*rr(6))
918 & - rr(4)*rr(4)) + rr(1)*rr(3)) - rr(2)*rr(2)) / 3.0
924 if( spur .le. zero ) goto 40
927 g = (spur*cof - det)/2.0 - spur*h
928 if( h .le. zero ) then
929 if( mode .eq. 0 ) then
937 if( d .lt. zero ) d = zero
938 d = datan2( dsqrt(d), -g ) / 3.0
939 cth = sqrth * dcos(d)
940 sth = sqrth*sqrt3*dsin(d)
941 e(1) = (spur + cth) + cth
942 e(2) = (spur - cth) + sth
943 e(3) = (spur - cth) - sth
945 if( mode .eq. 0 ) then
951 ss(1) = (d-rr(3)) * (d-rr(6)) - rr(5)*rr(5)
952 ss(2) = (d-rr(6)) * rr(2) + rr(4)*rr(5)
953 ss(3) = (d-rr(1)) * (d-rr(6)) - rr(4)*rr(4)
954 ss(4) = (d-rr(3)) * rr(4) + rr(2)*rr(5)
955 ss(5) = (d-rr(1)) * rr(5) + rr(2)*rr(4)
956 ss(6) = (d-rr(1)) * (d-rr(3)) - rr(2)*rr(2)
958 if( dabs(ss(1)) .ge. dabs(ss(3)) ) then
960 if( dabs(ss(1)) .lt. dabs(ss(6)) ) j = 3
961 else if( dabs(ss(3)) .ge. dabs(ss(6)) ) then
975 if( d .gt. zero ) d = 1.0 / dsqrt(d)
981 d = a(1,1)*a(1,3) + a(2,1)*a(2,3) + a(3,1)*a(3,3)
982 if ((e(1) - e(2)) .gt. (e(2) - e(3))) then
992 a(i,m1) = a(i,m1) - d*a(i,m)
995 if( p .le. tol ) then
998 if (p .lt. dabs(a(i,m))) goto 21
1004 p = dsqrt( a(k,m)**2 + a(l,m)**2 )
1005 if( p .le. tol ) goto 40
1016 a(1,2) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
1017 a(2,2) = a(3,3)*a(1,1) - a(3,1)*a(1,3)
1018 a(3,2) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
1023 b(i,l) = r(i,1)*a(1,l) + r(i,2)*a(2,l) + r(i,3)*a(3,l)
1026 if( d .gt. zero ) d = 1.0 / dsqrt(d)
1031 d = b(1,1)*b(1,2) + b(2,1)*b(2,2) + b(3,1)*b(3,2)
1035 b(i,2) = b(i,2) - d*b(i,1)
1038 if( p .le. tol ) then
1041 if(p.lt.dabs(b(i,1)))goto 22
1047 p = dsqrt( b(k,1)**2 + b(l,1)**2 )
1048 if( p .le. tol ) goto 40
1059 b(1,3) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
1060 b(2,3) = b(3,1)*b(1,2) - b(3,2)*b(1,1)
1061 b(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
1065 u(i,j) = b(i,1)*a(j,1) + b(i,2)*a(j,2) + b(i,3)*a(j,3)
1070 t(i) = ((yc(i) - u(i,1)*xc(1)) - u(i,2)*xc(2)) - u(i,3)*xc(3)
1073 if( e(i) .lt. zero ) e(i) = zero
1074 e(i) = dsqrt( e(i) )
1078 if( e(2) .le. (e(1) * 1.0d-05) ) ier = -1
1081 if( sigma .lt. 0.0 ) then
1083 if( (e(2) - e(3)) .le. (e(1) * 1.0d-05) ) ier = -1
1085 d = (d + e(2)) + e(1)
1087 if(mode .eq. 2.or.mode.eq.0) then ! need rmsd
1089 if( rms .lt. 0.0 ) rms = 0.0