1 FUNCTION CCuij(UU,VV) 2C =============================================================== 3C Correlation coefficient between two (3x3) anisotropic 4C displacement matrices Uij and Vij. 5C See Merritt (1999) Acta Crystallographica D55, 1997-2004. 6C Return 0 if any of the tensors involved are NPD 7C =============================================================== 8 REAL CCuij V3VINV 9 REAL UU(3,3), VV(3,3), WW(3,3), UI(3,3), VI(3,3), WI(3,3) 10 REAL DUI, DVI, DWI 11 DUI = 1. / V3INV(UI,UU) 12 DVI = 1. / V3INV(VI,VV) 13C Sylvester's criterion for NPD tensor 14 IF (DUI.LE.0 .OR. UU(1,1).LE.0 15 & .OR. (UU(1,1)*UU(2,2)-UU(1,2)*UU(2,1)).LE.0) THEN 16 CCuij = 0 17 RETURN 18 ENDIF 19 IF (DVI.LE.0 .OR. VV(1,1).LE.0 20 & .OR. (VV(1,1)*VV(2,2)-VV(1,2)*VV(2,1)).LE.0) THEN 21 CCuij = 0 22 RETURN 23 ENDIF 24 CALL M3ADD(WI,UI,VI) 25 DWI = V3INV(WW,WI) 26 IF (DWI.LE.0 .OR. WW(1,1).LE.0 27 & .OR. (WW(1,1)*WW(2,2)-WW(1,2)*WW(2,1)).LE.0) THEN 28 CCuij = 0 29 RETURN 30 ENDIF 31 CCuij = sqrt(sqrt(DUI*DVI)) 32 & / sqrt(0.125*DWI) 33 RETURN 34 END 35 36 37 FUNCTION CCuv(U,V) 38C =============================================================== 39C Wrapper for CCuij() 40C =============================================================== 41 REAL CCuv, CCuij 42 REAL U(6), V(6) 43 REAL UU(3,3), VV(3,3) 44 UU(1,1) = U(1) 45 UU(2,2) = U(2) 46 UU(3,3) = U(3) 47 UU(1,2) = U(4) 48 UU(2,1) = U(4) 49 UU(1,3) = U(5) 50 UU(3,1) = U(5) 51 UU(2,3) = U(6) 52 UU(3,2) = U(6) 53 VV(1,1) = V(1) 54 VV(2,2) = V(2) 55 VV(3,3) = V(3) 56 VV(1,2) = V(4) 57 VV(2,1) = V(4) 58 VV(1,3) = V(5) 59 VV(3,1) = V(5) 60 VV(2,3) = V(6) 61 VV(3,2) = V(6) 62 CCuv = CCuij(UU,VV) 63 RETURN 64 END 65 66 FUNCTION Suv(U,V) 67C =============================================================== 68C Normalized correlation coefficient ("similarity") between two 69C anisotropic displacement vectors Uij and Vij. 70C See Merritt (1999) Acta Crystallographica D55, 1997-2004. 71C =============================================================== 72 REAL Suv, CCuij 73 REAL U(6), V(6) 74 REAL UU(3,3), VV(3,3), Uiso(3,3), Viso(3,3), WW(3,3) 75 REAL Ueq, Veq 76 UU(1,1) = U(1) 77 UU(2,2) = U(2) 78 UU(3,3) = U(3) 79 UU(1,2) = U(4) 80 UU(2,1) = U(4) 81 UU(1,3) = U(5) 82 UU(3,1) = U(5) 83 UU(2,3) = U(6) 84 UU(3,2) = U(6) 85 VV(1,1) = V(1) 86 VV(2,2) = V(2) 87 VV(3,3) = V(3) 88 VV(1,2) = V(4) 89 VV(2,1) = V(4) 90 VV(1,3) = V(5) 91 VV(3,1) = V(5) 92 VV(2,3) = V(6) 93 VV(3,2) = V(6) 94 Ueq = (UU(1,1) + UU(2,2) + UU(3,3)) / 3. 95 Veq = (VV(1,1) + VV(2,2) + VV(3,3)) / 3. 96 DO I = 1,3 97 DO J = 1,3 98 WW(I,J) = VV(I,J) * Ueq / Veq 99 Uiso(I,J) = 0 100 Viso(I,J) = 0 101 ENDDO 102 ENDDO 103 DO I = 1,3 104 Uiso(I,I) = Ueq 105 Viso(I,I) = Veq 106 ENDDO 107 Suv_top = CCuij(UU,WW) 108 Suv_bot = (CCuij(UU,Uiso) * CCuij(VV,Viso)) 109 if (Suv_top.eq.0 .or. Suv_bot.eq.0) then 110 Suv = 0 111 else 112 Suv = Suv_top / Suv_bot 113 endif 114 RETURN 115 END 116 117C ======================================================= 118C The rest of the file is just generic 3x3 matrix algebra 119C ======================================================= 120 121 SUBROUTINE V3CROSS(B,C,A) 122C CROSS PRODUCT OF TWO VECTORS 123C ======================= 124 REAL A(3),B(3),C(3) 125 A(1)=B(2)*C(3)-C(2)*B(3) 126 A(2)=B(3)*C(1)-C(3)*B(1) 127 A(3)=B(1)*C(2)-C(1)*B(2) 128 RETURN 129 END 130C 131 FUNCTION V3DOT(A,B) 132C DOT PRODUCT OF TWO VECTORS 133C ================ 134 REAL A(3),B(3) 135 V3DOT=A(1)*B(1)+A(2)*B(2)+A(3)*B(3) 136 RETURN 137 END 138C 139 SUBROUTINE M3ADD(A,B,C) 140C A=B+C 141C ======================== 142 REAL A(3,3),B(3,3),C(3,3) 143 DO I=1,3 144 DO J=1,3 145 A(I,J) = B(I,J) + C(I,J) 146 ENDDO 147 ENDDO 148 RETURN 149 END 150C 151 SUBROUTINE M3MUL(A,B,C) 152C A=B*C 153C ======================== 154 REAL A(3,3),B(3,3),C(3,3),S 155 DO I=1,3 156 DO J=1,3 157 S=0 158 DO K=1,3 159 S=S+B(I,K)*C(K,J) 160 ENDDO 161 A(I,J)=S 162 ENDDO 163 ENDDO 164 RETURN 165 END 166C 167 FUNCTION V3INV(A,B) 168C INVERT A GENERAL 3X3 MATRIX AND RETURN DETERMINANT 169C A=(B)-1 170C ====================== 171 REAL A(3,3),B(3,3),C(3,3),D 172 CALL V3CROSS(B(1,2),B(1,3),C(1,1)) 173 CALL V3CROSS(B(1,3),B(1,1),C(1,2)) 174 CALL V3CROSS(B(1,1),B(1,2),C(1,3)) 175 D=V3DOT(B(1,1),C(1,1)) 176 IF (ABS(D).LE.1.E-30) THEN 177 V3INV=0.0 178 RETURN 179 ENDIF 180 DO I=1,3 181 DO J=1,3 182 A(I,J)=C(J,I)/D 183 ENDDO 184 ENDDO 185 V3INV=D 186 RETURN 187 END 188C 189