C *** B M F I T *** BEST MOLECULAR FIT *** S. C. NYBURG ** C ************************* UPDATE OCT 1981 *********** IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 TYPE CHARACTER*8 PTYPE INTEGER*2 ORD COMMON TYPE(2,99),N,J,DS,NA,NB COMMON D(99), U(3,3),W(99),X(2,99),Y(2,99),Z(2,99),PTYPE(99) DIMENSION CQ(3,3),A(9),R(9),ORD(99),XP(99),YP(99),ZP(99),EPX(99) DIMENSION CA(2),CB(2),CG(2),SG(2),CR(2),CP(2),EX(2,99) DIMENSION EY(2,99),EZ(2,99),FRX(2,99),FRY(2,99),FRZ(2,99),RLC(2,3) DIMENSION XS(2,99),YS(2,99),ZS(2,99),RAN(2,3),EPY(99),EPZ(99) DIMENSION TITLE(18) OPEN(1, FILE=' ') OPEN(2, FILE=' ') OPEN(7, FILE=' ') OPEN(8, FILE=' ') 1000 READ(7,79,END=999)TITLE 79 FORMAT(18A4) DO 100 NY=1,2 DO 100 NZ=1,50 X(NY,NZ)=0. Y(NY,NZ)=0. 100 Z(NY,NZ)=0. DS=0. J=0 C ITAP IS TH INPUT FILE HNIT FOR MOLECULE A IF=1,FOR B IF = 2 READ(7,*)N,NA,NB,NFIX,NFIT,NW,ITEMP,NDISCA,NDISCB,NCUTA,NCUTB,NNN C** 1 FORMAT(12I5) IF(NA-NB)101,101,103 101 NTOT=NB GO TO 10 103 NTOT=NA 10 CONTINUE DPR=57.2957795 DO 33 K=1,2 C ****** RAN IN EITHER COSINES OR DEGREES READ(7,* )RLC(K,1),RLC(K,2),RLC(K,3),RAN(K,1),RAN(K,2),RAN(K,3) C**12 FORMAT(6F10.6) IF (RAN(K,1).LT.1.0) GO TO 102 RAN(K,1)=COS(RAN(K,1)/DPR) RAN(K,2)=COS(RAN(K,2)/DPR) RAN(K,3)=COS(RAN(K,3)/DPR) 102 CA(K) = RAN(K,1) CB(K) = RAN(K,2) CG(K) = RAN(K,3) SG(K) = SQRT(1.0 - RAN(K,3)**2) CR(K)=(CA(K)-CB(K)*CG(K))/SG(K) CP(K)=SQRT(1.-CB(K)*CB(K)-CR(K)*CR(K)) IF(K-1)22,35,22 35 NOW=NA NOWW=NDISCA ITAP=1 IF(NCUTA.EQ.0)GO TO 2333 II=0 2033 II=II+1 READ(ITAP,*) C2133 FORMAT(A4) IF(II.EQ.NCUTA)GO TO 2333 GO TO 2033 2333 CONTINUE GO TO 25 22 NOW=NB NOWW=NDISCB ITAP=2 IF(NCUTB.EQ.0)GO TO 3333 II=0 3033 II=II+1 READ(ITAP,*) C3133 FORMAT(A4) IF(II.EQ.NCUTB)GO TO 3333 GO TO 3033 3333 CONTINUE 25 DO 33 KJ=1,NOW IF(KJ.GT.NOWW)GO TO 4033 READ(ITAP,*)TYPE(K,KJ),FRX(K,KJ),FRY(K,KJ),FRZ(K,KJ) 3 FORMAT (A6,21X,3F9.6 ) IF(ITEMP.NE.0) READ(ITAP,*) GO TO 4032 4033 READ(7,*)TYPE(K,KJ),FRX(K,KJ),FRY(K,KJ),FRZ(K,KJ) IF(ITEMP.NE.0) READ(7,3) 4032 CONTINUE X(K,KJ)=FRX(K,KJ)*RLC(K,1)+FRY(K,KJ)*RLC(K,2)*CG(K)+FRZ(K,KJ)* 1RLC(K,3)*CB(K) Y(K,KJ)=FRY(K,KJ)*RLC(K,2)*SG(K)+FRZ(K,KJ)*RLC(K,3)*CR(K) Z(K,KJ)=FRZ(K,KJ)*RLC(K,3)*CP(K) IF(NW)44,44,45 44 W(KJ)=1. GO TO 33 45 READ(7,* )EX(K,KJ),EY(K,KJ),EZ(K,KJ) EX(K,KJ)=EX(K,KJ)*RLC(K,1)+EY(K,KJ)*RLC(K,2)*CG(K)+EZ(K,KJ)* 1RLC(K,3)*CB(K) EY(K,KJ)=EY(K,KJ)*RLC(K,2)*SG(K)+EZ(K,KJ)*RLC(K,3)*CR(K) EZ(K,KJ)=EZ(K,KJ)*RLC(K,3)*CP(K) C**46 FORMAT( 27X,3F9.6) EX(K,KJ)=SQRT((EX(K,KJ)**2+EY(K,KJ)**2+EZ(K,KJ)**2)/3.) 33 CONTINUE IF(NNN.EQ.0)GO TO 7 READ(7,* )ORD C*105 FORMAT(24I3) DO 8 NORD=1,NB IF(ORD(NORD).EQ.0)ORD(NORD)=NORD PTYPE(NORD)=TYPE(2,ORD(NORD)) XP(NORD)=X(2,ORD(NORD)) YP(NORD)=Y(2,ORD(NORD)) ZP(NORD)=Z(2,ORD(NORD)) EPX(NORD)=EX(2,ORD(NORD)) EPY(NORD)=EY(2,ORD(NORD)) EPZ(NORD)=EZ(2,ORD(NORD)) 8 CONTINUE DO 9 M=1,NB TYPE(2,M)=PTYPE(M) X(2,M)=XP(M) Y(2,M)=YP(M) Z(2,M)=ZP(M) EX(2,M)=EPX(M) EY(2,M)=EPY(M) EZ(2,M)=EPZ(M) 9 CONTINUE C FIND WEIGHT FOR EACH MATCHED ATOM PAIR 7 IF(NW.EQ.0) GO TO 49 SES=0. DO 47 KJ=1,N EX(1,KJ)=SQRT((EX(1,KJ))**2+(EX(2,KJ))**2) 47 SES=SES+ (1./(EX(1,KJ)**2)) AN=N AN=AN/SES DO 48 KJ=1,N 48 W(KJ)=AN/(EX(1,KJ)**2) 49 WRITE(8,179)TITLE 179 FORMAT(1H1,20A4,//) IF(NW.EQ.0) GO TO 75 WRITE(8,76) 76 FORMAT(31X,'THIS IS A WEIGHTED FIT'//) GO TO 80 75 WRITE(8,77) 77 FORMAT(30X,'THIS IS AN UNWEIGHTED FIT'//) 80 WRITE(8,23) 23 FORMAT(72H ORTHOGONAL COORDINATES OF ATOMS 1 BEFORE MATCHING//) WRITE(8,24) 24 FORMAT(21X,'FIRST MOLECULE',31X,'SECOND MOLECULE'//) WRITE(8,20)RLC(1,1),RLC(1,2),RLC(1,3),RLC(2,1),RLC(2,2),RLC(2,3) 20 FORMAT(3X,'A=',F9.4,' B=',F9.4,' C=',F9.4,10X,'A=',F9.4,' B= 1',F9.4,' C=',F9.4) WRITE(8,43)RAN(1,1),RAN(1,2),RAN(1,3),RAN(2,1),RAN(2,2),RAN(2,3) 43 FORMAT(3X,'COSA=',F9.6,',COSB=',F9.6,',COSG=',F9.6,5X,'COSA=',F9.6 1,',COSB=',F9.6,',COSG=',F9.6//) CALL NATOMS C TYPE OF FORMAT 60 IF((NA-N).EQ.0)GO TO 50 IF((NB-N).EQ.0)GO TO 54 IF(NA-NB)51,52,53 C OUTPUT ROUTINE F5 51 NP=N+1 NAP=NA+1 IF(J.NE.0 )GO TO 84 WRITE(8,26)(MI,TYPE(1,MI),X(1,MI),Y(1,MI),Z(1,MI),TYPE(2,MI),X(2,M 1I),Y(2,MI),Z(2,MI),MI=NP,NA) 26 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6) WRITE(8,74)(MP,TYPE(2,MP),X(2,MP),Y(2,MP),Z(2,MP),MP=NAP,NB) 74 FORMAT(I4,39X,A6,3F10.6) GO TO 34 84 WRITE(8,82)(MM,TYPE(1,MM),X(1,MM),Y(1,MM),Z(1,MM),TYPE(2,MM),X(2,M 1M),Y(2,MM),Z(2,MM),D(MM),MM=NP,NA) 82 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) 85 WRITE(8,42)(MR,TYPE(2,MR),X(2,MR),Y(2,MR),Z(2,MR),MR=NAP,NB) 42 FORMAT(I4,39X,A6,3F10.6) 34 CONTINUE GO TO 55 C OUTPUT ROUTINE F6 52 NP=N+1 IF(J.NE.0 )GO TO 86 WRITE(8,27)(MS,TYPE(1,MS),X(1,MS),Y(1,MS),Z(1,MS),TYPE(2,MS),X(2,M 1S),Y(2,MS),Z(2,MS),MS=NP,NA) 27 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6) GO TO 87 86 WRITE(8,36)(MT,TYPE(1,MT),X(1,MT),Y(1,MT),Z(1,MT),TYPE(2,MT),X(2,M 1T),Y(2,MT),Z(2,MT),D(MT),MT=NP,NA) 36 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) 87 CONTINUE GO TO 55 C OUTPUT ROUTINE F4 53 NP=N+1 NBP=NB+1 IF(J.NE.0 )GO TO 81 WRITE(8,37)(MG,TYPE(1,MG),X(1,MG),Y(1,MG),Z(1,MG),TYPE(2,MG),X(2,M 1G),Y(2,MG),Z(2,MG),MG=NP,NB) 37 FORMAT(I4,A6,3F10.6,3X,A6,3F10.6) WRITE(8,58)(MH,TYPE(1,MH),X(1,MH),Y(1,MH),Z(1,MH),MH=NBP,NA) GO TO 135 81 WRITE(8,39)(MN,TYPE(1,MN),X(1,MN),Y(1,MN),Z(1,MN),TYPE(2,MN),X(2,M 1N),Y(2,MN),Z(2,MN),D(MN),MN=NP,NB) 39 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) 83 WRITE(8,58)(ML,TYPE(1,ML),X(1,ML),Y(1,ML),Z(1,ML),ML=NBP,NA) 58 FORMAT(I4,A6,3F10.6) 135 CONTINUE GO TO 55 50 IF(NB-N)66,65,66 65 WRITE(8,29) 29 FORMAT(33H NONE//) GO TO 55 C OUTPUT ROUTINE F3 66 NBP=N+1 WRITE(8,40)(MF,TYPE(2,MF),X(2,MF),Y(2,MF),Z(2,MF),MF=NBP,NB) 40 FORMAT(I4,39X,A6,3F10.6) GO TO 55 C OUTPUT ROUTINE F2 54 NAP=N+1 WRITE(8,41)(ME,TYPE(1,ME),X(1,ME),Y(1,ME),Z(1,ME),ME=NAP,NA) 41 FORMAT(I4,A6,3F10.6) 55 IF(J.NE.0)GO TO 61 WRITE(8,90) 90 FORMAT(1H0) WRITE(8,93) 93 FORMAT(69H CENTROID 1 CENTROID/) C FIND PAIR-WEIGHTED CENTROIDS OF MOLECULES IF(NFIX.NE.0)GO TO 62 XA=0. YA=0. ZA=0. XB=0. YB=0. ZB=0. DO 4 NF=1,N XA=XA+X(1,NF)*W(NF) YA=YA+Y(1,NF)*W(NF) ZA=ZA+Z(1,NF)*W(NF) XB=XB+X(2,NF)*W(NF) YB=YB+Y(2,NF)*W(NF) 4 ZB=ZB+Z(2,NF)*W(NF) AN=N XA=XA/AN YA=YA/AN ZA=ZA/AN XB=XB/AN YB=YB/AN ZB=ZB/AN GO TO 63 62 XA=X(1,NFIX) YA=Y(1,NFIX) ZA=Z(1,NFIX) XB=X(2,NFIX) YB=Y(2,NFIX) ZB=Z(2,NFIX) 63 WRITE(8,94)XA,YA,ZA,XB,YB,ZB 94 FORMAT(10X,3F10.6,9X,3F10.6,//) WRITE(8,78) 78 FORMAT(1H1) IF(NW.EQ.0) GO TO 175 WRITE(8,76) GO TO 180 175 WRITE(8,77) C MOVE CENTROIDS TO ORIGIN 180 DO 5 NG=1,NTOT X(1,NG)=X(1,NG)-XA Y(1,NG)=Y(1,NG)-YA Z(1,NG)=Z(1,NG)-ZA X(2,NG)=X(2,NG)-XB Y(2,NG)=Y(2,NG)-YB Z(2,NG)=Z(2,NG)-ZB XS(2,NG)=X(2,NG) YS(2,NG)=Y(2,NG) 5 ZS(2,NG)=Z(2,NG) DO 997 I=1,3 DO 997 JJJ=1,3 U(I,JJJ)=0. 997 CQ(I,JJJ)=0. DO 17 I=1,N CQ(1,1)=CQ(1,1)+X(1,I )*X(2,I )*W(I) CQ(1,2)=CQ(1,2)+Y(1,I )*X(2,I )*W(I) CQ(1,3)=CQ(1,3)+Z(1,I )*X(2,I )*W(I) CQ(2,1)=CQ(2,1)+X(1,I )*Y(2,I )*W(I) CQ(2,2)=CQ(2,2)+Y(1,I )*Y(2,I )*W(I) CQ(2,3)=CQ(2,3)+Z(1,I )*Y(2,I )*W(I) CQ(3,1)=CQ(3,1)+X(1,I )*Z(2,I )*W(I) CQ(3,2)=CQ(3,2)+Y(1,I )*Z(2,I )*W(I) CQ(3,3)=CQ(3,3)+Z(1,I )*Z(2,I )*W(I) 17 CONTINUE A(1) =CQ(1,1)*CQ(1,1)+CQ(1,2)*CQ(1,2)+CQ(1,3)*CQ(1,3) A(2) =CQ(1,1)*CQ(2,1)+CQ(1,2)*CQ(2,2)+CQ(1,3)*CQ(2,3) A(4) =CQ(1,1)*CQ(3,1)+CQ(1,2)*CQ(3,2)+CQ(1,3)*CQ(3,3) A(3) =CQ(2,1)*CQ(2,1)+CQ(2,2)*CQ(2,2)+CQ(2,3)*CQ(2,3) A(5) =CQ(2,1)*CQ(3,1)+CQ(2,2)*CQ(3,2)+CQ(2,3)*CQ(3,3) A(6) =CQ(3,1)*CQ(3,1)+CQ(3,2)*CQ(3,2)+CQ(3,3)*CQ(3,3) J=1 CALL EIGEN(A,R) NWR=0 NMIR=0 CALL VECTOR(A,R,CQ,NWR,NMIR,DETA) 213 IF(NWR.NE.0) GO TO 204 IF(NMIR.EQ.1) GO TO 209 IF(DETA.LT.0.) GO TO 203 WRITE(8,211) 211 FORMAT(20X,'THE BEST FIT IS GIVEN WITH THE SECOND MOLECULE UN-MIRR 1ORED',///) GO TO 207 203 WRITE(8,212) 212 FORMAT(20X,'THE BEST FIT IS GIVEN WITH THE SECOND MOLECULE MIRRORE 1D',///) GO TO 207 204 IF(DDET.LT.0.) GO TO 205 WRITE(8,206) 206 FORMAT(1H1,35X,'THIS IS THE MIRRORED FIT',///) GO TO 2 205 WRITE(8,208) 208 FORMAT(1H1,20X,'THIS IS THE FIT WITH THE SECOND MOLECULE UN-MIRROR 1ED',///) GO TO 2 209 WRITE(8,210) 210 FORMAT( 1X,'EITHER ONE OF THE MOLECULES IS FLAT OR THE FIT IS THE 1SAME WHETHER OR NOT THE 2ND MOLECULE IS MIRRORED',///) 207 DDET=DETA 2 DO 994 I=1,NB XN=A(1)*XS(2,I)+A(2)*YS(2,I)+A(3)*ZS(2,I) YN=A(4)*XS(2,I)+A(5)*YS(2,I)+A(6)*ZS(2,I) ZN=A(7)*XS(2,I)+A(8)*YS(2,I)+A(9)*ZS(2,I) X(2,I)=XN Y(2,I)=YN Z(2,I)=ZN 994 CONTINUE CALL DSUM WRITE(8,11) 11 FORMAT(87H 1 DELTA/) CALL NATOMS GO TO 60 61 CONTINUE C THIS FINAL SECTION PERFORMS THE TRANSFORMATION OF 2ND MOLECULE INT C THE FRACTIONAL COORDINATES OF THE FIRST. C IF THIS IS REQUIRED, PUT NFIT=1. OTHERWISE NFIT=0. IF(NFIT.EQ.0) GO TO 161 C BRING 2ND MOLECULE BACK TO SAME CRYSTAL FRAME AS 1ST MOLECULE C MOVE 2ND MOLECULE TO CENTROID OF 1ST MOLECULE IN ITS ORIGINAL POSI DO 16 NG=1,NB X(2,NG)=X(2,NG)+XA Y(2,NG)=Y(2,NG)+YA Z(2,NG)=Z(2,NG)+ZA 16 CONTINUE C INVERSE MATRIX TO BRING TO FRACTIONAL COORDS DO 30 K=1,NB X(2,K)=(X(2,K)-(Y(2,K)*CG(1)/SG(1)) 1 +(Z(2,K)/CP(1))*((CR(1)*CG(1)/SG(1)) -CB(1)))/RLC(1,1) Y(2,K)= (Y(2,K)/SG(1)-((Z(2,K)*CR(1))/(CP(1)*SG(1))))/RLC(1,2) 30 Z(2,K)=Z(2,K)/(CP(1)*RLC(1,3)) WRITE(8,250) 250 FORMAT ( 1H0,///,' NO', 4X,' ATOM ',6X,'FRACT. COORDINATES', 1 11X,' ATOM FRACT. COORDS. OF 2ND MOLECULE',/23X,'FIRST MOLECULE 2',26X,'IN CELL OF FIRST') IF(NA-NB)67,67,68 67 NUM=NA GO TO 70 68 NUM=NB 70 WRITE(8,59)(M,TYPE(1,M),FRX(1,M),FRY(1,M),FRZ(1,M),TYPE(2,M),X(2,M 1),Y(2,M),Z(2,M),M=1,NUM) 59 FORMAT(I4,5X,A6,3F10.6,5X,A6,3F10.6) IF(NA-NB)71,161,72 71 NAP=NA+1 WRITE(8,99)(M,TYPE(2,M),X(2,M),Y(2,M),Z(2,M),M=NAP,NB) 99 FORMAT(I4,46X,A6,3F10.6) GO TO 161 72 NBP=NB+1 WRITE(8,73)(M,TYPE(1,M),FRX(1,M),FRY(1,M),FRZ(1,M),M=NBP,NA) 73 FORMAT(I4,5X,A6,3F10.6) 161 IF(NMIR.EQ.1) GO TO 1000 NWR=NWR+1 IF(NWR.EQ.2)GO TO 1000 CALL VECTOR(A,R,CQ,NWR,NMIR,DETA) GO TO 204 999 STOP END SUBROUTINE EIGEN(A,R) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(6),R(9) RANGE=1.0E-6 IQ=-3 DO 20 J=1,3 IQ=IQ+3 DO 20 I=1,3 IJ=IQ+I R(IJ)=0.0 IF(I-J) 20,15,20 15 R(IJ)=1.0 20 CONTINUE 25 ANORM=0.0 DO 35 I=1,3 DO 35 J=I,3 IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM)165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*RANGE/3. IND=0 THR=ANORM 45 THR=THR/3. 50 L=1 55 M=L+1 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(ABS(A(LM))-THR) 130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X =0.5*(A(LL)-A(MM)) 68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y)))) SINX2=SINX*SINX 78 COSX=SQRT(1.0-SINX2) COSX2=COSX*COSX SINCS=SINX*COSX ILQ=3*(L-1) IMQ=3*(M-1) DO 125 I=1,3 IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GO TO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GO TO 110 105 IL=L+IQ 110 X =A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 ILR=ILQ+I IMR=IMQ+I X =R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X =2.0*A(LM)*SINCS Y =A(LL)*COSX2+A(MM)*SINX2-X X =A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X 130 IF(M-3) 135,140,135 135 M=M+1 GO TO 60 140 IF(L-2) 145,150,145 145 L=L+1 GO TO 55 150 IF(IND-1) 160,155,160 155 IND=0 GO TO 50 160 IF(THR-ANRMX) 165,165,45 165 IQ=-3 DO 185 I=1,3 IQ=IQ+3 LL=I+(I*I-I)/2 JQ=3*(I-2) DO 185 J=I,3 JQ=JQ+3 MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X =A(LL) A(LL)=A(MM) A(MM)=X DO 180 K=1,3 ILR=IQ+K IMR=JQ+K X =R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE RETURN END SUBROUTINE VECTOR(A,R,CQ,NWR,NMIR,DETA) IMPLICIT REAL*8 (A-H, O-Z) DIMENSION CQ(3,3) DIMENSION A(9),R(9) DIMENSION U(3,3) DIMENSION WV(3),ZZ(3,3),RC(3,3),B(3,3) DO 28 I=1,3 DO 28 J=1,3 28 U(I,J)=0. IF(NWR.NE.0) GO TO 52 WV(1)=A(1) WV(2)=A(3) WV(3)=A(6) ZZ(1,1)=R(1) ZZ(1,2)=R(2) ZZ(1,3)=R(3) ZZ(2,1)=R(4) ZZ(2,2)=R(5) ZZ(2,3)=R(6) ZZ(3,1)=R(7) ZZ(3,2)=R(8) ZZ(3,3)=R(9) DO 21 I=1,3 DO 21 J=1,3 B(I,J)=0. 21 RC(I,J)=CQ(J,I) C MAKE SURE DET OF RC IS NOT EXACTLY ZERO 999 FORMAT(/' VALUE OF DETERMINANT',E14.4) DETR= RC(1,1)*RC(2,2)*RC(3,3) $ -RC(1,1)*RC(2,3)*RC(3,2) $ +RC(1,2)*RC(2,3)*RC(3,1) $ -RC(1,2)*RC(2,1)*RC(3,3) $ +RC(1,3)*RC(2,1)*RC(3,2) $ -RC(1,3)*RC(2,2)*RC(3,1) WRITE(8,999)DETR IF(ABS(DETR).LT.0.0001 )GO TO 17 DO 19 K=1,3 DO 19 I=1,3 DO 19 J=1,3 B(K,I)=B(K,I)+RC(I,J)*ZZ(K,J) 19 CONTINUE GO TO 38 17 NMIR=1 IF(WV(2).LT.1.0E-6) GO TO 39 DO 40 K=1,2 DO 40 I=1,3 DO 40 J=1,3 B(K,I)=B(K,I)+RC(I,J)*ZZ(K,J) 40 CONTINUE B (3,1)=B (1,2)*B (2,3)-B (1,3)*B (2,2) B (3,2)=B (1,3)*B (2,1)-B (2,3)*B (1,1) B (3,3)=B (1,1)*B (2,2)-B (2,1)*B (1,2) 38 BS= B(1,1)**2+B(1,2)**2+B(1,3)**2 IF (BS.EQ.0.0) GO TO 100 FACS=1./(B(1,1)**2+B(1,2)**2+B(1,3)**2) DO 13 J=1,3 13 B(1,J)=SQRT(FACS)*B(1,J) 100 BS= B(2,1)**2+B(2,2)**2+B(2,3)**2 IF(BS.EQ.0.0) GO TO 101 FACS=1./(B(2,1)**2+B(2,2)**2+B(2,3)**2) DO 14 J=1,3 14 B(2,J)=SQRT(FACS)*B(2,J) 101 BS= B(3,1)**2+B(3,2)**2+B(3,3)**2 IF(BS.EQ.0.0) GO TO 102 FACS=1./(B(3,1)**2+B(3,2)**2+B(3,3)**2) DO 15 J=1,3 15 B(3,J)=SQRT(FACS)*B(3,J) GO TO 102 39 DO 41 I=1,3 DO 41 J=1,3 B(1,I)=B(1,I)+RC(I,J)*ZZ(1,J) 41 CONTINUE FACS=1./(B(1,1)**2+B(1,2)**2+B(1,3)**2) DO 42 J=1,3 B(1,J)=SQRT(FACS)*B(1,J) 42 CONTINUE IF(B(1,1).EQ.0.) GO TO 43 IF(B(1,2).EQ.0.) GO TO 37 B(2,2)=B(1,1)/SQRT(B(1,1)**2+B(1,2)**2) B(2,1)=-(B(1,2)*B(2,2))/B(1,1) B(2,3)=0. GO TO 44 43 B(2,1)=1. B(2,2)=0. B(2,3)=0. GO TO 44 37 B(2,1)=0. B(2,2)=1. B(2,3)=0. 44 B(3,1)=B(1,2)*B(2,3)-B(1,3)*B(2,2) B(3,2)=B(1,3)*B(2,1)-B(1,1)*B(2,3) B(3,3)=B(1,1)*B(2,2)-B(1,2)*B(2,1) GO TO 102 52 B(3,1)=-B(3,1) B(3,2)=-B(3,2) B(3,3)=-B(3,3) 102 DO 20 I=1,3 DO 20 J=1,3 DO 20 K=1,3 20 U(I,J)=U(I,J)+B(K,I)*ZZ(K,J) A(1)=U(1,1) A(2)=U(1,2) A(3)=U(1,3) A(4)=U(2,1) A(5)=U(2,2) A(6)=U(2,3) A(7)=U(3,1) A(8)=U(3,2) A(9)=U(3,3) C IS A MATRIX ORTHOGONAL? A11=A(1)*A(1)+A(2)*A(2)+A(3)*A(3) A12=A(1)*A(4)+A(2)*A(5)+A(3)*A(6) A13=A(1)*A(7)+A(2)*A(8)+A(3)*A(9) A21=A(4)*A(1)+A(5)*A(2)+A(6)*A(3) A22=A(4)*A(4)+A(5)*A(5)+A(6)*A(6) A23=A(4)*A(7)+A(5)*A(8)+A(6)*A(9) A31=A(7)*A(1)+A(8)*A(2)+A(9)*A(3) A32=A(7)*A(4)+A(8)*A(5)+A(9)*A(6) A33=A(7)*A(7)+A(8)*A(8)+A(9)*A(9) WRITE(8,998)A11,A12,A13 WRITE(8,998)A21,A22,A23 WRITE(8,998)A31,A32,A33 998 FORMAT(/20X,3E14.4,/) DETA=A(1)*A(5)*A(9)-A(1)*A(6)*A(8)+A(2)*A(6)*A(7)-A(2)*A(4)*A(9)+A 1(3)*A(4)*A(8)-A(3)*A(5)*A(7) RETURN END SUBROUTINE DSUM IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 TYPE COMMON TYPE(2,99),N,J,DS,NA,NB COMMON D(99), U(3,3),W(99),X(2,99),Y(2,99),Z(2,99) DIMENSION A(9),R(9) DS=0. IF(NA-NB)96,96,98 96 NSUM=NA GO TO 99 98 NSUM=NB 99 DO 60 MB=1,NSUM 60 D(MB)= (X(1,MB)-X(2,MB))**2+(Y(1,MB)-Y(2,MB))**2+(Z(1,MB 1)-Z(2,MB))**2 DO 160 MU=1,N 160 DS=DS+D(MU) DO 161 MV=1,NSUM 161 D(MV)=SQRT(D(MV)) RETURN END SUBROUTINE NATOMS IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 TYPE COMMON TYPE(2,99),N,J,DS,NA,NB COMMON D(99), U(3,3),W(99),X(2,99),Y(2,99),Z(2,99) IF(J.NE.0 )GO TO 56 WRITE(8,26)(MA,TYPE(1,MA),X(1,MA),Y(1,MA),Z(1,MA),TYPE(2,MA),X(2,M 1A),Y(2,MA),Z(2,MA),MA=1,N) 26 FORMAT(I4,A6,3F10.6,12X,A6,3F10.6) GO TO 57 56 WRITE(8,27)(MD,TYPE(1,MD),X(1,MD),Y(1,MD),Z(1,MD),TYPE(2,MD),X(2,M 1D),Y(2,MD),Z(2,MD),D(MD),MD=1,N) 27 FORMAT(I4,A6,3F10.6,3X,A6,4F10.6) WRITE(8,164) DS 164 FORMAT(/,' SUM OF DIFFERENCES SQUAR 1ED ',F10.6,//) 57 WRITE (8,29) 29 FORMAT(2H0 ,///) WRITE (8,28) 28 FORMAT(56H UNMATCHED ATOMS 1//) RETURN END