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