Line data Source code
1 : !!****m* ABINIT/m_contract
2 : !! NAME
3 : !! m_contract
4 : !!
5 : !! FUNCTION
6 : !! Low-level procedeures used in nonlop_pl to contract tensors
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, MT, GZ)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_contract
23 :
24 : use defs_basis
25 : use m_errors
26 :
27 : implicit none
28 :
29 : private
30 : !!***
31 :
32 : public :: cont22cso
33 : public :: cont22so
34 : public :: cont24
35 : public :: cont33cso
36 : public :: cont33so
37 : public :: cont35
38 : public :: cont22
39 : public :: cont3
40 : public :: cont13
41 : public :: metcon
42 : public :: metcon_so
43 : public :: metric_so
44 : !!***
45 :
46 : contains
47 : !!***
48 :
49 : !!****f* m_contract/cont13
50 : !! NAME
51 : !! cont13
52 : !!
53 : !! FUNCTION
54 : !! Contract rank1 tensor with rank3 symmetric tensor to
55 : !! produce symmetric rank2 tensor
56 : !!
57 : !! INPUTS
58 : !! rank1(2,3)=rank 1 complex tensor (vector of length 3)
59 : !! rank3(2,10)=rank 3 complex tensor (symmetric storage)
60 : !!
61 : !! OUTPUT
62 : !! rank2(6)=rank 2 real tensor (symmetric storage)
63 : !!
64 : !! NOTES
65 : !! Tensors are in "symmetric" storage mode.
66 : !! For rank1 this is 1, 2, 3;
67 : !! for rank2 this is 11, 22, 33, 32, 31, 21;
68 : !! for rank3 this is 111, 221, 331, 321, 311, 211, 222, 332, 322, 333.
69 : !! rank1 and rank3 are complex; rank2 is real.
70 : !! Want $2 Re[contraction]$.
71 : !! $rank2(a,b)=2 Re[rank1(i)^"*" rank3(a,b,i)]$.
72 : !! In typical usage the input rank1 tensor is actually
73 : !! $rank1(i)=gmet(i,j) gxa(j)$
74 : !!
75 : !! SOURCE
76 :
77 867129 : subroutine cont13(rank1,rank3,rank2)
78 :
79 : !Arguments ------------------------------------
80 : !arrays
81 : real(dp),intent(in) :: rank1(2,3),rank3(2,10)
82 : real(dp),intent(out) :: rank2(6)
83 :
84 : !Local variables-------------------------------
85 : !scalars
86 : integer,parameter :: im=2,re=1
87 :
88 : ! *************************************************************************
89 :
90 : !Simply write out index summations
91 : !a=1, b=1 in rank2(a,b) --> maps to index 1
92 : rank2(1)=2.0d0*(&
93 : & (rank1(re,1)*rank3(re,1)+rank1(im,1)*rank3(im,1))+&
94 : & (rank1(re,2)*rank3(re,6)+rank1(im,2)*rank3(im,6))+&
95 867129 : & (rank1(re,3)*rank3(re,5)+rank1(im,3)*rank3(im,5)))
96 :
97 : !a=2, b=2 in rank2(a,b) --> maps to index 2
98 : rank2(2)=2.0d0*(&
99 : & (rank1(re,1)*rank3(re,2)+rank1(im,1)*rank3(im,2))+&
100 : & (rank1(re,2)*rank3(re,7)+rank1(im,2)*rank3(im,7))+&
101 867129 : & (rank1(re,3)*rank3(re,9)+rank1(im,3)*rank3(im,9)))
102 :
103 : !a=3, b=3 in rank2(a,b) --> maps to index 3
104 : rank2(3)=2.0d0*(&
105 : & (rank1(re,1)*rank3(re,3)+rank1(im,1)*rank3(im,3))+&
106 : & (rank1(re,2)*rank3(re,8)+rank1(im,2)*rank3(im,8))+&
107 867129 : & (rank1(re,3)*rank3(re,10)+rank1(im,3)*rank3(im,10)))
108 :
109 : !a=3, b=2 in rank2(a,b) --> maps to index 4
110 : rank2(4)=2.0d0*(&
111 : & (rank1(re,1)*rank3(re,4)+rank1(im,1)*rank3(im,4))+&
112 : & (rank1(re,2)*rank3(re,9)+rank1(im,2)*rank3(im,9))+&
113 867129 : & (rank1(re,3)*rank3(re,8)+rank1(im,3)*rank3(im,8)))
114 :
115 : !a=3, b=1 in rank2(a,b) --> maps to index 5
116 : rank2(5)=2.0d0*(&
117 : & (rank1(re,1)*rank3(re,5)+rank1(im,1)*rank3(im,5))+&
118 : & (rank1(re,2)*rank3(re,4)+rank1(im,2)*rank3(im,4))+&
119 867129 : & (rank1(re,3)*rank3(re,3)+rank1(im,3)*rank3(im,3)))
120 :
121 : !a=2, b=1 in rank2(a,b) --> maps to index 6
122 : rank2(6)=2.0d0*(&
123 : & (rank1(re,1)*rank3(re,6)+rank1(im,1)*rank3(im,6))+&
124 : & (rank1(re,2)*rank3(re,2)+rank1(im,2)*rank3(im,2))+&
125 867129 : & (rank1(re,3)*rank3(re,4)+rank1(im,3)*rank3(im,4)))
126 :
127 867129 : end subroutine cont13
128 : !!***
129 :
130 :
131 : !!****f* m_contract/cont22
132 : !! NAME
133 : !! cont22
134 : !!
135 : !! FUNCTION
136 : !! Contract symmetric rank 2 tensor gxa with itself using gmet to
137 : !! produce symmetric rank 2 tensor.
138 : !!
139 : !! INPUTS
140 : !! gxa(2,6)=rank 2 complex tensor
141 : !! gmet(3,3)=real symmetric metric tensor (full storage)
142 : !!
143 : !! OUTPUT
144 : !! rank2(6)=rank 2 real tensor (symmetric storage)
145 : !!
146 : !! NOTES
147 : !! Symmetric gxa is stored as 11 22 33 32 31 21;
148 : !! gmet(3,3) is symmetric but stored fully (9 elements);
149 : !! output rank2 is stored as 11 22 33 32 31 21.
150 : !! Want $2 Re[contraction]$.
151 : !! $rank2(a,b)=2 Re[gxa(i,a)^"*" gmet(i,j) gxa(j,b)]$.
152 : !!
153 : !! SOURCE
154 :
155 491121 : subroutine cont22(gxa,gmet,rank2)
156 :
157 : !Arguments ------------------------------------
158 : !arrays
159 : real(dp),intent(in) :: gmet(3,3),gxa(2,6)
160 : real(dp),intent(out) :: rank2(6)
161 :
162 : !Local variables-------------------------------
163 : !scalars
164 : integer,parameter :: im=2,re=1
165 :
166 : ! *************************************************************************
167 :
168 : !Simply write out index summations
169 : !a=1, b=1 in rank2(a,b) --> maps to index 1
170 : rank2(1)=2.0d0*(&
171 : & gmet(1,1)*(gxa(re,1)*gxa(re,1)+gxa(im,1)*gxa(im,1))+&
172 : & gmet(2,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6))+&
173 : & gmet(3,3)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+&
174 : & 2.0d0*(&
175 : & gmet(3,2)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6))+&
176 : & gmet(3,1)*(gxa(re,5)*gxa(re,1)+gxa(im,5)*gxa(im,1))+&
177 491121 : & gmet(2,1)*(gxa(re,6)*gxa(re,1)+gxa(im,6)*gxa(im,1))))
178 :
179 : !a=2, b=2 in rank2(a,b) --> maps to index 2
180 : rank2(2)=2.0d0*(&
181 : & gmet(1,1)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6))+&
182 : & gmet(2,2)*(gxa(re,2)*gxa(re,2)+gxa(im,2)*gxa(im,2))+&
183 : & gmet(3,3)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+&
184 : & 2.0d0*(&
185 : & gmet(3,2)*(gxa(re,4)*gxa(re,2)+gxa(im,4)*gxa(im,2))+&
186 : & gmet(3,1)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
187 491121 : & gmet(2,1)*(gxa(re,2)*gxa(re,6)+gxa(im,2)*gxa(im,6))))
188 :
189 : !a=3, b=3 in rank2(a,b) --> maps to index 3
190 : rank2(3)=2.0d0*(&
191 : & gmet(1,1)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+&
192 : & gmet(2,2)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+&
193 : & gmet(3,3)*(gxa(re,3)*gxa(re,3)+gxa(im,3)*gxa(im,3))+&
194 : & 2.0d0*(&
195 : & gmet(3,2)*(gxa(re,4)*gxa(re,3)+gxa(im,4)*gxa(im,3))+&
196 : & gmet(3,1)*(gxa(re,5)*gxa(re,3)+gxa(im,5)*gxa(im,3))+&
197 491121 : & gmet(2,1)*(gxa(re,5)*gxa(re,4)+gxa(im,5)*gxa(im,4))))
198 :
199 : !a=3, b=2 in rank2(a,b) --> maps to index 4
200 : rank2(4)=2.0d0*(&
201 : & gmet(1,1)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6))+&
202 : & gmet(2,2)*(gxa(re,4)*gxa(re,2)+gxa(im,4)*gxa(im,2))+&
203 : & gmet(3,3)*(gxa(re,3)*gxa(re,4)+gxa(im,3)*gxa(im,4))+&
204 : & gmet(3,2)*(gxa(re,3)*gxa(re,2)+gxa(im,3)*gxa(im,2))+&
205 : & gmet(3,1)*(gxa(re,3)*gxa(re,6)+gxa(im,3)*gxa(im,6))+&
206 : & gmet(2,1)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
207 : & gmet(2,3)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+&
208 : & gmet(1,3)*(gxa(re,5)*gxa(re,4)+gxa(im,5)*gxa(im,4))+&
209 491121 : & gmet(1,2)*(gxa(re,5)*gxa(re,2)+gxa(im,5)*gxa(im,2)))
210 :
211 : !a=3, b=1 in rank2(a,b) --> maps to index 5
212 : rank2(5)=2.0d0*(&
213 : & gmet(1,1)*(gxa(re,5)*gxa(re,1)+gxa(im,5)*gxa(im,1))+&
214 : & gmet(2,2)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
215 : & gmet(3,3)*(gxa(re,3)*gxa(re,5)+gxa(im,3)*gxa(im,5))+&
216 : & gmet(3,2)*(gxa(re,3)*gxa(re,6)+gxa(im,3)*gxa(im,6))+&
217 : & gmet(3,1)*(gxa(re,3)*gxa(re,1)+gxa(im,3)*gxa(im,1))+&
218 : & gmet(2,1)*(gxa(re,4)*gxa(re,1)+gxa(im,4)*gxa(im,1))+&
219 : & gmet(2,3)*(gxa(re,4)*gxa(re,5)+gxa(im,4)*gxa(im,5))+&
220 : & gmet(1,3)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+&
221 491121 : & gmet(1,2)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6)))
222 :
223 : !a=2, b=1 in rank2(a,b) --> maps to index 6
224 : rank2(6)=2.0d0*(&
225 : & gmet(1,1)*(gxa(re,6)*gxa(re,1)+gxa(im,6)*gxa(im,1))+&
226 : & gmet(2,2)*(gxa(re,2)*gxa(re,6)+gxa(im,2)*gxa(im,6))+&
227 : & gmet(3,3)*(gxa(re,4)*gxa(re,5)+gxa(im,4)*gxa(im,5))+&
228 : & gmet(3,2)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
229 : & gmet(3,1)*(gxa(re,4)*gxa(re,1)+gxa(im,4)*gxa(im,1))+&
230 : & gmet(2,1)*(gxa(re,2)*gxa(re,1)+gxa(im,2)*gxa(im,1))+&
231 : & gmet(2,3)*(gxa(re,2)*gxa(re,5)+gxa(im,2)*gxa(im,5))+&
232 : & gmet(1,3)*(gxa(re,6)*gxa(re,5)+gxa(im,6)*gxa(im,5))+&
233 491121 : & gmet(1,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6)))
234 :
235 491121 : end subroutine cont22
236 : !!***
237 :
238 : !!****f* m_contract/cont22cso
239 : !! NAME
240 : !! cont22cso
241 : !!
242 : !! FUNCTION
243 : !! Contract symmetric rank 2 tensor gxa1 with symmetric rank 2 tensor
244 : !! gxa2 using metric tensor gmet to produce rank 2 complex tensor.
245 : !!
246 : !! INPUTS
247 : !! gxa1(2,10)=rank 2 complex symmetric tensor
248 : !! gxa2(2,10)=rank 2 complex symmetric tensor
249 : !! gmet(3,3)=usual metric tensor (symmetric, real)
250 : !!
251 : !! OUTPUT
252 : !! rank2c(2,6)=rank 2 complex tensor (pseudo-symmetric storage)
253 : !!
254 : !! NOTES
255 : !! This contraction is used for spin-orbit correction in non-local
256 : !! contribution to stresses.
257 : !!
258 : !! Symmetric gxa1, gxa2 are stored as 11 22 33 32 31 21;
259 : !! gmet(3,3) is symmetric but stored fully (9 elements);
260 : !! Output rank2c is not symmetric but since
261 : !! $rank2c_{gxa1,gxa2}(a,b)=conjg(rank2c_{gxa2,gxa1}(b,a))$
262 : !! it is stored as 11 22 33 32 31 21.
263 : !!
264 : !! rank2c(1,1), rank2c(2,2), rank3c(3,3) are not needed;
265 : !! They are not calculated.
266 : !!
267 : !!{{\ \begin{equation}
268 : !! rank2c(a,b)=3 conjg(gxa1(i,a)) gmet(i,j) gxa2(j,b)
269 : !!\end{equation} }}
270 : !!
271 : !! SOURCE
272 :
273 285896 : subroutine cont22cso(gxa1,gxa2,gmet,rank2c)
274 :
275 : !Arguments ------------------------------------
276 : !arrays
277 : real(dp),intent(in) :: gmet(3,3),gxa1(2,6),gxa2(2,6)
278 : real(dp),intent(out) :: rank2c(2,6)
279 :
280 : !Local variables-------------------------------
281 : !scalars
282 : integer,parameter :: im=2,re=1
283 : !arrays
284 : real(dp) :: r2(2,3,3)
285 :
286 : ! *************************************************************************
287 :
288 : !Initialize output tensor
289 285896 : rank2c(:,:)=0.d0
290 :
291 : !First compute r2(i,j) = gmet(i,k) gxa2(j,k)
292 857688 : r2(:,1,1)=gmet(1,1)*gxa2(:,1)+gmet(1,2)*gxa2(:,6)+gmet(1,3)*gxa2(:,5)
293 857688 : r2(:,1,2)=gmet(1,1)*gxa2(:,6)+gmet(1,2)*gxa2(:,2)+gmet(1,3)*gxa2(:,4)
294 : r2(:,1,3)=gmet(1,1)*gxa2(:,5)+gmet(1,2)*gxa2(:,4)+gmet(1,3)*gxa2(:,3)
295 857688 : r2(:,2,1)=gmet(2,1)*gxa2(:,1)+gmet(2,2)*gxa2(:,6)+gmet(2,3)*gxa2(:,5)
296 857688 : r2(:,2,2)=gmet(2,1)*gxa2(:,6)+gmet(2,2)*gxa2(:,2)+gmet(2,3)*gxa2(:,4)
297 : r2(:,2,3)=gmet(2,1)*gxa2(:,5)+gmet(2,2)*gxa2(:,4)+gmet(2,3)*gxa2(:,3)
298 857688 : r2(:,3,1)=gmet(3,1)*gxa2(:,1)+gmet(3,2)*gxa2(:,6)+gmet(3,3)*gxa2(:,5)
299 857688 : r2(:,3,2)=gmet(3,1)*gxa2(:,6)+gmet(3,2)*gxa2(:,2)+gmet(3,3)*gxa2(:,4)
300 : r2(:,3,3)=gmet(3,1)*gxa2(:,5)+gmet(3,2)*gxa2(:,4)+gmet(3,3)*gxa2(:,3)
301 :
302 : !Then compute rank2c(a,b) = 3 conjg(gxa1(a,i)) r2(i,b)
303 : !stored as 11 22 33 32 31 21
304 : !rank2c(re,1) = 3.d0*(gxa1(re,1)*r2(re,1,1)+gxa1(im,1)*r2(im,1,1)&
305 : !& +gxa1(re,6)*r2(re,2,1)+gxa1(im,6)*r2(im,2,1)&
306 : !& +gxa1(re,5)*r2(re,3,1)+gxa1(im,5)*r2(im,3,1))
307 : !rank2c(re,2) = 3.d0*(gxa1(re,6)*r2(re,1,2)+gxa1(im,6)*r2(im,1,2)&
308 : !& +gxa1(re,2)*r2(re,2,2)+gxa1(im,2)*r2(im,2,2)&
309 : !& +gxa1(re,4)*r2(re,3,2)+gxa1(im,4)*r2(im,3,2))
310 : !rank2c(re,3) = 3.d0*(gxa1(re,5)*r2(re,1,3)+gxa1(im,5)*r2(im,1,3)&
311 : !& +gxa1(re,4)*r2(re,2,3)+gxa1(im,4)*r2(im,2,3)&
312 : !& +gxa1(re,3)*r2(re,3,3)+gxa1(im,3)*r2(im,3,3))
313 : rank2c(re,4) = 3.d0*(gxa1(re,5)*r2(re,1,2)+gxa1(im,5)*r2(im,1,2)&
314 : & +gxa1(re,4)*r2(re,2,2)+gxa1(im,4)*r2(im,2,2)&
315 285896 : & +gxa1(re,3)*r2(re,3,2)+gxa1(im,3)*r2(im,3,2))
316 : rank2c(re,5) = 3.d0*(gxa1(re,5)*r2(re,1,1)+gxa1(im,5)*r2(im,1,1)&
317 : & +gxa1(re,4)*r2(re,2,1)+gxa1(im,4)*r2(im,2,1)&
318 285896 : & +gxa1(re,3)*r2(re,3,1)+gxa1(im,3)*r2(im,3,1))
319 : rank2c(re,6) = 3.d0*(gxa1(re,6)*r2(re,1,1)+gxa1(im,6)*r2(im,1,1)&
320 : & +gxa1(re,2)*r2(re,2,1)+gxa1(im,2)*r2(im,2,1)&
321 285896 : & +gxa1(re,4)*r2(re,3,1)+gxa1(im,4)*r2(im,3,1))
322 : !rank2c(im,1) = 3.d0*(gxa1(re,1)*r2(im,1,1)-gxa1(im,1)*r2(re,1,1)&
323 : !& +gxa1(re,6)*r2(im,2,1)-gxa1(im,6)*r2(re,2,1)&
324 : !& +gxa1(re,5)*r2(im,3,1)-gxa1(im,5)*r2(re,3,1))
325 : !rank2c(im,2) = 3.d0*(gxa1(re,6)*r2(im,1,2)-gxa1(im,6)*r2(re,1,2)&
326 : !& +gxa1(re,2)*r2(im,2,2)-gxa1(im,2)*r2(re,2,2)&
327 : !& +gxa1(re,4)*r2(im,3,2)-gxa1(im,4)*r2(re,3,2))
328 : !rank2c(im,3) = 3.d0*(gxa1(re,5)*r2(im,1,3)-gxa1(im,5)*r2(re,1,3)&
329 : !& +gxa1(re,4)*r2(im,2,3)-gxa1(im,4)*r2(re,2,3)&
330 : !& +gxa1(re,3)*r2(im,3,3)-gxa1(im,3)*r2(re,3,3))
331 : rank2c(im,4) = 3.d0*(gxa1(re,5)*r2(im,1,2)-gxa1(im,5)*r2(re,1,2)&
332 : & +gxa1(re,4)*r2(im,2,2)-gxa1(im,4)*r2(re,2,2)&
333 285896 : & +gxa1(re,3)*r2(im,3,2)-gxa1(im,3)*r2(re,3,2))
334 : rank2c(im,5) = 3.d0*(gxa1(re,5)*r2(im,1,1)-gxa1(im,5)*r2(re,1,1)&
335 : & +gxa1(re,4)*r2(im,2,1)-gxa1(im,4)*r2(re,2,1)&
336 285896 : & +gxa1(re,3)*r2(im,3,1)-gxa1(im,3)*r2(re,3,1))
337 : rank2c(im,6) = 3.d0*(gxa1(re,6)*r2(im,1,1)-gxa1(im,6)*r2(re,1,1)&
338 : & +gxa1(re,2)*r2(im,2,1)-gxa1(im,2)*r2(re,2,1)&
339 285896 : & +gxa1(re,4)*r2(im,3,1)-gxa1(im,4)*r2(re,3,1))
340 :
341 285896 : end subroutine cont22cso
342 : !!***
343 :
344 :
345 : !!****f* m_contract/cont22so
346 : !! NAME
347 : !! cont22so
348 : !!
349 : !! FUNCTION
350 : !! Contract symmetric rank 2 tensor gxa1 with symmetric rank 2 tensor
351 : !! gxa2 using antisymmetric tensor amet to produce rank 2 real tensor.
352 : !!
353 : !! INPUTS
354 : !! gxa1(2,6)=rank 2 complex symmetric tensor
355 : !! gxa2(2,6)=rank 2 complex symmetric tensor
356 : !! amet(2,3,3)=antisymmetric complex tensor used for spin-orbit
357 : !!
358 : !! OUTPUT
359 : !! rank2(6)=rank 2 real tensor (pseudo-symmetric storage)
360 : !!
361 : !! NOTES
362 : !! This contraction is used for spin-orbit correction in non-local
363 : !! contribution to stresses.
364 : !!
365 : !! Symmetric gxa1, gxa2 are stored as 11 22 33 32 31 21;
366 : !! amet(3,3) is antisymmetric but stored fully (9 elements);
367 : !! Output rank2 is not symmetric but since
368 : !! $rank2_{gxa1,gxa2}(a,b)=conjg(rank2_{gxa2,gxa1}(b,a))$
369 : !! it is stored as 11 22 33 32 31 21.
370 : !! Want 2*Re[contraction].
371 : !!
372 : !!{{\ \begin{equation}
373 : !! rank2(a,b)=2 Re[conjg(gxa1(i,a)) amet(i,j) gxa2(j,b)]
374 : !!\end{equation} }}
375 : !!
376 : !! Note that, since amet is antisymmetric, amet(i,i)=0
377 : !!
378 : !! SOURCE
379 :
380 285896 : subroutine cont22so(gxa1,gxa2,amet,rank2)
381 :
382 : !Arguments ------------------------------------
383 : !arrays
384 : real(dp),intent(in) :: amet(2,3,3),gxa1(2,6),gxa2(2,6)
385 : real(dp),intent(out) :: rank2(6)
386 :
387 : !Local variables-------------------------------
388 : !scalars
389 : integer,parameter :: im=2,re=1
390 :
391 : ! *************************************************************************
392 :
393 : !Simply write out index summations
394 : !a=1, b=1 in rank2(a,b) --> maps to index 1
395 : rank2(1)=2.0d0*(&
396 : & amet(1,3,2)*(gxa1(re,5)*gxa2(re,6)+gxa1(im,5)*gxa2(im,6))-&
397 : & amet(2,3,2)*(gxa1(re,5)*gxa2(im,6)-gxa1(im,5)*gxa2(re,6))+&
398 : & amet(1,3,1)*(gxa1(re,5)*gxa2(re,1)+gxa1(im,5)*gxa2(im,1))-&
399 : & amet(2,3,1)*(gxa1(re,5)*gxa2(im,1)-gxa1(im,5)*gxa2(re,1))+&
400 : & amet(1,2,1)*(gxa1(re,6)*gxa2(re,1)+gxa1(im,6)*gxa2(im,1))-&
401 : & amet(2,2,1)*(gxa1(re,6)*gxa2(im,1)-gxa1(im,6)*gxa2(re,1))+&
402 : & amet(1,2,3)*(gxa1(re,6)*gxa2(re,5)+gxa1(im,6)*gxa2(im,5))-&
403 : & amet(2,2,3)*(gxa1(re,6)*gxa2(im,5)-gxa1(im,6)*gxa2(re,5))+&
404 : & amet(1,1,3)*(gxa1(re,1)*gxa2(re,5)+gxa1(im,1)*gxa2(im,5))-&
405 : & amet(2,1,3)*(gxa1(re,1)*gxa2(im,5)-gxa1(im,1)*gxa2(re,5))+&
406 : & amet(1,1,2)*(gxa1(re,1)*gxa2(re,6)+gxa1(im,1)*gxa2(im,6))-&
407 285896 : & amet(2,1,2)*(gxa1(re,1)*gxa2(im,6)-gxa1(im,1)*gxa2(re,6)))
408 :
409 : !a=2, b=2 in rank2(a,b) --> maps to index 2
410 : rank2(2)=2.0d0*(&
411 : & amet(1,3,2)*(gxa1(re,4)*gxa2(re,2)+gxa1(im,4)*gxa2(im,2))-&
412 : & amet(2,3,2)*(gxa1(re,4)*gxa2(im,2)-gxa1(im,4)*gxa2(re,2))+&
413 : & amet(1,3,1)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-&
414 : & amet(2,3,1)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+&
415 : & amet(1,2,1)*(gxa1(re,2)*gxa2(re,6)+gxa1(im,2)*gxa2(im,6))-&
416 : & amet(2,2,1)*(gxa1(re,2)*gxa2(im,6)-gxa1(im,2)*gxa2(re,6))+&
417 : & amet(1,2,3)*(gxa1(re,2)*gxa2(re,4)+gxa1(im,2)*gxa2(im,4))-&
418 : & amet(2,2,3)*(gxa1(re,2)*gxa2(im,4)-gxa1(im,2)*gxa2(re,4))+&
419 : & amet(1,1,3)*(gxa1(re,6)*gxa2(re,4)+gxa1(im,6)*gxa2(im,4))-&
420 : & amet(2,1,3)*(gxa1(re,6)*gxa2(im,4)-gxa1(im,6)*gxa2(re,4))+&
421 : & amet(1,1,2)*(gxa1(re,6)*gxa2(re,2)+gxa1(im,6)*gxa2(im,2))-&
422 285896 : & amet(2,1,2)*(gxa1(re,6)*gxa2(im,2)-gxa1(im,6)*gxa2(re,2)))
423 :
424 : !a=3, b=3 in rank2(a,b) --> maps to index 3
425 : rank2(3)=2.0d0*(&
426 : & amet(1,3,2)*(gxa1(re,3)*gxa2(re,4)+gxa1(im,3)*gxa2(im,4))-&
427 : & amet(2,3,2)*(gxa1(re,3)*gxa2(im,4)-gxa1(im,3)*gxa2(re,4))+&
428 : & amet(1,3,1)*(gxa1(re,3)*gxa2(re,5)+gxa1(im,3)*gxa2(im,5))-&
429 : & amet(2,3,1)*(gxa1(re,3)*gxa2(im,5)-gxa1(im,3)*gxa2(re,5))+&
430 : & amet(1,2,1)*(gxa1(re,4)*gxa2(re,5)+gxa1(im,4)*gxa2(im,5))-&
431 : & amet(2,2,1)*(gxa1(re,4)*gxa2(im,5)-gxa1(im,4)*gxa2(re,5))+&
432 : & amet(1,2,3)*(gxa1(re,4)*gxa2(re,3)+gxa1(im,4)*gxa2(im,3))-&
433 : & amet(2,2,3)*(gxa1(re,4)*gxa2(im,3)-gxa1(im,4)*gxa2(re,3))+&
434 : & amet(1,1,3)*(gxa1(re,5)*gxa2(re,3)+gxa1(im,5)*gxa2(im,3))-&
435 : & amet(2,1,3)*(gxa1(re,5)*gxa2(im,3)-gxa1(im,5)*gxa2(re,3))+&
436 : & amet(1,1,2)*(gxa1(re,5)*gxa2(re,4)+gxa1(im,5)*gxa2(im,4))-&
437 285896 : & amet(2,1,2)*(gxa1(re,5)*gxa2(im,4)-gxa1(im,5)*gxa2(re,4)))
438 :
439 : !a=3, b=2 in rank2(a,b) --> maps to index 4
440 : rank2(4)=2.0d0*(&
441 : & amet(1,3,2)*(gxa1(re,3)*gxa2(re,2)+gxa1(im,3)*gxa2(im,2))-&
442 : & amet(2,3,2)*(gxa1(re,3)*gxa2(im,2)-gxa1(im,3)*gxa2(re,2))+&
443 : & amet(1,3,1)*(gxa1(re,3)*gxa2(re,6)+gxa1(im,3)*gxa2(im,6))-&
444 : & amet(2,3,1)*(gxa1(re,3)*gxa2(im,6)-gxa1(im,3)*gxa2(re,6))+&
445 : & amet(1,2,1)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-&
446 : & amet(2,2,1)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+&
447 : & amet(1,2,3)*(gxa1(re,4)*gxa2(re,4)+gxa1(im,4)*gxa2(im,4))-&
448 : & amet(2,2,3)*(gxa1(re,4)*gxa2(im,4)-gxa1(im,4)*gxa2(re,4))+&
449 : & amet(1,1,3)*(gxa1(re,5)*gxa2(re,4)+gxa1(im,5)*gxa2(im,4))-&
450 : & amet(2,1,3)*(gxa1(re,5)*gxa2(im,4)-gxa1(im,5)*gxa2(re,4))+&
451 : & amet(1,1,2)*(gxa1(re,5)*gxa2(re,2)+gxa1(im,5)*gxa2(im,2))-&
452 285896 : & amet(2,1,2)*(gxa1(re,5)*gxa2(im,2)-gxa1(im,5)*gxa2(re,2)))
453 :
454 : !a=3, b=1 in rank2(a,b) --> maps to index 5
455 : rank2(5)=2.0d0*(&
456 : & amet(1,3,2)*(gxa1(re,3)*gxa2(re,6)+gxa1(im,3)*gxa2(im,6))-&
457 : & amet(2,3,2)*(gxa1(re,3)*gxa2(im,6)-gxa1(im,3)*gxa2(re,6))+&
458 : & amet(1,3,1)*(gxa1(re,3)*gxa2(re,1)+gxa1(im,3)*gxa2(im,1))-&
459 : & amet(2,3,1)*(gxa1(re,3)*gxa2(im,1)-gxa1(im,3)*gxa2(re,1))+&
460 : & amet(1,2,1)*(gxa1(re,4)*gxa2(re,1)+gxa1(im,4)*gxa2(im,1))-&
461 : & amet(2,2,1)*(gxa1(re,4)*gxa2(im,1)-gxa1(im,4)*gxa2(re,1))+&
462 : & amet(1,2,3)*(gxa1(re,4)*gxa2(re,5)+gxa1(im,4)*gxa2(im,5))-&
463 : & amet(2,2,3)*(gxa1(re,4)*gxa2(im,5)-gxa1(im,4)*gxa2(re,5))+&
464 : & amet(1,1,3)*(gxa1(re,5)*gxa2(re,5)+gxa1(im,5)*gxa2(im,5))-&
465 : & amet(2,1,3)*(gxa1(re,5)*gxa2(im,5)-gxa1(im,5)*gxa2(re,5))+&
466 : & amet(1,1,2)*(gxa1(re,5)*gxa2(re,6)+gxa1(im,5)*gxa2(im,6))-&
467 285896 : & amet(2,1,2)*(gxa1(re,5)*gxa2(im,6)-gxa1(im,5)*gxa2(re,6)))
468 :
469 : !a=2, b=1 in rank2(a,b) --> maps to index 6
470 : rank2(6)=2.0d0*(&
471 : & amet(1,3,2)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-&
472 : & amet(2,3,2)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+&
473 : & amet(1,3,1)*(gxa1(re,4)*gxa2(re,1)+gxa1(im,4)*gxa2(im,1))-&
474 : & amet(2,3,1)*(gxa1(re,4)*gxa2(im,1)-gxa1(im,4)*gxa2(re,1))+&
475 : & amet(1,2,1)*(gxa1(re,2)*gxa2(re,1)+gxa1(im,2)*gxa2(im,1))-&
476 : & amet(2,2,1)*(gxa1(re,2)*gxa2(im,1)-gxa1(im,2)*gxa2(re,1))+&
477 : & amet(1,2,3)*(gxa1(re,2)*gxa2(re,5)+gxa1(im,2)*gxa2(im,5))-&
478 : & amet(2,2,3)*(gxa1(re,2)*gxa2(im,5)-gxa1(im,2)*gxa2(re,5))+&
479 : & amet(1,1,3)*(gxa1(re,6)*gxa2(re,5)+gxa1(im,6)*gxa2(im,5))-&
480 : & amet(2,1,3)*(gxa1(re,6)*gxa2(im,5)-gxa1(im,6)*gxa2(re,5))+&
481 : & amet(1,1,2)*(gxa1(re,6)*gxa2(re,6)+gxa1(im,6)*gxa2(im,6))-&
482 285896 : & amet(2,1,2)*(gxa1(re,6)*gxa2(im,6)-gxa1(im,6)*gxa2(re,6)))
483 :
484 285896 : end subroutine cont22so
485 : !!***
486 :
487 :
488 : !!****f* m_contract/cont24
489 : !! NAME
490 : !! cont24
491 : !!
492 : !! FUNCTION
493 : !! Contract symmetric rank2 tensor gxa with rank4 symmetric tensor to
494 : !! produce symmetric rank2 tensor.
495 : !!
496 : !! INPUTS
497 : !! gxa(2,6)=rank 2 symmetric complex tensor in order 11 22 33 32 31 21
498 : !! rank4(2,15)=rank 4 complex tensor (symmetric storage)
499 : !!
500 : !! OUTPUT
501 : !! rank2(6)=rank 2 real tensor (symmetric storage) 11 22 33 32 31 21.
502 : !!
503 : !! NOTES
504 : !! Tensors are in "symmetric" storage mode.
505 : !! for gxa and rank2 this is 11, 22, 33, 32, 31, 21;
506 : !! for the rank 4 tensor rank4 this is
507 : !! 1111 2211 3311 3211 3111 2111 2221 3321 3221 3331 2222 3322 3222 3332 3333.
508 : !! gxa and rank4 are complex; rank2 is real.
509 : !! Want $2 Re[contraction]$.
510 : !! $rank2(a,b)=2 Re[gxa(i,j)^"*" rank4(a,b,i,j)]$.
511 : !!
512 : !! Note that the input gxa is typically the result of
513 : !! $gxa(i,j)=[{3 \over 2} gmet(i,l) gmet(j,m) - {1 \over 2} gmet(i,j) gmet(l,m)] gxa_old(l,m)$
514 : !! where the subroutine "metcon" already includes weights in the
515 : !! definition of gxa for off-diagonal elements (weight of 2 for
516 : !! symmetry).
517 : !! Components 4, 5, and 6 of gxa have already been multiplied by 2
518 : !! so the expressions below do not carry the 2.
519 : !!
520 : !! SOURCE
521 :
522 634069 : subroutine cont24(gxa,rank4,rank2)
523 :
524 : !Arguments ------------------------------------
525 : !arrays
526 : real(dp),intent(in) :: gxa(2,6),rank4(2,15)
527 : real(dp),intent(out) :: rank2(6)
528 :
529 : !Local variables-------------------------------
530 : !scalars
531 : integer,parameter :: im=2,re=1
532 :
533 : ! *************************************************************************
534 :
535 : !Simply write out index summations
536 :
537 : !a=1, b=1 in rank2(a,b) --> maps to index 1
538 : rank2(1)=2.0d0*(&
539 : & gxa(re,1)*rank4(re, 1)+gxa(im,1)*rank4(im, 1)+&
540 : & gxa(re,2)*rank4(re, 2)+gxa(im,2)*rank4(im, 2)+&
541 : & gxa(re,3)*rank4(re, 3)+gxa(im,3)*rank4(im, 3)+&
542 : & gxa(re,4)*rank4(re, 4)+gxa(im,4)*rank4(im, 4)+&
543 : & gxa(re,5)*rank4(re, 5)+gxa(im,5)*rank4(im, 5)+&
544 634069 : & gxa(re,6)*rank4(re, 6)+gxa(im,6)*rank4(im, 6))
545 :
546 : !a=2, b=2 in rank2(a,b) --> maps to index 2
547 : rank2(2)=2.0d0*(&
548 : & gxa(re,1)*rank4(re, 2)+gxa(im,1)*rank4(im, 2)+&
549 : & gxa(re,2)*rank4(re,11)+gxa(im,2)*rank4(im,11)+&
550 : & gxa(re,3)*rank4(re,12)+gxa(im,3)*rank4(im,12)+&
551 : & gxa(re,4)*rank4(re,13)+gxa(im,4)*rank4(im,13)+&
552 : & gxa(re,5)*rank4(re, 9)+gxa(im,5)*rank4(im, 9)+&
553 634069 : & gxa(re,6)*rank4(re, 7)+gxa(im,6)*rank4(im, 7))
554 :
555 : !a=3, b=3 in rank2(a,b) --> maps to index 3
556 : rank2(3)=2.0d0*(&
557 : & gxa(re,1)*rank4(re, 3)+gxa(im,1)*rank4(im, 3)+&
558 : & gxa(re,2)*rank4(re,12)+gxa(im,2)*rank4(im,12)+&
559 : & gxa(re,3)*rank4(re,15)+gxa(im,3)*rank4(im,15)+&
560 : & gxa(re,4)*rank4(re,14)+gxa(im,4)*rank4(im,14)+&
561 : & gxa(re,5)*rank4(re,10)+gxa(im,5)*rank4(im,10)+&
562 634069 : & gxa(re,6)*rank4(re, 8)+gxa(im,6)*rank4(im, 8))
563 :
564 : !a=3, b=2 in rank2(a,b) --> maps to index 4
565 : rank2(4)=2.0d0*(&
566 : & gxa(re,1)*rank4(re, 4)+gxa(im,1)*rank4(im, 4)+&
567 : & gxa(re,2)*rank4(re,13)+gxa(im,2)*rank4(im,13)+&
568 : & gxa(re,3)*rank4(re,14)+gxa(im,3)*rank4(im,14)+&
569 : & gxa(re,4)*rank4(re,12)+gxa(im,4)*rank4(im,12)+&
570 : & gxa(re,5)*rank4(re, 8)+gxa(im,5)*rank4(im, 8)+&
571 634069 : & gxa(re,6)*rank4(re, 9)+gxa(im,6)*rank4(im, 9))
572 :
573 : !a=3, b=1 in rank2(a,b) --> maps to index 5
574 : rank2(5)=2.0d0*(&
575 : & gxa(re,1)*rank4(re, 5)+gxa(im,1)*rank4(im, 5)+&
576 : & gxa(re,2)*rank4(re, 9)+gxa(im,2)*rank4(im, 9)+&
577 : & gxa(re,3)*rank4(re,10)+gxa(im,3)*rank4(im,10)+&
578 : & gxa(re,4)*rank4(re, 8)+gxa(im,4)*rank4(im, 8)+&
579 : & gxa(re,5)*rank4(re, 3)+gxa(im,5)*rank4(im, 3)+&
580 634069 : & gxa(re,6)*rank4(re, 4)+gxa(im,6)*rank4(im, 4))
581 :
582 : !a=2, b=1 in rank2(a,b) --> maps to index 6
583 : rank2(6)=2.0d0*(&
584 : & gxa(re,1)*rank4(re, 6)+gxa(im,1)*rank4(im, 6)+&
585 : & gxa(re,2)*rank4(re, 7)+gxa(im,2)*rank4(im, 7)+&
586 : & gxa(re,3)*rank4(re, 8)+gxa(im,3)*rank4(im, 8)+&
587 : & gxa(re,4)*rank4(re, 9)+gxa(im,4)*rank4(im, 9)+&
588 : & gxa(re,5)*rank4(re, 4)+gxa(im,5)*rank4(im, 4)+&
589 634069 : & gxa(re,6)*rank4(re, 2)+gxa(im,6)*rank4(im, 2))
590 :
591 634069 : end subroutine cont24
592 : !!***
593 :
594 : !!****f* m_contract/cont3
595 : !! NAME
596 : !! cont3
597 : !!
598 : !! FUNCTION
599 : !! Compute several specialized contractions needed for the
600 : !! l=3 part of the stress tensor.
601 : !!
602 : !! INPUTS
603 : !! gxa(2,10)=complex symmetric rank 3 tensor
604 : !! gmet(3,3)=usual metric tensor, a symmetric matrix stored in
605 : !! full storage mode (bohr^-2)
606 : !!
607 : !! OUTPUT
608 : !! rank2(6)=2*Re[contraction] given by
609 : !! 2*Re[(15/2)*r3(a,i,j)*r3(b,j,i)-3*r1(i)*r3(a,b,i)-(3/2)*r1(a)*r1(b)]
610 : !! where r3(a,i,j)=gmet(j,k) gxa(a,i,k) and r1(a)=gmet(i,j) gxa(i,j,a).
611 : !! rank2 is stored in the compressed form 11 22 33 32 31 21.
612 : !!
613 : !! NOTES
614 : !! Input gxa is a completely symmetric rank 3 tensor (complex)
615 : !! in compressed storage: 111 221 331 321 311 211 222 332 322 333.
616 : !! The output tensor is completely symmetric rank 2, real, and is given by
617 : !! $2 Re[{15 \over 2} r3(a,i,j) r3(b,j,i) - 3 r1(i) r3(a,b,i) - {3 \over 2} r1(a) r1(b)]$
618 : !! where $r3(a,i,j)=gmet(j,k) gxa(a,i,k)$ and $r1(a)=gmet(i,j) gxa(i,j,a)$.
619 : !! rank2 is stored in the compressed form 11 22 33 32 31 21.
620 : !!
621 : !! SOURCE
622 :
623 37366 : subroutine cont3(gxa,gmet,rank2)
624 :
625 : !Arguments ------------------------------------
626 : !arrays
627 : real(dp),intent(in) :: gmet(3,3),gxa(2,10)
628 : real(dp),intent(out) :: rank2(6)
629 :
630 : !Local variables-------------------------------
631 : !scalars
632 : integer,parameter :: im=2,re=1
633 : integer :: ii
634 : !arrays
635 : real(dp) :: r1(2,3),r3(2,18),s13(6),s33(6)
636 :
637 : ! *************************************************************************
638 :
639 : !Compute r1(a) = gmet(i,j) gxa(i,j,a)
640 :
641 : !Write out components for 3 distinct terms, Re and Im
642 112098 : do ii=1,2
643 : r1(ii,1)=gmet(1,1)*gxa(ii,1)+gmet(2,2)*gxa(ii,2)+&
644 : & gmet(3,3)*gxa(ii,3)+2.d0*(&
645 : & gmet(3,2)*gxa(ii,4)+gmet(3,1)*gxa(ii,5)+&
646 74732 : & gmet(2,1)*gxa(ii,6))
647 : r1(ii,2)=gmet(1,1)*gxa(ii,6)+gmet(2,2)*gxa(ii,7)+&
648 : & gmet(3,3)*gxa(ii,8)+2.d0*(&
649 : & gmet(3,2)*gxa(ii,9)+gmet(3,1)*gxa(ii,4)+&
650 74732 : & gmet(2,1)*gxa(ii,2))
651 : r1(ii,3)=gmet(1,1)*gxa(ii,5)+gmet(2,2)*gxa(ii,9)+&
652 : & gmet(3,3)*gxa(ii,10)+2.d0*(&
653 : & gmet(3,2)*gxa(ii,8)+gmet(3,1)*gxa(ii,3)+&
654 112098 : & gmet(2,1)*gxa(ii,4))
655 : end do
656 :
657 : !Compute r3(a,b,k)=gmet(k,n) gxa(a,b,n)
658 :
659 : !Write out components for 18 distinct terms, Re and Im
660 : !(symmetric in first two indices, not in all permutations)
661 : !store as 111 221 331 321 311 211
662 : !112 222 332 322 312 212
663 : !113 223 333 323 313 213
664 112098 : do ii=1,2
665 : r3(ii, 1)=gmet(1,1)*gxa(ii,1)+gmet(2,1)*gxa(ii,6)+&
666 74732 : & gmet(3,1)*gxa(ii,5)
667 : r3(ii, 2)=gmet(1,1)*gxa(ii,2)+gmet(2,1)*gxa(ii,7)+&
668 74732 : & gmet(3,1)*gxa(ii,9)
669 : r3(ii, 3)=gmet(1,1)*gxa(ii,3)+gmet(2,1)*gxa(ii,8)+&
670 74732 : & gmet(3,1)*gxa(ii,10)
671 : r3(ii, 4)=gmet(1,1)*gxa(ii,4)+gmet(2,1)*gxa(ii,9)+&
672 74732 : & gmet(3,1)*gxa(ii,8)
673 : r3(ii, 5)=gmet(1,1)*gxa(ii,5)+gmet(2,1)*gxa(ii,4)+&
674 74732 : & gmet(3,1)*gxa(ii,3)
675 : r3(ii, 6)=gmet(1,1)*gxa(ii,6)+gmet(2,1)*gxa(ii,2)+&
676 74732 : & gmet(3,1)*gxa(ii,4)
677 : r3(ii, 7)=gmet(2,1)*gxa(ii,1)+gmet(2,2)*gxa(ii,6)+&
678 74732 : & gmet(3,2)*gxa(ii,5)
679 : r3(ii, 8)=gmet(2,1)*gxa(ii,2)+gmet(2,2)*gxa(ii,7)+&
680 74732 : & gmet(3,2)*gxa(ii,9)
681 : r3(ii, 9)=gmet(2,1)*gxa(ii,3)+gmet(2,2)*gxa(ii,8)+&
682 74732 : & gmet(3,2)*gxa(ii,10)
683 : r3(ii,10)=gmet(2,1)*gxa(ii,4)+gmet(2,2)*gxa(ii,9)+&
684 74732 : & gmet(3,2)*gxa(ii,8)
685 : r3(ii,11)=gmet(2,1)*gxa(ii,5)+gmet(2,2)*gxa(ii,4)+&
686 74732 : & gmet(3,2)*gxa(ii,3)
687 : r3(ii,12)=gmet(2,1)*gxa(ii,6)+gmet(2,2)*gxa(ii,2)+&
688 74732 : & gmet(3,2)*gxa(ii,4)
689 : r3(ii,13)=gmet(3,1)*gxa(ii,1)+gmet(3,2)*gxa(ii,6)+&
690 74732 : & gmet(3,3)*gxa(ii,5)
691 : r3(ii,14)=gmet(3,1)*gxa(ii,2)+gmet(3,2)*gxa(ii,7)+&
692 74732 : & gmet(3,3)*gxa(ii,9)
693 : r3(ii,15)=gmet(3,1)*gxa(ii,3)+gmet(3,2)*gxa(ii,8)+&
694 74732 : & gmet(3,3)*gxa(ii,10)
695 : r3(ii,16)=gmet(3,1)*gxa(ii,4)+gmet(3,2)*gxa(ii,9)+&
696 74732 : & gmet(3,3)*gxa(ii,8)
697 : r3(ii,17)=gmet(3,1)*gxa(ii,5)+gmet(3,2)*gxa(ii,4)+&
698 74732 : & gmet(3,3)*gxa(ii,3)
699 : r3(ii,18)=gmet(3,1)*gxa(ii,6)+gmet(3,2)*gxa(ii,2)+&
700 112098 : & gmet(3,3)*gxa(ii,4)
701 :
702 : end do
703 :
704 : !Now need
705 : !2*Re[(15/2)*r3(a,i,j)*r3(b,j,i)-3*r1(i)*r3(a,b,i)-(3/2)*r1(a)*r1(b)].
706 :
707 : !Write out s33(a,b)=2*Re[r3(a,i,j)*r3(b,j,i)]
708 :
709 : s33(1)=2.d0*(r3(re, 1)*r3(re, 1)+r3(im, 1)*r3(im, 1)+&
710 : & r3(re,12)*r3(re,12)+r3(im,12)*r3(im,12)+&
711 : & r3(re,17)*r3(re,17)+r3(im,17)*r3(im,17)+&
712 : & r3(re,11)*r3(re,18)+r3(im,11)*r3(im,18)+&
713 : & r3(re,18)*r3(re,11)+r3(im,18)*r3(im,11)+&
714 : & r3(re, 5)*r3(re,13)+r3(im, 5)*r3(im,13)+&
715 : & r3(re,13)*r3(re, 5)+r3(im,13)*r3(im, 5)+&
716 : & r3(re, 6)*r3(re, 7)+r3(im, 6)*r3(im, 7)+&
717 37366 : & r3(re, 7)*r3(re, 6)+r3(im, 7)*r3(im, 6))
718 :
719 : s33(2)=2.d0*(r3(re, 6)*r3(re, 6)+r3(im, 6)*r3(im, 6)+&
720 : & r3(re, 8)*r3(re, 8)+r3(im, 8)*r3(im, 8)+&
721 : & r3(re,16)*r3(re,16)+r3(im,16)*r3(im,16)+&
722 : & r3(re,10)*r3(re,14)+r3(im,10)*r3(im,14)+&
723 : & r3(re,14)*r3(re,10)+r3(im,14)*r3(im,10)+&
724 : & r3(re, 4)*r3(re,18)+r3(im, 4)*r3(im,18)+&
725 : & r3(re,18)*r3(re, 4)+r3(im,18)*r3(im, 4)+&
726 : & r3(re, 2)*r3(re,12)+r3(im, 2)*r3(im,12)+&
727 37366 : & r3(re,12)*r3(re, 2)+r3(im,12)*r3(im, 2))
728 :
729 : s33(3)=2.d0*(r3(re, 5)*r3(re, 5)+r3(im, 5)*r3(im, 5)+&
730 : & r3(re,10)*r3(re,10)+r3(im,10)*r3(im,10)+&
731 : & r3(re,15)*r3(re,15)+r3(im,15)*r3(im,15)+&
732 : & r3(re, 9)*r3(re,16)+r3(im, 9)*r3(im,16)+&
733 : & r3(re,16)*r3(re, 9)+r3(im,16)*r3(im, 9)+&
734 : & r3(re, 3)*r3(re,17)+r3(im, 3)*r3(im,17)+&
735 : & r3(re,17)*r3(re, 3)+r3(im,17)*r3(im, 3)+&
736 : & r3(re, 4)*r3(re,11)+r3(im, 4)*r3(im,11)+&
737 37366 : & r3(re,11)*r3(re, 4)+r3(im,11)*r3(im, 4))
738 :
739 : s33(4)=2.d0*(r3(re, 5)*r3(re, 6)+r3(im, 5)*r3(im, 6)+&
740 : & r3(re,10)*r3(re, 8)+r3(im,10)*r3(im, 8)+&
741 : & r3(re,15)*r3(re,16)+r3(im,15)*r3(im,16)+&
742 : & r3(re, 9)*r3(re,14)+r3(im, 9)*r3(im,14)+&
743 : & r3(re,16)*r3(re,10)+r3(im,16)*r3(im,10)+&
744 : & r3(re, 3)*r3(re,18)+r3(im, 3)*r3(im,18)+&
745 : & r3(re,17)*r3(re, 4)+r3(im,17)*r3(im, 4)+&
746 : & r3(re, 4)*r3(re,12)+r3(im, 4)*r3(im,12)+&
747 37366 : & r3(re,11)*r3(re, 2)+r3(im,11)*r3(im, 2))
748 :
749 : s33(5)=2.d0*(r3(re, 5)*r3(re, 1)+r3(im, 5)*r3(im, 1)+&
750 : & r3(re,10)*r3(re,12)+r3(im,10)*r3(im,12)+&
751 : & r3(re,15)*r3(re,17)+r3(im,15)*r3(im,17)+&
752 : & r3(re, 9)*r3(re,18)+r3(im, 9)*r3(im,18)+&
753 : & r3(re,16)*r3(re,11)+r3(im,16)*r3(im,11)+&
754 : & r3(re, 3)*r3(re,13)+r3(im, 3)*r3(im,13)+&
755 : & r3(re,17)*r3(re, 5)+r3(im,17)*r3(im, 5)+&
756 : & r3(re, 4)*r3(re, 7)+r3(im, 4)*r3(im, 7)+&
757 37366 : & r3(re,11)*r3(re, 6)+r3(im,11)*r3(im, 6))
758 :
759 : s33(6)=2.d0*(r3(re, 6)*r3(re, 1)+r3(im, 6)*r3(im, 1)+&
760 : & r3(re, 8)*r3(re,12)+r3(im, 8)*r3(im,12)+&
761 : & r3(re,16)*r3(re,17)+r3(im,16)*r3(im,17)+&
762 : & r3(re,10)*r3(re,18)+r3(im,10)*r3(im,18)+&
763 : & r3(re,14)*r3(re,11)+r3(im,14)*r3(im,11)+&
764 : & r3(re, 4)*r3(re,13)+r3(im, 4)*r3(im,13)+&
765 : & r3(re,18)*r3(re, 5)+r3(im,18)*r3(im, 5)+&
766 : & r3(re, 2)*r3(re, 7)+r3(im, 2)*r3(im, 7)+&
767 37366 : & r3(re,12)*r3(re, 6)+r3(im,12)*r3(im, 6))
768 :
769 :
770 : !Write out s13(a,b)=2*Re[r1(i)*r3(a,b,i)]
771 :
772 : s13(1)=2.d0*(r1(re,1)*r3(re, 1)+r1(im,1)*r3(im, 1)+&
773 : & r1(re,2)*r3(re, 7)+r1(im,2)*r3(im, 7)+&
774 37366 : & r1(re,3)*r3(re,13)+r1(im,3)*r3(im,13))
775 : s13(2)=2.d0*(r1(re,1)*r3(re, 2)+r1(im,1)*r3(im, 2)+&
776 : & r1(re,2)*r3(re, 8)+r1(im,2)*r3(im, 8)+&
777 37366 : & r1(re,3)*r3(re,14)+r1(im,3)*r3(im,14))
778 : s13(3)=2.d0*(r1(re,1)*r3(re, 3)+r1(im,1)*r3(im, 3)+&
779 : & r1(re,2)*r3(re, 9)+r1(im,2)*r3(im, 9)+&
780 37366 : & r1(re,3)*r3(re,15)+r1(im,3)*r3(im,15))
781 : s13(4)=2.d0*(r1(re,1)*r3(re, 4)+r1(im,1)*r3(im, 4)+&
782 : & r1(re,2)*r3(re,10)+r1(im,2)*r3(im,10)+&
783 37366 : & r1(re,3)*r3(re,16)+r1(im,3)*r3(im,16))
784 : s13(5)=2.d0*(r1(re,1)*r3(re, 5)+r1(im,1)*r3(im, 5)+&
785 : & r1(re,2)*r3(re,11)+r1(im,2)*r3(im,11)+&
786 37366 : & r1(re,3)*r3(re,17)+r1(im,3)*r3(im,17))
787 : s13(6)=2.d0*(r1(re,1)*r3(re, 6)+r1(im,1)*r3(im, 6)+&
788 : & r1(re,2)*r3(re,12)+r1(im,2)*r3(im,12)+&
789 37366 : & r1(re,3)*r3(re,18)+r1(im,3)*r3(im,18))
790 :
791 : !Finally, write out the six terms as final answer
792 : !rank2(a,b)=(15/2)*s33(a,b)-3*s13(a,b)-(3/2)*2*Re[r1(a)*r1(b)]
793 :
794 : rank2(1)=7.5d0*s33(1)-3.d0*s13(1)&
795 37366 : & -3.d0*(r1(re,1)*r1(re,1)+r1(im,1)*r1(im,1))
796 : rank2(2)=7.5d0*s33(2)-3.d0*s13(2)&
797 37366 : & -3.d0*(r1(re,2)*r1(re,2)+r1(im,2)*r1(im,2))
798 : rank2(3)=7.5d0*s33(3)-3.d0*s13(3)&
799 37366 : & -3.d0*(r1(re,3)*r1(re,3)+r1(im,3)*r1(im,3))
800 : rank2(4)=7.5d0*s33(4)-3.d0*s13(4)&
801 37366 : & -3.d0*(r1(re,3)*r1(re,2)+r1(im,3)*r1(im,2))
802 : rank2(5)=7.5d0*s33(5)-3.d0*s13(5)&
803 37366 : & -3.d0*(r1(re,3)*r1(re,1)+r1(im,3)*r1(im,1))
804 : rank2(6)=7.5d0*s33(6)-3.d0*s13(6)&
805 37366 : & -3.d0*(r1(re,2)*r1(re,1)+r1(im,2)*r1(im,1))
806 :
807 37366 : end subroutine cont3
808 : !!***
809 :
810 : !!****f* m_contract/cont33cso
811 : !! NAME
812 : !! cont33cso
813 : !!
814 : !! FUNCTION
815 : !! Contract symmetric rank 3 tensor gxa1 with symmetric rank 3 tensor
816 : !! gxa2 using metric tensor gmet to produce rank 2 complex tensor.
817 : !!
818 : !! INPUTS
819 : !! gxa1(2,10)=rank 3 complex symmetric tensor
820 : !! gxa2(2,10)=rank 3 complex symmetric tensor
821 : !! gmet(3,3)=usual metric tensor (symmetric, real)
822 : !!
823 : !! OUTPUT
824 : !! rank2c(2,6)=rank 2 complex tensor (pseudo-symmetric storage)
825 : !!
826 : !! NOTES
827 : !! This contraction is used for spin-orbit correction in non-local
828 : !! contribution to stresses.
829 : !!
830 : !! Symmetric gxa1, gxa2 are stored as
831 : !! 111 221 331 321 311 211 222 332 322 333;
832 : !! gmet(3,3) is symmetric but stored fully (9 elements);
833 : !! Output rank2c is not symmetric but since
834 : !! $rank2c_{gxa1,gxa2}(a,b)=conjg(rank2c_{gxa2,gxa1}(b,a))$
835 : !! it is stored as 11 22 33 32 31 21.
836 : !!
837 : !! rank2c(1,1), rank2c(2,2), rank3c(3,3) are not needed;
838 : !! They are not calculated.
839 : !!
840 : !!{{\ \begin{equation}
841 : !! rank2c(a,b)=7.5 conjg(gxa1(a,i,j))*r_3(i,j,b) - 1.5 r_{11}(a)*r_{12}(b)
842 : !!\end{equation} }}
843 : !! where:
844 : !!{{\ \begin{eqnarray}
845 : !! r_3(i,j,b) & = & gxa2(b,l,m) gmet(i,l) gmet(j,m) \nonumber
846 : !! r_{11}(a) & = & conjg(gxa1(a,l,m)) gmet(l,m) \nonumber
847 : !! r_{12}(b) & = & gxa2(b,l,m) gmet(l,m)
848 : !! \end{eqnarray} }}
849 : !!
850 : !! SOURCE
851 :
852 65360 : subroutine cont33cso(gxa1,gxa2,gmet,rank2c)
853 :
854 : !Arguments ------------------------------------
855 : !arrays
856 : real(dp),intent(in) :: gmet(3,3),gxa1(2,10),gxa2(2,10)
857 : real(dp),intent(out) :: rank2c(2,6)
858 :
859 : !Local variables-------------------------------
860 : !scalars
861 : integer,parameter :: im=2,re=1
862 : !arrays
863 : real(dp) :: r11(2,3),r12(2,3),r3(2,6,3),r3a(2,6,3)
864 :
865 : ! *************************************************************************
866 :
867 : !Initialize output tensor
868 65360 : rank2c(:,:)=0.d0
869 :
870 : !Compute r3(i,j,b)=gxa2(b,l,m)*gmet(i,l)*gmet(j,m)
871 : !stored as 11 22 33 32 31 21 for (i,j)
872 : !First compute r3a(b,l,j)=gxa2(b,l,m)*gmet(j,m)
873 : !stored as r3a(re/im,(b,l),j)
874 196080 : r3a(:,1,1)=gxa2(:,1 )*gmet(1,1)+gxa2(:,6 )*gmet(1,2)+gxa2(:,5 )*gmet(1,3)
875 196080 : r3a(:,2,1)=gxa2(:,2 )*gmet(1,1)+gxa2(:,7 )*gmet(1,2)+gxa2(:,9 )*gmet(1,3)
876 : r3a(:,3,1)=gxa2(:,3 )*gmet(1,1)+gxa2(:,8 )*gmet(1,2)+gxa2(:,10)*gmet(1,3)
877 196080 : r3a(:,4,1)=gxa2(:,4 )*gmet(1,1)+gxa2(:,9 )*gmet(1,2)+gxa2(:,8 )*gmet(1,3)
878 196080 : r3a(:,5,1)=gxa2(:,5 )*gmet(1,1)+gxa2(:,4 )*gmet(1,2)+gxa2(:,3 )*gmet(1,3)
879 196080 : r3a(:,6,1)=gxa2(:,6 )*gmet(1,1)+gxa2(:,2 )*gmet(1,2)+gxa2(:,4 )*gmet(1,3)
880 196080 : r3a(:,1,2)=gxa2(:,1 )*gmet(2,1)+gxa2(:,6 )*gmet(2,2)+gxa2(:,5 )*gmet(2,3)
881 196080 : r3a(:,2,2)=gxa2(:,2 )*gmet(2,1)+gxa2(:,7 )*gmet(2,2)+gxa2(:,9 )*gmet(2,3)
882 : r3a(:,3,2)=gxa2(:,3 )*gmet(2,1)+gxa2(:,8 )*gmet(2,2)+gxa2(:,10)*gmet(2,3)
883 196080 : r3a(:,4,2)=gxa2(:,4 )*gmet(2,1)+gxa2(:,9 )*gmet(2,2)+gxa2(:,8 )*gmet(2,3)
884 196080 : r3a(:,5,2)=gxa2(:,5 )*gmet(2,1)+gxa2(:,4 )*gmet(2,2)+gxa2(:,3 )*gmet(2,3)
885 196080 : r3a(:,6,2)=gxa2(:,6 )*gmet(2,1)+gxa2(:,2 )*gmet(2,2)+gxa2(:,4 )*gmet(2,3)
886 196080 : r3a(:,1,3)=gxa2(:,1 )*gmet(3,1)+gxa2(:,6 )*gmet(3,2)+gxa2(:,5 )*gmet(3,3)
887 196080 : r3a(:,2,3)=gxa2(:,2 )*gmet(3,1)+gxa2(:,7 )*gmet(3,2)+gxa2(:,9 )*gmet(3,3)
888 : r3a(:,3,3)=gxa2(:,3 )*gmet(3,1)+gxa2(:,8 )*gmet(3,2)+gxa2(:,10)*gmet(3,3)
889 196080 : r3a(:,4,3)=gxa2(:,4 )*gmet(3,1)+gxa2(:,9 )*gmet(3,2)+gxa2(:,8 )*gmet(3,3)
890 196080 : r3a(:,5,3)=gxa2(:,5 )*gmet(3,1)+gxa2(:,4 )*gmet(3,2)+gxa2(:,3 )*gmet(3,3)
891 196080 : r3a(:,6,3)=gxa2(:,6 )*gmet(3,1)+gxa2(:,2 )*gmet(3,2)+gxa2(:,4 )*gmet(3,3)
892 :
893 : !Then compute r3(i,j,b)=r3a(b,l,j)*gmet(i,l)
894 : !stored as r3(re/im,(i,j),b)
895 196080 : r3(:,1,1)=r3a(:,1,1)*gmet(1,1)+r3a(:,6,1)*gmet(1,2)+r3a(:,5,1)*gmet(1,3)
896 196080 : r3(:,1,2)=r3a(:,6,1)*gmet(1,1)+r3a(:,2,1)*gmet(1,2)+r3a(:,4,1)*gmet(1,3)
897 : r3(:,1,3)=r3a(:,5,1)*gmet(1,1)+r3a(:,4,1)*gmet(1,2)+r3a(:,3,1)*gmet(1,3)
898 196080 : r3(:,2,1)=r3a(:,1,2)*gmet(2,1)+r3a(:,6,2)*gmet(2,2)+r3a(:,5,2)*gmet(2,3)
899 196080 : r3(:,2,2)=r3a(:,6,2)*gmet(2,1)+r3a(:,2,2)*gmet(2,2)+r3a(:,4,2)*gmet(2,3)
900 : r3(:,2,3)=r3a(:,5,2)*gmet(2,1)+r3a(:,4,2)*gmet(2,2)+r3a(:,3,2)*gmet(2,3)
901 196080 : r3(:,3,1)=r3a(:,1,3)*gmet(3,1)+r3a(:,6,3)*gmet(3,2)+r3a(:,5,3)*gmet(3,3)
902 196080 : r3(:,3,2)=r3a(:,6,3)*gmet(3,1)+r3a(:,2,3)*gmet(3,2)+r3a(:,4,3)*gmet(3,3)
903 : r3(:,3,3)=r3a(:,5,3)*gmet(3,1)+r3a(:,4,3)*gmet(3,2)+r3a(:,3,3)*gmet(3,3)
904 196080 : r3(:,4,1)=r3a(:,1,2)*gmet(3,1)+r3a(:,6,2)*gmet(3,2)+r3a(:,5,2)*gmet(3,3)
905 196080 : r3(:,4,2)=r3a(:,6,2)*gmet(3,1)+r3a(:,2,2)*gmet(3,2)+r3a(:,4,2)*gmet(3,3)
906 : r3(:,4,3)=r3a(:,5,2)*gmet(3,1)+r3a(:,4,2)*gmet(3,2)+r3a(:,3,2)*gmet(3,3)
907 196080 : r3(:,5,1)=r3a(:,1,1)*gmet(3,1)+r3a(:,6,1)*gmet(3,2)+r3a(:,5,1)*gmet(3,3)
908 196080 : r3(:,5,2)=r3a(:,6,1)*gmet(3,1)+r3a(:,2,1)*gmet(3,2)+r3a(:,4,1)*gmet(3,3)
909 : r3(:,5,3)=r3a(:,5,1)*gmet(3,1)+r3a(:,4,1)*gmet(3,2)+r3a(:,3,1)*gmet(3,3)
910 196080 : r3(:,6,1)=r3a(:,1,1)*gmet(2,1)+r3a(:,6,1)*gmet(2,2)+r3a(:,5,1)*gmet(2,3)
911 196080 : r3(:,6,2)=r3a(:,6,1)*gmet(2,1)+r3a(:,2,1)*gmet(2,2)+r3a(:,4,1)*gmet(2,3)
912 : r3(:,6,3)=r3a(:,5,1)*gmet(2,1)+r3a(:,4,1)*gmet(2,2)+r3a(:,3,1)*gmet(2,3)
913 :
914 : !Compute r11(a)=conjg(gxa1(a,l,m))*gmet(l,m)
915 : r11(:,1)=gxa1(:,1)*gmet(1,1)+gxa1(:,2)*gmet(2,2)+gxa1(:,3 )*gmet(3,3)&
916 : & +2.d0*(gxa1(:,4)*gmet(3,2)+gxa1(:,5)*gmet(3,1)+gxa1(:,6 )*gmet(2,1))
917 : r11(:,2)=gxa1(:,6)*gmet(1,1)+gxa1(:,7)*gmet(2,2)+gxa1(:,8 )*gmet(3,3)&
918 196080 : & +2.d0*(gxa1(:,9)*gmet(3,2)+gxa1(:,4)*gmet(3,1)+gxa1(:,2 )*gmet(2,1))
919 : r11(:,3)=gxa1(:,5)*gmet(1,1)+gxa1(:,9)*gmet(2,2)+gxa1(:,10)*gmet(3,3)&
920 196080 : & +2.d0*(gxa1(:,8)*gmet(3,2)+gxa1(:,3)*gmet(3,1)+gxa1(:,4 )*gmet(2,1))
921 65360 : r11(im,1)=-r11(im,1);r11(im,2)=-r11(im,2);r11(im,3)=-r11(im,3)
922 :
923 : !Compute r12(b)=gxa2(b,l,m)*gmet(l,m)
924 : r12(:,1)=gxa2(:,1)*gmet(1,1)+gxa2(:,2)*gmet(2,2)+gxa2(:,3 )*gmet(3,3)&
925 196080 : & +2.d0*(gxa2(:,4)*gmet(3,2)+gxa2(:,5)*gmet(3,1)+gxa2(:,6 )*gmet(2,1))
926 : r12(:,2)=gxa2(:,6)*gmet(1,1)+gxa2(:,7)*gmet(2,2)+gxa2(:,8 )*gmet(3,3)&
927 196080 : & +2.d0*(gxa2(:,9)*gmet(3,2)+gxa2(:,4)*gmet(3,1)+gxa2(:,2 )*gmet(2,1))
928 : r12(:,3)=gxa2(:,5)*gmet(1,1)+gxa2(:,9)*gmet(2,2)+gxa2(:,10)*gmet(3,3)&
929 : & +2.d0*(gxa2(:,8)*gmet(3,2)+gxa2(:,3)*gmet(3,1)+gxa2(:,4 )*gmet(2,1))
930 :
931 : !Finally compute rank2c(a,b)=7.5*conjg(gxa1(a,i,j))*r3(i,j,b) - 1.5*r11(a)*r12(b)
932 : !rank2c(re,1)=7.5d0*(gxa1(re,1 )*r3(re,1,1)+gxa1(im,1 )*r3(im,1,1)&
933 : !& +gxa1(re,2 )*r3(re,2,1)+gxa1(im,2 )*r3(im,2,1)&
934 : !& +gxa1(re,3 )*r3(re,3,1)+gxa1(im,3 )*r3(im,3,1))&
935 : !& +15.d0*(gxa1(re,4 )*r3(re,4,1)+gxa1(im,4 )*r3(im,4,1)&
936 : !& +gxa1(re,5 )*r3(re,5,1)+gxa1(im,5 )*r3(im,5,1)&
937 : !& +gxa1(re,6 )*r3(re,6,1)+gxa1(im,6 )*r3(im,6,1))&
938 : !& -1.5d0*(r11(re,1)*r12(re,1)-r11(im,1)*r12(im,1))
939 : !rank2c(re,2)=7.5d0*(gxa1(re,6 )*r3(re,1,2)+gxa1(im,6 )*r3(im,1,2)&
940 : !& +gxa1(re,7 )*r3(re,2,2)+gxa1(im,7 )*r3(im,2,2)&
941 : !& +gxa1(re,8 )*r3(re,3,2)+gxa1(im,8 )*r3(im,3,2))&
942 : !& +15.d0*(gxa1(re,9 )*r3(re,4,2)+gxa1(im,9 )*r3(im,4,2)&
943 : !& +gxa1(re,4 )*r3(re,5,2)+gxa1(im,4 )*r3(im,5,2)&
944 : !& +gxa1(re,2 )*r3(re,6,2)+gxa1(im,2 )*r3(im,6,2))&
945 : !& -1.5d0*(r11(re,2)*r12(re,2)-r11(im,2)*r12(im,2))
946 : !rank2c(re,3)=7.5d0*(gxa1(re,5 )*r3(re,1,3)+gxa1(im,5 )*r3(im,1,3)&
947 : !& +gxa1(re,9 )*r3(re,2,3)+gxa1(im,9 )*r3(im,2,3)&
948 : !& +gxa1(re,10)*r3(re,3,3)+gxa1(im,10)*r3(im,3,3))&
949 : !& +15.d0*(gxa1(re,8 )*r3(re,4,3)+gxa1(im,8 )*r3(im,4,3)&
950 : !& +gxa1(re,3 )*r3(re,5,3)+gxa1(im,3 )*r3(im,5,3)&
951 : !& +gxa1(re,4 )*r3(re,6,3)+gxa1(im,4 )*r3(im,6,3))&
952 : !& -1.5d0*(r11(re,3)*r12(re,3)-r11(im,3)*r12(im,3))
953 : rank2c(re,4)=7.5d0*(gxa1(re,5 )*r3(re,1,2)+gxa1(im,5 )*r3(im,1,2)&
954 : & +gxa1(re,9 )*r3(re,2,2)+gxa1(im,9 )*r3(im,2,2)&
955 : & +gxa1(re,10)*r3(re,3,2)+gxa1(im,10)*r3(im,3,2))&
956 : & +15.d0*(gxa1(re,8 )*r3(re,4,2)+gxa1(im,8 )*r3(im,4,2)&
957 : & +gxa1(re,3 )*r3(re,5,2)+gxa1(im,3 )*r3(im,5,2)&
958 : & +gxa1(re,4 )*r3(re,6,2)+gxa1(im,4 )*r3(im,6,2))&
959 65360 : & -1.5d0*(r11(re,3)*r12(re,2)-r11(im,3)*r12(im,2))
960 : rank2c(re,5)=7.5d0*(gxa1(re,5 )*r3(re,1,1)+gxa1(im,5 )*r3(im,1,1)&
961 : & +gxa1(re,9 )*r3(re,2,1)+gxa1(im,9 )*r3(im,2,1)&
962 : & +gxa1(re,10)*r3(re,3,1)+gxa1(im,10)*r3(im,3,1))&
963 : & +15.d0*(gxa1(re,8 )*r3(re,4,1)+gxa1(im,8 )*r3(im,4,1)&
964 : & +gxa1(re,3 )*r3(re,5,1)+gxa1(im,3 )*r3(im,5,1)&
965 : & +gxa1(re,4 )*r3(re,6,1)+gxa1(im,4 )*r3(im,6,1))&
966 65360 : & -1.5d0*(r11(re,3)*r12(re,1)-r11(im,3)*r12(im,1))
967 : rank2c(re,6)=7.5d0*(gxa1(re,6 )*r3(re,1,1)+gxa1(im,6 )*r3(im,1,1)&
968 : & +gxa1(re,7 )*r3(re,2,1)+gxa1(im,7 )*r3(im,2,1)&
969 : & +gxa1(re,8 )*r3(re,3,1)+gxa1(im,8 )*r3(im,3,1))&
970 : & +15.d0*(gxa1(re,9 )*r3(re,4,1)+gxa1(im,9 )*r3(im,4,1)&
971 : & +gxa1(re,4 )*r3(re,5,1)+gxa1(im,4 )*r3(im,5,1)&
972 : & +gxa1(re,2 )*r3(re,6,1)+gxa1(im,2 )*r3(im,6,1))&
973 65360 : & -1.5d0*(r11(re,2)*r12(re,1)-r11(im,2)*r12(im,1))
974 : !rank2c(im,1)=7.5d0*(gxa1(re,1 )*r3(im,1,1)-gxa1(im,1 )*r3(re,1,1)&
975 : !& +gxa1(re,2 )*r3(im,2,1)-gxa1(im,2 )*r3(re,2,1)&
976 : !& +gxa1(re,3 )*r3(im,3,1)-gxa1(im,3 )*r3(re,3,1))&
977 : !& +15.d0*(gxa1(re,4 )*r3(im,4,1)-gxa1(im,4 )*r3(re,4,1)&
978 : !& +gxa1(re,5 )*r3(im,5,1)-gxa1(im,5 )*r3(re,5,1)&
979 : !& +gxa1(re,6 )*r3(im,6,1)-gxa1(im,6 )*r3(re,6,1))&
980 : !& -1.5d0*(r11(re,1)*r12(im,1)+r11(im,1)*r12(re,1))
981 : !rank2c(im,2)=7.5d0*(gxa1(re,6 )*r3(im,1,2)-gxa1(im,6 )*r3(re,1,2)&
982 : !& +gxa1(re,7 )*r3(im,2,2)-gxa1(im,7 )*r3(re,2,2)&
983 : !& +gxa1(re,8 )*r3(im,3,2)-gxa1(im,8 )*r3(re,3,2))&
984 : !& +15.d0*(gxa1(re,9 )*r3(im,4,2)-gxa1(im,9 )*r3(re,4,2)&
985 : !& +gxa1(re,4 )*r3(im,5,2)-gxa1(im,4 )*r3(re,5,2)&
986 : !& +gxa1(re,2 )*r3(im,6,2)-gxa1(im,2 )*r3(re,6,2))&
987 : !& -1.5d0*(r11(re,2)*r12(im,2)+r11(im,2)*r12(re,2))
988 : !rank2c(im,3)=7.5d0*(gxa1(re,5 )*r3(im,1,3)-gxa1(im,5 )*r3(re,1,3)&
989 : !& +gxa1(re,9 )*r3(im,2,3)-gxa1(im,9 )*r3(re,2,3)&
990 : !& +gxa1(re,10)*r3(im,3,3)-gxa1(im,10)*r3(re,3,3))&
991 : !& +15.d0*(gxa1(re,8 )*r3(im,4,3)-gxa1(im,8 )*r3(re,4,3)&
992 : !& +gxa1(re,3 )*r3(im,5,3)-gxa1(im,3 )*r3(re,5,3)&
993 : !& +gxa1(re,4 )*r3(im,6,3)-gxa1(im,4 )*r3(re,6,3))&
994 : !& -1.5d0*(r11(re,3)*r12(im,3)+r11(im,3)*r12(re,3))
995 : rank2c(im,4)=7.5d0*(gxa1(re,5 )*r3(im,1,2)-gxa1(im,5 )*r3(re,1,2)&
996 : & +gxa1(re,9 )*r3(im,2,2)-gxa1(im,9 )*r3(re,2,2)&
997 : & +gxa1(re,10)*r3(im,3,2)-gxa1(im,10)*r3(re,3,2))&
998 : & +15.d0*(gxa1(re,8 )*r3(im,4,2)-gxa1(im,8 )*r3(re,4,2)&
999 : & +gxa1(re,3 )*r3(im,5,2)-gxa1(im,3 )*r3(re,5,2)&
1000 : & +gxa1(re,4 )*r3(im,6,2)-gxa1(im,4 )*r3(re,6,2))&
1001 65360 : & -1.5d0*(r11(re,3)*r12(im,2)+r11(im,3)*r12(re,2))
1002 : rank2c(im,5)=7.5d0*(gxa1(re,5 )*r3(im,1,1)-gxa1(im,5 )*r3(re,1,1)&
1003 : & +gxa1(re,9 )*r3(im,2,1)-gxa1(im,9 )*r3(re,2,1)&
1004 : & +gxa1(re,10)*r3(im,3,1)-gxa1(im,10)*r3(re,3,1))&
1005 : & +15.d0*(gxa1(re,8 )*r3(im,4,1)-gxa1(im,8 )*r3(re,4,1)&
1006 : & +gxa1(re,3 )*r3(im,5,1)-gxa1(im,3 )*r3(re,5,1)&
1007 : & +gxa1(re,4 )*r3(im,6,1)-gxa1(im,4 )*r3(re,6,1))&
1008 65360 : & -1.5d0*(r11(re,3)*r12(im,1)+r11(im,3)*r12(re,1))
1009 : rank2c(im,6)=7.5d0*(gxa1(re,6 )*r3(im,1,1)-gxa1(im,6 )*r3(re,1,1)&
1010 : & +gxa1(re,7 )*r3(im,2,1)-gxa1(im,7 )*r3(re,2,1)&
1011 : & +gxa1(re,8 )*r3(im,3,1)-gxa1(im,8 )*r3(re,3,1))&
1012 : & +15.d0*(gxa1(re,9 )*r3(im,4,1)-gxa1(im,9 )*r3(re,4,1)&
1013 : & +gxa1(re,4 )*r3(im,5,1)-gxa1(im,4 )*r3(re,5,1)&
1014 : & +gxa1(re,2 )*r3(im,6,1)-gxa1(im,2 )*r3(re,6,1))&
1015 65360 : & -1.5d0*(r11(re,2)*r12(im,1)+r11(im,2)*r12(re,1))
1016 :
1017 65360 : end subroutine cont33cso
1018 : !!***
1019 :
1020 :
1021 : !!****f* m_contract/cont33so
1022 : !! NAME
1023 : !! cont33so
1024 : !!
1025 : !! FUNCTION
1026 : !! Contract symmetric rank 3 tensor gxa1 with symmetric rank 3 tensor
1027 : !! gxa2 using metric tensor gmet and antisymmetric tensor amet to
1028 : !! produce rank 2 real tensor.
1029 : !!
1030 : !! INPUTS
1031 : !! gxa1(2,10)=rank 3 complex symmetric tensor
1032 : !! gxa2(2,10)=rank 3 complex symmetric tensor
1033 : !! gmet(3,3)=usual metric tensor (symmetric, real)
1034 : !! amet(2,3,3)=antisymmetric complex tensor used for spin-orbit
1035 : !!
1036 : !! OUTPUT
1037 : !! rank2(6)=rank 2 real tensor (pseudo-symmetric storage)
1038 : !!
1039 : !! NOTES
1040 : !! This contraction is used for spin-orbit correction in non-local
1041 : !! contribution to stresses.
1042 : !!
1043 : !! Symmetric gxa1, gxa2 are stored as
1044 : !! 111 221 331 321 311 211 222 332 322 333;
1045 : !! gmet(3,3) is symmetric but stored fully (9 elements);
1046 : !! amet(3,3) is antisymmetric but stored fully (9 elements);
1047 : !! Output rank2 is not symmetric but since
1048 : !! $rank2_{gxa1,gxa2}(a,b)=conjg(rank2_{gxa2,gxa1}(b,a))$
1049 : !! it is stored as 11 22 33 32 31 21.
1050 : !! Want 2*Re[contraction].
1051 : !!
1052 : !!{{\ \begin{equation}
1053 : !! rank2(a,b)=2 Re[15 r_{3A}(a,i,j) r_{3G}(b,j,i)-3 r_{3A}(a,b,i) r_1(i)]
1054 : !!\end{equation} }}
1055 : !! where:
1056 : !!{{\ \begin{eqnarray}
1057 : !! r_1(i) & = & gxa2(i,j,k) gmet(j,k) \nonumber
1058 : !! r_{3A}(i,j,k) & = & conjg(gxa1(p,i,j)) amet(p,k) \nonumber
1059 : !! r_{3G}(i,j,k) & = & gxa2(p,i,j) gmet(p,k)
1060 : !! \end{eqnarray} }}
1061 : !!
1062 : !! SOURCE
1063 :
1064 65360 : subroutine cont33so(gxa1,gxa2,gmet,amet,rank2)
1065 :
1066 : !Arguments ------------------------------------
1067 : !arrays
1068 : real(dp),intent(in) :: amet(2,3,3),gmet(3,3),gxa1(2,10),gxa2(2,10)
1069 : real(dp),intent(out) :: rank2(6)
1070 :
1071 : !Local variables-------------------------------
1072 : !scalars
1073 : integer,parameter :: im=2,re=1
1074 : integer :: ii
1075 : !arrays
1076 : real(dp) :: r1(2,3),r3A(2,18),r3G(2,18),s31(6),s33(6)
1077 :
1078 : ! *************************************************************************
1079 :
1080 : !Compute r1(i)=gxa2(i,j,k)*gmet(j,k)
1081 :
1082 196080 : do ii=1,2
1083 : r1(ii,1)=gmet(1,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,2)+&
1084 : & gmet(3,3)*gxa2(ii,3)+2.d0*(&
1085 : & gmet(3,2)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,5)+&
1086 130720 : & gmet(2,1)*gxa2(ii,6))
1087 : r1(ii,2)=gmet(1,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,7)+&
1088 : & gmet(3,3)*gxa2(ii,8)+2.d0*(&
1089 : & gmet(3,2)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,4)+&
1090 130720 : & gmet(2,1)*gxa2(ii,2))
1091 : r1(ii,3)=gmet(1,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,9)+&
1092 : & gmet(3,3)*gxa2(ii,10)+2.d0*(&
1093 : & gmet(3,2)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,3)+&
1094 196080 : & gmet(2,1)*gxa2(ii,4))
1095 : end do
1096 :
1097 : !Compute r3G(i,j,k)=gxa2(p,i,j)*gmet(p,k)
1098 : !Write out components for 18 distinct terms, Re and Im
1099 : !(r3G is symmetric in first two indices, not in all permutations)
1100 : !Store as 111 221 331 321 311 211
1101 : !112 222 332 322 312 212
1102 : !113 223 333 323 313 213
1103 196080 : do ii=1,2
1104 130720 : r3G(ii, 1)=gmet(1,1)*gxa2(ii,1)+gmet(2,1)*gxa2(ii,6)+gmet(3,1)*gxa2(ii,5)
1105 130720 : r3G(ii, 2)=gmet(1,1)*gxa2(ii,2)+gmet(2,1)*gxa2(ii,7)+gmet(3,1)*gxa2(ii,9)
1106 130720 : r3G(ii, 3)=gmet(1,1)*gxa2(ii,3)+gmet(2,1)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,10)
1107 130720 : r3G(ii, 4)=gmet(1,1)*gxa2(ii,4)+gmet(2,1)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,8)
1108 130720 : r3G(ii, 5)=gmet(1,1)*gxa2(ii,5)+gmet(2,1)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,3)
1109 130720 : r3G(ii, 6)=gmet(1,1)*gxa2(ii,6)+gmet(2,1)*gxa2(ii,2)+gmet(3,1)*gxa2(ii,4)
1110 130720 : r3G(ii, 7)=gmet(2,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,5)
1111 130720 : r3G(ii, 8)=gmet(2,1)*gxa2(ii,2)+gmet(2,2)*gxa2(ii,7)+gmet(3,2)*gxa2(ii,9)
1112 130720 : r3G(ii, 9)=gmet(2,1)*gxa2(ii,3)+gmet(2,2)*gxa2(ii,8)+gmet(3,2)*gxa2(ii,10)
1113 130720 : r3G(ii,10)=gmet(2,1)*gxa2(ii,4)+gmet(2,2)*gxa2(ii,9)+gmet(3,2)*gxa2(ii,8)
1114 130720 : r3G(ii,11)=gmet(2,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,3)
1115 130720 : r3G(ii,12)=gmet(2,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,4)
1116 130720 : r3G(ii,13)=gmet(3,1)*gxa2(ii,1)+gmet(3,2)*gxa2(ii,6)+gmet(3,3)*gxa2(ii,5)
1117 130720 : r3G(ii,14)=gmet(3,1)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,7)+gmet(3,3)*gxa2(ii,9)
1118 130720 : r3G(ii,15)=gmet(3,1)*gxa2(ii,3)+gmet(3,2)*gxa2(ii,8)+gmet(3,3)*gxa2(ii,10)
1119 130720 : r3G(ii,16)=gmet(3,1)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,9)+gmet(3,3)*gxa2(ii,8)
1120 130720 : r3G(ii,17)=gmet(3,1)*gxa2(ii,5)+gmet(3,2)*gxa2(ii,4)+gmet(3,3)*gxa2(ii,3)
1121 196080 : r3G(ii,18)=gmet(3,1)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,2)+gmet(3,3)*gxa2(ii,4)
1122 : end do
1123 :
1124 : !Compute r3A(i,j,k)=conjg(gxa1(p,i,j))*amet(p,k)
1125 : !Write out components for 18 distinct terms, Re and Im
1126 : !(r3A is symmetric in first two indices, not in all permutations)
1127 : !Store as 111 221 331 321 311 211
1128 : !112 222 332 322 312 212
1129 : !113 223 333 323 313 213
1130 : !Note that, since amet is antisymmetric, amet(i,i)=0
1131 :
1132 : r3A(re, 1)=amet(re,2,1)*gxa1(re,6 )+amet(im,2,1)*gxa1(im,6 )+&
1133 65360 : & amet(re,3,1)*gxa1(re,5 )+amet(im,3,1)*gxa1(im,5 )
1134 : r3A(re, 2)=amet(re,2,1)*gxa1(re,7 )+amet(im,2,1)*gxa1(im,7 )+&
1135 65360 : & amet(re,3,1)*gxa1(re,9 )+amet(im,3,1)*gxa1(im,9 )
1136 : r3A(re, 3)=amet(re,2,1)*gxa1(re,8 )+amet(im,2,1)*gxa1(im,8 )+&
1137 65360 : & amet(re,3,1)*gxa1(re,10)+amet(im,3,1)*gxa1(im,10)
1138 : r3A(re, 4)=amet(re,2,1)*gxa1(re,9 )+amet(im,2,1)*gxa1(im,9 )+&
1139 65360 : & amet(re,3,1)*gxa1(re,8 )+amet(im,3,1)*gxa1(im,8 )
1140 : r3A(re, 5)=amet(re,2,1)*gxa1(re,4 )+amet(im,2,1)*gxa1(im,4 )+&
1141 65360 : & amet(re,3,1)*gxa1(re,3 )+amet(im,3,1)*gxa1(im,3 )
1142 : r3A(re, 6)=amet(re,2,1)*gxa1(re,2 )+amet(im,2,1)*gxa1(im,2 )+&
1143 65360 : & amet(re,3,1)*gxa1(re,4 )+amet(im,3,1)*gxa1(im,4 )
1144 : r3A(re, 7)=amet(re,1,2)*gxa1(re,1 )+amet(im,1,2)*gxa1(im,1 )+&
1145 65360 : & amet(re,3,2)*gxa1(re,5 )+amet(im,3,2)*gxa1(im,5 )
1146 : r3A(re, 8)=amet(re,1,2)*gxa1(re,2 )+amet(im,1,2)*gxa1(im,2 )+&
1147 65360 : & amet(re,3,2)*gxa1(re,9 )+amet(im,3,2)*gxa1(im,9 )
1148 : r3A(re, 9)=amet(re,1,2)*gxa1(re,3 )+amet(im,1,2)*gxa1(im,3 )+&
1149 65360 : & amet(re,3,2)*gxa1(re,10)+amet(im,3,2)*gxa1(im,10)
1150 : r3A(re,10)=amet(re,1,2)*gxa1(re,4 )+amet(im,1,2)*gxa1(im,4 )+&
1151 65360 : & amet(re,3,2)*gxa1(re,8 )+amet(im,3,2)*gxa1(im,8 )
1152 : r3A(re,11)=amet(re,1,2)*gxa1(re,5 )+amet(im,1,2)*gxa1(im,5 )+&
1153 65360 : & amet(re,3,2)*gxa1(re,3 )+amet(im,3,2)*gxa1(im,3 )
1154 : r3A(re,12)=amet(re,1,2)*gxa1(re,6 )+amet(im,1,2)*gxa1(im,6 )+&
1155 65360 : & amet(re,3,2)*gxa1(re,4 )+amet(im,3,2)*gxa1(im,4 )
1156 : r3A(re,13)=amet(re,1,3)*gxa1(re,1 )+amet(im,1,3)*gxa1(im,1 )+&
1157 65360 : & amet(re,2,3)*gxa1(re,6 )+amet(im,2,3)*gxa1(im,6 )
1158 : r3A(re,14)=amet(re,1,3)*gxa1(re,2 )+amet(im,1,3)*gxa1(im,2 )+&
1159 65360 : & amet(re,2,3)*gxa1(re,7 )+amet(im,2,3)*gxa1(im,7 )
1160 : r3A(re,15)=amet(re,1,3)*gxa1(re,3 )+amet(im,1,3)*gxa1(im,3 )+&
1161 65360 : & amet(re,2,3)*gxa1(re,8 )+amet(im,2,3)*gxa1(im,8 )
1162 : r3A(re,16)=amet(re,1,3)*gxa1(re,4 )+amet(im,1,3)*gxa1(im,4 )+&
1163 65360 : & amet(re,2,3)*gxa1(re,9 )+amet(im,2,3)*gxa1(im,9 )
1164 : r3A(re,17)=amet(re,1,3)*gxa1(re,5 )+amet(im,1,3)*gxa1(im,5 )+&
1165 65360 : & amet(re,2,3)*gxa1(re,4 )+amet(im,2,3)*gxa1(im,4 )
1166 : r3A(re,18)=amet(re,1,3)*gxa1(re,6 )+amet(im,1,3)*gxa1(im,6 )+&
1167 65360 : & amet(re,2,3)*gxa1(re,2 )+amet(im,2,3)*gxa1(im,2 )
1168 :
1169 : r3A(im, 1)=amet(im,2,1)*gxa1(re,6 )-amet(re,2,1)*gxa1(im,6 )+&
1170 65360 : & amet(im,3,1)*gxa1(re,5 )-amet(re,3,1)*gxa1(im,5 )
1171 : r3A(im, 2)=amet(im,2,1)*gxa1(re,7 )-amet(re,2,1)*gxa1(im,7 )+&
1172 65360 : & amet(im,3,1)*gxa1(re,9 )-amet(re,3,1)*gxa1(im,9 )
1173 : r3A(im, 3)=amet(im,2,1)*gxa1(re,8 )-amet(re,2,1)*gxa1(im,8 )+&
1174 65360 : & amet(im,3,1)*gxa1(re,10)-amet(re,3,1)*gxa1(im,10)
1175 : r3A(im, 4)=amet(im,2,1)*gxa1(re,9 )-amet(re,2,1)*gxa1(im,9 )+&
1176 65360 : & amet(im,3,1)*gxa1(re,8 )-amet(re,3,1)*gxa1(im,8 )
1177 : r3A(im, 5)=amet(im,2,1)*gxa1(re,4 )-amet(re,2,1)*gxa1(im,4 )+&
1178 65360 : & amet(im,3,1)*gxa1(re,3 )-amet(re,3,1)*gxa1(im,3 )
1179 : r3A(im, 6)=amet(im,2,1)*gxa1(re,2 )-amet(re,2,1)*gxa1(im,2 )+&
1180 65360 : & amet(im,3,1)*gxa1(re,4 )-amet(re,3,1)*gxa1(im,4 )
1181 : r3A(im, 7)=amet(im,1,2)*gxa1(re,1 )-amet(re,1,2)*gxa1(im,1 )+&
1182 65360 : & amet(im,3,2)*gxa1(re,5 )-amet(re,3,2)*gxa1(im,5 )
1183 : r3A(im, 8)=amet(im,1,2)*gxa1(re,2 )-amet(re,1,2)*gxa1(im,2 )+&
1184 65360 : & amet(im,3,2)*gxa1(re,9 )-amet(re,3,2)*gxa1(im,9 )
1185 : r3A(im, 9)=amet(im,1,2)*gxa1(re,3 )-amet(re,1,2)*gxa1(im,3 )+&
1186 65360 : & amet(im,3,2)*gxa1(re,10)-amet(re,3,2)*gxa1(im,10)
1187 : r3A(im,10)=amet(im,1,2)*gxa1(re,4 )-amet(re,1,2)*gxa1(im,4 )+&
1188 65360 : & amet(im,3,2)*gxa1(re,8 )-amet(re,3,2)*gxa1(im,8 )
1189 : r3A(im,11)=amet(im,1,2)*gxa1(re,5 )-amet(re,1,2)*gxa1(im,5 )+&
1190 65360 : & amet(im,3,2)*gxa1(re,3 )-amet(re,3,2)*gxa1(im,3 )
1191 : r3A(im,12)=amet(im,1,2)*gxa1(re,6 )-amet(re,1,2)*gxa1(im,6 )+&
1192 65360 : & amet(im,3,2)*gxa1(re,4 )-amet(re,3,2)*gxa1(im,4 )
1193 : r3A(im,13)=amet(im,1,3)*gxa1(re,1 )-amet(re,1,3)*gxa1(im,1 )+&
1194 65360 : & amet(im,2,3)*gxa1(re,6 )-amet(re,2,3)*gxa1(im,6 )
1195 : r3A(im,14)=amet(im,1,3)*gxa1(re,2 )-amet(re,1,3)*gxa1(im,2 )+&
1196 65360 : & amet(im,2,3)*gxa1(re,7 )-amet(re,2,3)*gxa1(im,7 )
1197 : r3A(im,15)=amet(im,1,3)*gxa1(re,3 )-amet(re,1,3)*gxa1(im,3 )+&
1198 65360 : & amet(im,2,3)*gxa1(re,8 )-amet(re,2,3)*gxa1(im,8 )
1199 : r3A(im,16)=amet(im,1,3)*gxa1(re,4 )-amet(re,1,3)*gxa1(im,4 )+&
1200 65360 : & amet(im,2,3)*gxa1(re,9 )-amet(re,2,3)*gxa1(im,9 )
1201 : r3A(im,17)=amet(im,1,3)*gxa1(re,5 )-amet(re,1,3)*gxa1(im,5 )+&
1202 65360 : & amet(im,2,3)*gxa1(re,4 )-amet(re,2,3)*gxa1(im,4 )
1203 : r3A(im,18)=amet(im,1,3)*gxa1(re,6 )-amet(re,1,3)*gxa1(im,6 )+&
1204 65360 : & amet(im,2,3)*gxa1(re,2 )-amet(re,2,3)*gxa1(im,2 )
1205 :
1206 : !Compute s33(a,b)=2*Re[r3A(a,i,j)*r3G(b,j,i)]
1207 :
1208 : s33(1)=2.d0*(r3A(re, 1)*r3G(re, 1)-r3A(im, 1)*r3G(im, 1)+&
1209 : & r3A(re,12)*r3G(re,12)-r3A(im,12)*r3G(im,12)+&
1210 : & r3A(re,17)*r3G(re,17)-r3A(im,17)*r3G(im,17)+&
1211 : & r3A(re,11)*r3G(re,18)-r3A(im,11)*r3G(im,18)+&
1212 : & r3A(re,18)*r3G(re,11)-r3A(im,18)*r3G(im,11)+&
1213 : & r3A(re, 5)*r3G(re,13)-r3A(im, 5)*r3G(im,13)+&
1214 : & r3A(re,13)*r3G(re, 5)-r3A(im,13)*r3G(im, 5)+&
1215 : & r3A(re, 6)*r3G(re, 7)-r3A(im, 6)*r3G(im, 7)+&
1216 65360 : & r3A(re, 7)*r3G(re, 6)-r3A(im, 7)*r3G(im, 6))
1217 : s33(2)=2.d0*(r3A(re, 6)*r3G(re, 6)-r3A(im, 6)*r3G(im, 6)+&
1218 : & r3A(re, 8)*r3G(re, 8)-r3A(im, 8)*r3G(im, 8)+&
1219 : & r3A(re,16)*r3G(re,16)-r3A(im,16)*r3G(im,16)+&
1220 : & r3A(re,10)*r3G(re,14)-r3A(im,10)*r3G(im,14)+&
1221 : & r3A(re,14)*r3G(re,10)-r3A(im,14)*r3G(im,10)+&
1222 : & r3A(re, 4)*r3G(re,18)-r3A(im, 4)*r3G(im,18)+&
1223 : & r3A(re,18)*r3G(re, 4)-r3A(im,18)*r3G(im, 4)+&
1224 : & r3A(re, 2)*r3G(re,12)-r3A(im, 2)*r3G(im,12)+&
1225 65360 : & r3A(re,12)*r3G(re, 2)-r3A(im,12)*r3G(im, 2))
1226 : s33(3)=2.d0*(r3A(re, 5)*r3G(re, 5)-r3A(im, 5)*r3G(im, 5)+&
1227 : & r3A(re,10)*r3G(re,10)-r3A(im,10)*r3G(im,10)+&
1228 : & r3A(re,15)*r3G(re,15)-r3A(im,15)*r3G(im,15)+&
1229 : & r3A(re, 9)*r3G(re,16)-r3A(im, 9)*r3G(im,16)+&
1230 : & r3A(re,16)*r3G(re, 9)-r3A(im,16)*r3G(im, 9)+&
1231 : & r3A(re, 3)*r3G(re,17)-r3A(im, 3)*r3G(im,17)+&
1232 : & r3A(re,17)*r3G(re, 3)-r3A(im,17)*r3G(im, 3)+&
1233 : & r3A(re, 4)*r3G(re,11)-r3A(im, 4)*r3G(im,11)+&
1234 65360 : & r3A(re,11)*r3G(re, 4)-r3A(im,11)*r3G(im, 4))
1235 : s33(4)=2.d0*(r3A(re, 5)*r3G(re, 6)-r3A(im, 5)*r3G(im, 6)+&
1236 : & r3A(re,10)*r3G(re, 8)-r3A(im,10)*r3G(im, 8)+&
1237 : & r3A(re,15)*r3G(re,16)-r3A(im,15)*r3G(im,16)+&
1238 : & r3A(re, 9)*r3G(re,14)-r3A(im, 9)*r3G(im,14)+&
1239 : & r3A(re,16)*r3G(re,10)-r3A(im,16)*r3G(im,10)+&
1240 : & r3A(re, 3)*r3G(re,18)-r3A(im, 3)*r3G(im,18)+&
1241 : & r3A(re,17)*r3G(re, 4)-r3A(im,17)*r3G(im, 4)+&
1242 : & r3A(re, 4)*r3G(re,12)-r3A(im, 4)*r3G(im,12)+&
1243 65360 : & r3A(re,11)*r3G(re, 2)-r3A(im,11)*r3G(im, 2))
1244 : s33(5)=2.d0*(r3A(re, 5)*r3G(re, 1)-r3A(im, 5)*r3G(im, 1)+&
1245 : & r3A(re,10)*r3G(re,12)-r3A(im,10)*r3G(im,12)+&
1246 : & r3A(re,15)*r3G(re,17)-r3A(im,15)*r3G(im,17)+&
1247 : & r3A(re, 9)*r3G(re,18)-r3A(im, 9)*r3G(im,18)+&
1248 : & r3A(re,16)*r3G(re,11)-r3A(im,16)*r3G(im,11)+&
1249 : & r3A(re, 3)*r3G(re,13)-r3A(im, 3)*r3G(im,13)+&
1250 : & r3A(re,17)*r3G(re, 5)-r3A(im,17)*r3G(im, 5)+&
1251 : & r3A(re, 4)*r3G(re, 7)-r3A(im, 4)*r3G(im, 7)+&
1252 65360 : & r3A(re,11)*r3G(re, 6)-r3A(im,11)*r3G(im, 6))
1253 : s33(6)=2.d0*(r3A(re, 6)*r3G(re, 1)-r3A(im, 6)*r3G(im, 1)+&
1254 : & r3A(re, 8)*r3G(re,12)-r3A(im, 8)*r3G(im,12)+&
1255 : & r3A(re,16)*r3G(re,17)-r3A(im,16)*r3G(im,17)+&
1256 : & r3A(re,10)*r3G(re,18)-r3A(im,10)*r3G(im,18)+&
1257 : & r3A(re,14)*r3G(re,11)-r3A(im,14)*r3G(im,11)+&
1258 : & r3A(re, 4)*r3G(re,13)-r3A(im, 4)*r3G(im,13)+&
1259 : & r3A(re,18)*r3G(re, 5)-r3A(im,18)*r3G(im, 5)+&
1260 : & r3A(re, 2)*r3G(re, 7)-r3A(im, 2)*r3G(im, 7)+&
1261 65360 : & r3A(re,12)*r3G(re, 6)-r3A(im,12)*r3G(im, 6))
1262 :
1263 : !Compute s31(a,b)=2*Re[r3A(a,b,i)*r1(i)]
1264 :
1265 : s31(1)=2.d0*(r1(re,1)*r3A(re, 1)-r1(im,1)*r3A(im, 1)+&
1266 : & r1(re,2)*r3A(re, 7)-r1(im,2)*r3A(im, 7)+&
1267 65360 : & r1(re,3)*r3A(re,13)-r1(im,3)*r3A(im,13))
1268 : s31(2)=2.d0*(r1(re,1)*r3A(re, 2)-r1(im,1)*r3A(im, 2)+&
1269 : & r1(re,2)*r3A(re, 8)-r1(im,2)*r3A(im, 8)+&
1270 65360 : & r1(re,3)*r3A(re,14)-r1(im,3)*r3A(im,14))
1271 : s31(3)=2.d0*(r1(re,1)*r3A(re, 3)-r1(im,1)*r3A(im, 3)+&
1272 : & r1(re,2)*r3A(re, 9)-r1(im,2)*r3A(im, 9)+&
1273 65360 : & r1(re,3)*r3A(re,15)-r1(im,3)*r3A(im,15))
1274 : s31(4)=2.d0*(r1(re,1)*r3A(re, 4)-r1(im,1)*r3A(im, 4)+&
1275 : & r1(re,2)*r3A(re,10)-r1(im,2)*r3A(im,10)+&
1276 65360 : & r1(re,3)*r3A(re,16)-r1(im,3)*r3A(im,16))
1277 : s31(5)=2.d0*(r1(re,1)*r3A(re, 5)-r1(im,1)*r3A(im, 5)+&
1278 : & r1(re,2)*r3A(re,11)-r1(im,2)*r3A(im,11)+&
1279 65360 : & r1(re,3)*r3A(re,17)-r1(im,3)*r3A(im,17))
1280 : s31(6)=2.d0*(r1(re,1)*r3A(re, 6)-r1(im,1)*r3A(im, 6)+&
1281 : & r1(re,2)*r3A(re,12)-r1(im,2)*r3A(im,12)+&
1282 65360 : & r1(re,3)*r3A(re,18)-r1(im,3)*r3A(im,18))
1283 :
1284 : !Finally, compute rank2(a,b)=-15*s33(a,b)+3*s31(a,b)
1285 :
1286 457520 : rank2(:)=15.d0*s33(:)-3.d0*s31(:)
1287 :
1288 65360 : end subroutine cont33so
1289 : !!***
1290 :
1291 : !!****f* m_contract/cont35
1292 : !! NAME
1293 : !! cont35
1294 : !!
1295 : !! FUNCTION
1296 : !! Contract symmetric rank3 tensor gxa with rank5 symmetric tensor to
1297 : !! produce symmetric rank2 tensor.
1298 : !!
1299 : !! INPUTS
1300 : !! gxa(2,10)=rank 3 symmetric complex tensor in order
1301 : !! rank5(2,21)=rank 5 complex tensor (symmetric storage)
1302 : !!
1303 : !! OUTPUT
1304 : !! rank2(6)=rank 2 real tensor (symmetric storage) 11 22 33 32 31 21.
1305 : !!
1306 : !! NOTES
1307 : !! Tensors are in "symmetric" storage mode.
1308 : !! For rank 3 tensor gxa this is
1309 : !! 111 221 331 321 311 211 222 332 322 333;
1310 : !! For rank 5 tensor rank5 this is
1311 : !! 11111 22111 33111 32111 31111 21111 22211 33211 32211 33311
1312 : !! 22221 33221 32221 33321 33331 22222 33222 32222 33322 33332 33333;
1313 : !! For rank 2 tensor rank2 this is 11, 22, 33, 32, 31, 21;
1314 : !! gxa and rank5 are complex; rank2 is real.
1315 : !! Want $2 Re[contraction]$.
1316 : !! $rank2(a,b)=2 Re[gxa(i,j,k)^"*" rank5(a,b,i,j,k)]$.
1317 : !!
1318 : !! Note that the input gxa is typically the result of
1319 : !!{{\ \begin{equation}
1320 : !! gxa(i,j,k)=[{5 \over 2} gmet(i,l) gmet(j,m) gmet(k,n) - {3 \over 2} gmet(i,j) gmet(l,m) gmet(k,n)] gxa_old(l,m)
1321 : !!\end{equation} }}
1322 : !! where the subroutine "metcon" already includes weights in the definition
1323 : !! of gxa for off-diagonal elements.
1324 : !!
1325 : !! SOURCE
1326 :
1327 70046 : subroutine cont35(gxa,rank5,rank2)
1328 :
1329 : !Arguments ------------------------------------
1330 : !arrays
1331 : real(dp),intent(in) :: gxa(2,10),rank5(2,21)
1332 : real(dp),intent(out) :: rank2(6)
1333 :
1334 : !Local variables-------------------------------
1335 : !scalars
1336 : integer,parameter :: im=2,re=1
1337 :
1338 : ! *************************************************************************
1339 :
1340 : !Simply write out index summations
1341 :
1342 : !a=1, b=1 in rank2(a,b) --> maps to index 1
1343 : rank2(1)=2.0d0*(&
1344 : & gxa(re, 1)*rank5(re, 1)+gxa(im, 1)*rank5(im, 1)+&
1345 : & gxa(re, 7)*rank5(re, 7)+gxa(im, 7)*rank5(im, 7)+&
1346 : & gxa(re,10)*rank5(re,10)+gxa(im,10)*rank5(im,10)+&
1347 : & gxa(re, 2)*rank5(re, 2)+gxa(im, 2)*rank5(im, 2)+&
1348 : & gxa(re, 3)*rank5(re, 3)+gxa(im, 3)*rank5(im, 3)+&
1349 : & gxa(re, 5)*rank5(re, 5)+gxa(im, 5)*rank5(im, 5)+&
1350 : & gxa(re, 6)*rank5(re, 6)+gxa(im, 6)*rank5(im, 6)+&
1351 : & gxa(re, 8)*rank5(re, 8)+gxa(im, 8)*rank5(im, 8)+&
1352 : & gxa(re, 9)*rank5(re, 9)+gxa(im, 9)*rank5(im, 9)+&
1353 70046 : & gxa(re, 4)*rank5(re, 4)+gxa(im, 4)*rank5(im, 4))
1354 :
1355 :
1356 : !a=2, b=2 in rank2(a,b) --> maps to index 2
1357 : rank2(2)=2.0d0*(&
1358 : & gxa(re, 1)*rank5(re, 2)+gxa(im, 1)*rank5(im, 2)+&
1359 : & gxa(re, 7)*rank5(re,16)+gxa(im, 7)*rank5(im,16)+&
1360 : & gxa(re,10)*rank5(re,19)+gxa(im,10)*rank5(im,19)+&
1361 : & gxa(re, 2)*rank5(re,11)+gxa(im, 2)*rank5(im,11)+&
1362 : & gxa(re, 3)*rank5(re,12)+gxa(im, 3)*rank5(im,12)+&
1363 : & gxa(re, 5)*rank5(re, 9)+gxa(im, 5)*rank5(im, 9)+&
1364 : & gxa(re, 6)*rank5(re, 7)+gxa(im, 6)*rank5(im, 7)+&
1365 : & gxa(re, 8)*rank5(re,17)+gxa(im, 8)*rank5(im,17)+&
1366 : & gxa(re, 9)*rank5(re,18)+gxa(im, 9)*rank5(im,18)+&
1367 70046 : & gxa(re, 4)*rank5(re,13)+gxa(im, 4)*rank5(im,13))
1368 :
1369 : !a=3, b=3 in rank2(a,b) --> maps to index 3
1370 : rank2(3)=2.0d0*(&
1371 : & gxa(re, 1)*rank5(re, 3)+gxa(im, 1)*rank5(im, 3)+&
1372 : & gxa(re, 7)*rank5(re,17)+gxa(im, 7)*rank5(im,17)+&
1373 : & gxa(re,10)*rank5(re,21)+gxa(im,10)*rank5(im,21)+&
1374 : & gxa(re, 2)*rank5(re,12)+gxa(im, 2)*rank5(im,12)+&
1375 : & gxa(re, 3)*rank5(re,15)+gxa(im, 3)*rank5(im,15)+&
1376 : & gxa(re, 5)*rank5(re,10)+gxa(im, 5)*rank5(im,10)+&
1377 : & gxa(re, 6)*rank5(re, 8)+gxa(im, 6)*rank5(im, 8)+&
1378 : & gxa(re, 8)*rank5(re,20)+gxa(im, 8)*rank5(im,20)+&
1379 : & gxa(re, 9)*rank5(re,19)+gxa(im, 9)*rank5(im,19)+&
1380 70046 : & gxa(re, 4)*rank5(re,14)+gxa(im, 4)*rank5(im,14))
1381 :
1382 : !a=3, b=2 in rank2(a,b) --> maps to index 4
1383 : rank2(4)=2.0d0*(&
1384 : & gxa(re, 1)*rank5(re, 4)+gxa(im, 1)*rank5(im, 4)+&
1385 : & gxa(re, 7)*rank5(re,18)+gxa(im, 7)*rank5(im,18)+&
1386 : & gxa(re,10)*rank5(re,20)+gxa(im,10)*rank5(im,20)+&
1387 : & gxa(re, 2)*rank5(re,13)+gxa(im, 2)*rank5(im,13)+&
1388 : & gxa(re, 3)*rank5(re,14)+gxa(im, 3)*rank5(im,14)+&
1389 : & gxa(re, 5)*rank5(re, 8)+gxa(im, 5)*rank5(im, 8)+&
1390 : & gxa(re, 6)*rank5(re, 9)+gxa(im, 6)*rank5(im, 9)+&
1391 : & gxa(re, 8)*rank5(re,19)+gxa(im, 8)*rank5(im,19)+&
1392 : & gxa(re, 9)*rank5(re,17)+gxa(im, 9)*rank5(im,17)+&
1393 70046 : & gxa(re, 4)*rank5(re,12)+gxa(im, 4)*rank5(im,12))
1394 :
1395 : !a=3, b=1 in rank2(a,b) --> maps to index 5
1396 : rank2(5)=2.0d0*(&
1397 : & gxa(re, 1)*rank5(re, 5)+gxa(im, 1)*rank5(im, 5)+&
1398 : & gxa(re, 7)*rank5(re,13)+gxa(im, 7)*rank5(im,13)+&
1399 : & gxa(re,10)*rank5(re,15)+gxa(im,10)*rank5(im,15)+&
1400 : & gxa(re, 2)*rank5(re, 9)+gxa(im, 2)*rank5(im, 9)+&
1401 : & gxa(re, 3)*rank5(re,10)+gxa(im, 3)*rank5(im,10)+&
1402 : & gxa(re, 5)*rank5(re, 3)+gxa(im, 5)*rank5(im, 3)+&
1403 : & gxa(re, 6)*rank5(re, 4)+gxa(im, 6)*rank5(im, 4)+&
1404 : & gxa(re, 8)*rank5(re,14)+gxa(im, 8)*rank5(im,14)+&
1405 : & gxa(re, 9)*rank5(re,12)+gxa(im, 9)*rank5(im,12)+&
1406 70046 : & gxa(re, 4)*rank5(re, 8)+gxa(im, 4)*rank5(im, 8))
1407 :
1408 : !a=2, b=1 in rank2(a,b) --> maps to index 6
1409 : rank2(6)=2.0d0*(&
1410 : & gxa(re, 1)*rank5(re, 6)+gxa(im, 1)*rank5(im, 6)+&
1411 : & gxa(re, 7)*rank5(re,11)+gxa(im, 7)*rank5(im,11)+&
1412 : & gxa(re,10)*rank5(re,14)+gxa(im,10)*rank5(im,14)+&
1413 : & gxa(re, 2)*rank5(re, 7)+gxa(im, 2)*rank5(im, 7)+&
1414 : & gxa(re, 3)*rank5(re, 8)+gxa(im, 3)*rank5(im, 8)+&
1415 : & gxa(re, 5)*rank5(re, 4)+gxa(im, 5)*rank5(im, 4)+&
1416 : & gxa(re, 6)*rank5(re, 2)+gxa(im, 6)*rank5(im, 2)+&
1417 : & gxa(re, 8)*rank5(re,12)+gxa(im, 8)*rank5(im,12)+&
1418 : & gxa(re, 9)*rank5(re,13)+gxa(im, 9)*rank5(im,13)+&
1419 70046 : & gxa(re, 4)*rank5(re, 9)+gxa(im, 4)*rank5(im, 9))
1420 :
1421 70046 : end subroutine cont35
1422 : !!***
1423 :
1424 : !!****f* m_contract/metcon
1425 : !! NAME
1426 : !! metcon
1427 : !!
1428 : !! FUNCTION
1429 : !! Carries out specialized metric tensor contractions needed for
1430 : !! l=0,1,2,3 nonlocal Kleinman-Bylander pseudopotential operation.
1431 : !! Full advantage is taken of the full permutational symmetry of these
1432 : !! tensors.
1433 : !!
1434 : !! INPUTS
1435 : !! rank=0,1,2, or 3 = rank of input tensor aa
1436 : !! gmet(3,3)=metric tensor (array is symmetric but stored as 3x3)
1437 : !! aa(2,(rank+1)*(rank+2)/2)=unique elements of complex input tensor
1438 : !!
1439 : !! OUTPUT
1440 : !! bb(2,(rank+1)*(rank+2)/2)=unique elements of complex output tensor
1441 : !!
1442 : !! NOTES
1443 : !! All tensors are stored in a compressed storage mode defined below;
1444 : !! input and output conform to this scheme.
1445 : !! When tensor elements occur repeatedly due to symmetry, the
1446 : !! WEIGHT IS INCLUDED in the output tensor element to simplify later
1447 : !! contractions with other tensors of the same rank and form, i.e. the
1448 : !! next contraction is then simply a dot product over the unique elements.
1449 : !!
1450 : !! Definitions of the contractions:
1451 : !!
1452 : !! rank=0: bb = aa (simply copy the scalar value, no use of gmet)
1453 : !!
1454 : !! rank=1: bb(i)= $gmet(i,l) aa(l)$ (3 elements in, 3 elements out)
1455 : !!
1456 : !! rank=2: bb(i,j)= $[{3 \over 2} gmet(i,l) gmet(j,m) - {1 \over 2} gmet(i,j) gmet(l,m)] aa(l,m)$
1457 : !! (6 elements in, 6 elements out)
1458 : !!
1459 : !! rank=3: bb(i,j,k)= $[{5 \over 2} g(i,l) g(j,m) g(k,n) - {3 \over 2} g(i,j) g(l,m) g(k,n)] aa(l,m,n)$
1460 : !! (10 elements in, 10 elements out)
1461 : !! In this rank 3 case, the second term is NOT symmetric in all
1462 : !! permutations of i,j,k, but the final tensor b(ijk) may be
1463 : !! symmetrized over all permutations because it will be
1464 : !! contracted with a completely symmetric tensor.
1465 : !!
1466 : !! The compressed storage scheme is based on storing a symmetric 3x3 matrix as
1467 : !! (1 . .)
1468 : !! (6 2 .)
1469 : !! (5 4 3)
1470 : !! which leads to the following mappings for all ranks
1471 : !! where the compressed storage index is to the right of the arrow:
1472 : !! rank=0 1->1 (only a scalar)
1473 : !! rank=1 1->1 2->2 3->3 (standard vector, no compression)
1474 : !! rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1475 : !! weights 1 1 1 2 2 2
1476 : !! rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
1477 : !! weights 1 3 3 6 3 3 1 3 3 1
1478 : !!
1479 : !! SOURCE
1480 :
1481 292310138 : subroutine metcon(rank,gmet,aa,bb)
1482 :
1483 : !Arguments ------------------------------------
1484 : !scalars
1485 : integer,intent(in) :: rank
1486 : !arrays
1487 : real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),gmet(3,3)
1488 : real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2)
1489 :
1490 : !Local variables-------------------------------
1491 : !scalars
1492 : integer :: ii,jj
1493 : real(dp) :: scalar,tmpiii,tmpijk
1494 : character(len=500) :: message
1495 : !arrays
1496 : real(dp) :: vector(3)
1497 :
1498 : ! *************************************************************************
1499 :
1500 : !This statement function defines the l=3 contraction in
1501 : !terms of the free indices of the contracted tensor (re and im)
1502 : ! coniii(ii,i1,i2,i3)=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+&
1503 : !& gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+&
1504 : !& gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10)
1505 : ! conijk(ii,i1,i2,i3)=aa(ii,4)*&
1506 : !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+&
1507 : !& gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+&
1508 : !& gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+&
1509 : !& gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+&
1510 : !& gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+&
1511 : !& gmet(i1,2)*gmet(i2,1)*gmet(i3,3))
1512 : ! con(ii,i1,i2,i3)=coniii(ii,i1,i2,i3)+conijk(ii,i1,i2,i3)+&
1513 : !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+&
1514 : !& gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+&
1515 : !& gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+&
1516 : !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+&
1517 : !& gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+&
1518 : !& gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+&
1519 : !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+&
1520 : !& gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+&
1521 : !& gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+&
1522 : !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+&
1523 : !& gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+&
1524 : !& gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+&
1525 : !& (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+&
1526 : !& gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+&
1527 : !& gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+&
1528 : !& (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+&
1529 : !& gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+&
1530 : !& gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8)
1531 :
1532 : !DEBUG
1533 : !write(std_out,*)' metcon : enter '
1534 : !stop
1535 : !ENDDEBUG
1536 292310138 : if (rank==0) then
1537 : ! Simply copy scalar, re and im
1538 115183099 : bb(1,1)=aa(1,1)
1539 115183099 : bb(2,1)=aa(2,1)
1540 :
1541 177127039 : else if (rank==1) then
1542 : ! Apply gmet to input vector, re and im
1543 325552500 : do ii=1,2
1544 976657500 : do jj=1,3
1545 868140000 : bb(ii,jj)=gmet(jj,1)*aa(ii,1)+gmet(jj,2)*aa(ii,2)+gmet(jj,3)*aa(ii,3)
1546 : end do
1547 : end do
1548 :
1549 :
1550 68609539 : else if (rank==2) then
1551 : ! Apply rank2 expression, re and im
1552 203505462 : do ii=1,2
1553 : ! Carry out g(c,d)*aa(c,d) contraction to get scalar
1554 : scalar=gmet(1,1)*aa(ii,1)+gmet(2,2)*aa(ii,2)+&
1555 : & gmet(3,3)*aa(ii,3)+2.0d0*(gmet(2,1)*aa(ii,6)+&
1556 135670308 : & gmet(3,1)*aa(ii,5)+gmet(3,2)*aa(ii,4) )
1557 : ! Write out components of contraction
1558 : ! (1,1)->1
1559 : bb(ii,1)=0.5d0*(3.0d0*(gmet(1,1)*gmet(1,1)*aa(ii,1)+&
1560 : & gmet(1,2)*gmet(1,2)*aa(ii,2)+gmet(1,3)*gmet(1,3)*aa(ii,3)+&
1561 : & (gmet(1,2)*gmet(1,1)+gmet(1,1)*gmet(1,2))*aa(ii,6)+&
1562 : & (gmet(1,3)*gmet(1,1)+gmet(1,1)*gmet(1,3))*aa(ii,5)+&
1563 : & (gmet(1,3)*gmet(1,2)+gmet(1,2)*gmet(1,3))*aa(ii,4) ) &
1564 135670308 : & - gmet(1,1)*scalar)
1565 : ! (2,2)->2
1566 : bb(ii,2)=0.5d0*(3.0d0*(gmet(2,1)*gmet(2,1)*aa(ii,1)+&
1567 : & gmet(2,2)*gmet(2,2)*aa(ii,2)+gmet(2,3)*gmet(2,3)*aa(ii,3)+&
1568 : & (gmet(2,2)*gmet(2,1)+gmet(2,1)*gmet(2,2))*aa(ii,6)+&
1569 : & (gmet(2,3)*gmet(2,1)+gmet(2,1)*gmet(2,3))*aa(ii,5)+&
1570 : & (gmet(2,3)*gmet(2,2)+gmet(2,2)*gmet(2,3))*aa(ii,4) )&
1571 135670308 : & - gmet(2,2)*scalar)
1572 : ! (3,3)->3
1573 : bb(ii,3)=0.5d0*(3.0d0*(gmet(3,1)*gmet(3,1)*aa(ii,1)+&
1574 : & gmet(3,2)*gmet(3,2)*aa(ii,2)+gmet(3,3)*gmet(3,3)*aa(ii,3)+&
1575 : & (gmet(3,2)*gmet(3,1)+gmet(3,1)*gmet(3,2))*aa(ii,6)+&
1576 : & (gmet(3,3)*gmet(3,1)+gmet(3,1)*gmet(3,3))*aa(ii,5)+&
1577 : & (gmet(3,3)*gmet(3,2)+gmet(3,2)*gmet(3,3))*aa(ii,4) )&
1578 135670308 : & - gmet(3,3)*scalar)
1579 : ! (3,2)->4
1580 : bb(ii,4)=0.5d0*(3.0d0*(gmet(3,1)*gmet(2,1)*aa(ii,1)+&
1581 : & gmet(3,2)*gmet(2,2)*aa(ii,2)+gmet(3,3)*gmet(2,3)*aa(ii,3)+&
1582 : & (gmet(3,2)*gmet(2,1)+gmet(3,1)*gmet(2,2))*aa(ii,6)+&
1583 : & (gmet(3,3)*gmet(2,1)+gmet(3,1)*gmet(2,3))*aa(ii,5)+&
1584 : & (gmet(3,3)*gmet(2,2)+gmet(3,2)*gmet(2,3))*aa(ii,4) )&
1585 135670308 : & - gmet(3,2)*scalar)
1586 : ! (3,1)->5
1587 : bb(ii,5)=0.5d0*(3.0d0*(gmet(3,1)*gmet(1,1)*aa(ii,1)+&
1588 : & gmet(3,2)*gmet(1,2)*aa(ii,2)+gmet(3,3)*gmet(1,3)*aa(ii,3)+&
1589 : & (gmet(3,2)*gmet(1,1)+gmet(3,1)*gmet(1,2))*aa(ii,6)+&
1590 : & (gmet(3,3)*gmet(1,1)+gmet(3,1)*gmet(1,3))*aa(ii,5)+&
1591 : & (gmet(3,3)*gmet(1,2)+gmet(3,2)*gmet(1,3))*aa(ii,4) )&
1592 135670308 : & - gmet(3,1)*scalar)
1593 : ! (2,1)->6
1594 : bb(ii,6)=0.5d0*(3.0d0*(gmet(2,1)*gmet(1,1)*aa(ii,1)+&
1595 : & gmet(2,2)*gmet(1,2)*aa(ii,2)+gmet(2,3)*gmet(1,3)*aa(ii,3)+&
1596 : & (gmet(2,2)*gmet(1,1)+gmet(2,1)*gmet(1,2))*aa(ii,6)+&
1597 : & (gmet(2,3)*gmet(1,1)+gmet(2,1)*gmet(1,3))*aa(ii,5)+&
1598 : & (gmet(2,3)*gmet(1,2)+gmet(2,2)*gmet(1,3))*aa(ii,4) ) &
1599 135670308 : & - gmet(2,1)*scalar)
1600 : ! Include appropriate weights for multiplicity
1601 135670308 : bb(ii,4)=2.d0*bb(ii,4)
1602 135670308 : bb(ii,5)=2.d0*bb(ii,5)
1603 203505462 : bb(ii,6)=2.d0*bb(ii,6)
1604 : end do
1605 :
1606 774385 : else if (rank==3) then
1607 : ! Apply rank2 expression, re and im
1608 2323155 : do ii=1,2
1609 : ! Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j)
1610 6195080 : do jj=1,3
1611 : tmpiii= gmet(1,1)*gmet(jj,1)*aa(ii,1)+&
1612 : & gmet(2,2)*gmet(jj,2)*aa(ii,7)+&
1613 4646310 : & gmet(3,3)*gmet(jj,3)*aa(ii,10)
1614 : tmpijk= (gmet(1,2)*gmet(jj,3)+&
1615 : & gmet(3,1)*gmet(jj,2)+&
1616 : & gmet(2,3)*gmet(jj,1)+&
1617 : & gmet(3,2)*gmet(jj,1)+&
1618 : & gmet(1,3)*gmet(jj,2)+&
1619 4646310 : & gmet(2,1)*gmet(jj,3)) *aa(ii,4)
1620 : vector(jj)=tmpiii + tmpijk +&
1621 : & (gmet(1,2)*gmet(jj,1)+&
1622 : & gmet(2,1)*gmet(jj,1)+&
1623 : & gmet(1,1)*gmet(jj,2)) *aa(ii,6)+&
1624 : & (gmet(1,2)*gmet(jj,2)+&
1625 : & gmet(2,1)*gmet(jj,2)+&
1626 : & gmet(2,2)*gmet(jj,1)) *aa(ii,2)+&
1627 : & (gmet(1,3)*gmet(jj,1)+&
1628 : & gmet(3,1)*gmet(jj,1)+&
1629 : & gmet(1,1)*gmet(jj,3)) *aa(ii,5)+&
1630 : & (gmet(1,3)*gmet(jj,3)+&
1631 : & gmet(3,1)*gmet(jj,3)+&
1632 : & gmet(3,3)*gmet(jj,1)) *aa(ii,3)+&
1633 : & (gmet(2,3)*gmet(jj,2)+&
1634 : & gmet(3,2)*gmet(jj,2)+&
1635 : & gmet(2,2)*gmet(jj,3)) *aa(ii,9)+&
1636 : & (gmet(2,3)*gmet(jj,3)+&
1637 : & gmet(3,2)*gmet(jj,3)+&
1638 6195080 : & gmet(3,3)*gmet(jj,2)) *aa(ii,8)
1639 : end do
1640 : ! Write out components of contraction
1641 : ! (111)->1
1642 1548770 : bb(ii,1) =2.5d0*con_met(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1))
1643 : ! (221)->2
1644 : bb(ii,2) =2.5d0*con_met(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+&
1645 1548770 : & gmet(1,2)*vector(2)+gmet(2,2)*vector(1))
1646 : ! (331)->3
1647 : bb(ii,3) =2.5d0*con_met(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+&
1648 1548770 : & gmet(1,3)*vector(3)+gmet(3,3)*vector(1))
1649 : ! (321)->4
1650 : bb(ii,4) =2.5d0*con_met(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+&
1651 1548770 : & gmet(1,2)*vector(3)+gmet(3,2)*vector(1))
1652 : ! (311)->5
1653 : bb(ii,5) =2.5d0*con_met(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+&
1654 1548770 : & gmet(1,1)*vector(3)+gmet(3,1)*vector(1))
1655 : ! (211)->6
1656 : bb(ii,6) =2.5d0*con_met(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+&
1657 1548770 : & gmet(1,1)*vector(2)+gmet(2,1)*vector(1))
1658 : ! (222)->7
1659 1548770 : bb(ii,7) =2.5d0*con_met(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2))
1660 :
1661 : ! (332)->8
1662 : bb(ii,8) =2.5d0*con_met(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+&
1663 1548770 : & gmet(2,3)*vector(3)+gmet(3,3)*vector(2))
1664 : ! (322)->9
1665 : bb(ii,9) =2.5d0*con_met(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+&
1666 1548770 : & gmet(2,2)*vector(3)+gmet(3,2)*vector(2))
1667 : ! (333)->10
1668 1548770 : bb(ii,10)=2.5d0*con_met(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3))
1669 : ! Include appropriate weights for multiplicity
1670 1548770 : bb(ii,2)=3.d0*bb(ii,2)
1671 1548770 : bb(ii,3)=3.d0*bb(ii,3)
1672 1548770 : bb(ii,4)=6.d0*bb(ii,4)
1673 1548770 : bb(ii,5)=3.d0*bb(ii,5)
1674 1548770 : bb(ii,6)=3.d0*bb(ii,6)
1675 1548770 : bb(ii,8)=3.d0*bb(ii,8)
1676 2323155 : bb(ii,9)=3.d0*bb(ii,9)
1677 : end do
1678 :
1679 : else
1680 : write(message, '(a,i0,a,a,a)' )&
1681 0 : & 'Input rank=',rank,' not allowed.',ch10,&
1682 0 : & 'Possible values are 0,1,2,3 only.'
1683 0 : ABI_BUG(message)
1684 : end if
1685 :
1686 : contains
1687 :
1688 15487700 : function con_met(ii,i1,i2,i3)
1689 :
1690 : real(dp) :: con_met
1691 : integer :: ii,i1,i2,i3
1692 : real(dp)::coniii,conijk
1693 :
1694 : coniii=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+&
1695 : & gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+&
1696 15487700 : & gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10)
1697 : conijk=aa(ii,4)*&
1698 : & (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+&
1699 : & gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+&
1700 : & gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+&
1701 : & gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+&
1702 : & gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+&
1703 15487700 : & gmet(i1,2)*gmet(i2,1)*gmet(i3,3))
1704 : con_met=coniii+conijk+&
1705 : & (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+&
1706 : & gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+&
1707 : & gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+&
1708 : & (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+&
1709 : & gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+&
1710 : & gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+&
1711 : & (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+&
1712 : & gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+&
1713 : & gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+&
1714 : & (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+&
1715 : & gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+&
1716 : & gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+&
1717 : & (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+&
1718 : & gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+&
1719 : & gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+&
1720 : & (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+&
1721 : & gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+&
1722 15487700 : & gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8)
1723 :
1724 15487700 : end function con_met
1725 :
1726 : end subroutine metcon
1727 : !!***
1728 :
1729 : !!****f* m_contract/metcon_so
1730 : !! NAME
1731 : !! metcon_so
1732 : !!
1733 : !! FUNCTION
1734 : !! Carries out specialized metric tensor contractions needed for
1735 : !! l=0,1,2,3 nonlocal Kleinman-Bylander pseudopotential operation
1736 : !! in the spin-orbit case.
1737 : !! Full advantage is taken of the full permutational symmetry of these tensors.
1738 : !!
1739 : !! INPUTS
1740 : !! rank=0,1,2, or 3 = rank of input tensor aa
1741 : !! gmet(3,3)=metric tensor (array is symmetric but stored as 3x3)
1742 : !! amet(3,3)=real or imaginary part of one spin matrix element of the
1743 : !! "spin metric" tensor
1744 : !! aa(2,(rank+1)*(rank+2)/2)=unique elements of complex input tensor
1745 : !!
1746 : !! OUTPUT
1747 : !! bb(2,(rank+1)*(rank+2)/2)=unique elements of complex output tensor
1748 : !!
1749 : !! NOTES
1750 : !! All tensors are stored in a compressed storage mode defined below;
1751 : !! input and output conform to this scheme.
1752 : !! When tensor elements occur repeatedly due to symmetry, the
1753 : !! WEIGHT IS INCLUDED in the output tensor element to simplify later
1754 : !! contractions with other tensors of the same rank and form, i.e. the
1755 : !! next contraction is then simply a dot product over the unique elements.
1756 : !!
1757 : !! Definitions of the contractions:
1758 : !!
1759 : !! rank=0: bb=0
1760 : !!
1761 : !! rank=1: bb(i)= $amet(i,l) aa(l)$ (3 elements in, 3 elements out)
1762 : !!
1763 : !! rank=2: bb(i,j)= $[3 gmet(i,l) amet(j,m)] aa(l,m)$
1764 : !! (6 elements in, 6 elements out)
1765 : !!
1766 : !! rank=3: bb(i,j,k)= $[{15 \over 2} g(i,l) g(j,m) a(k,n) - {3 \over 2} g(i,j) g(l,m) a(k,n)] aa(l,m,n)$
1767 : !! (10 elements in, 10 elements out)
1768 : !! In this rank 3 case, the second term is NOT symmetric in all
1769 : !! permutations of i,j,k, but the final tensor b(ijk) may be
1770 : !! symmetrized over all permutations because it will be
1771 : !! contracted with a completely symmetric tensor.
1772 : !!
1773 : !! The compressed storage scheme is based on storing
1774 : !! a symmetric 3x3 matrix as
1775 : !! (1 . .)
1776 : !! (6 2 .)
1777 : !! (5 4 3)
1778 : !! which leads to the following mappings for all ranks
1779 : !! where the compressed storage index is to the right of the arrow:
1780 : !! rank=0 1->1 (only a scalar)
1781 : !! rank=1 1->1 2->2 3->3 (standard vector, no compression)
1782 : !! rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1783 : !! weights 1 1 1 2 2 2
1784 : !! rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
1785 : !! weights 1 3 3 6 3 3 1 3 3 1
1786 : !!
1787 : !! SOURCE
1788 :
1789 276604864 : subroutine metcon_so(rank,gmet,amet,aa,bb)
1790 :
1791 : !Arguments ------------------------------------
1792 : !scalars
1793 : integer,intent(in) :: rank
1794 : !arrays
1795 : real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),amet(3,3),gmet(3,3)
1796 : real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2)
1797 :
1798 : !Local variables-------------------------------
1799 : !scalars
1800 : integer :: ii,jj
1801 : real(dp) :: tmpiii,tmpijk
1802 : character(len=500) :: message
1803 : !arrays
1804 : real(dp) :: vector(3)
1805 :
1806 : ! *************************************************************************
1807 :
1808 276604864 : if (rank==0) then
1809 : ! Simply copy scalar, re and im
1810 0 : bb(1,1)=0.d0
1811 0 : bb(2,1)=0.d0
1812 : !
1813 276604864 : else if (rank==1) then
1814 : ! Apply gmet to input vector, re and im
1815 449117088 : do ii=1,2
1816 1347351264 : do jj=1,3
1817 1197645568 : bb(ii,jj)=amet(jj,1)*aa(ii,1)+amet(jj,2)*aa(ii,2)+amet(jj,3)*aa(ii,3)
1818 : end do
1819 : end do
1820 :
1821 126899168 : else if (rank==2) then
1822 : ! Apply rank2 expression, re and im
1823 377305056 : do ii=1,2
1824 : ! Write out components of contraction
1825 : ! (1,1)->1
1826 : bb(ii,1)=3.0d0*(gmet(1,1)*amet(1,1)*aa(ii,1)+&
1827 : & gmet(1,2)*amet(1,2)*aa(ii,2)+gmet(1,3)*amet(1,3)*aa(ii,3)+&
1828 : & (gmet(1,2)*amet(1,1)+gmet(1,1)*amet(1,2))*aa(ii,6)+&
1829 : & (gmet(1,3)*amet(1,1)+gmet(1,1)*amet(1,3))*aa(ii,5)+&
1830 251536704 : & (gmet(1,3)*amet(1,2)+gmet(1,2)*amet(1,3))*aa(ii,4))
1831 : ! (2,2)->2
1832 : bb(ii,2)=3.0d0*(gmet(2,1)*amet(2,1)*aa(ii,1)+&
1833 : & gmet(2,2)*amet(2,2)*aa(ii,2)+gmet(2,3)*amet(2,3)*aa(ii,3)+&
1834 : & (gmet(2,2)*amet(2,1)+gmet(2,1)*amet(2,2))*aa(ii,6)+&
1835 : & (gmet(2,3)*amet(2,1)+gmet(2,1)*amet(2,3))*aa(ii,5)+&
1836 251536704 : & (gmet(2,3)*amet(2,2)+gmet(2,2)*amet(2,3))*aa(ii,4) )
1837 : ! (3,3)->3
1838 : bb(ii,3)=3.0d0*(gmet(3,1)*amet(3,1)*aa(ii,1)+&
1839 : & gmet(3,2)*amet(3,2)*aa(ii,2)+gmet(3,3)*amet(3,3)*aa(ii,3)+&
1840 : & (gmet(3,2)*amet(3,1)+gmet(3,1)*amet(3,2))*aa(ii,6)+&
1841 : & (gmet(3,3)*amet(3,1)+gmet(3,1)*amet(3,3))*aa(ii,5)+&
1842 251536704 : & (gmet(3,3)*amet(3,2)+gmet(3,2)*amet(3,3))*aa(ii,4) )
1843 : ! (3,2)->4
1844 : bb(ii,4)=3.0d0*(gmet(3,1)*amet(2,1)*aa(ii,1)+&
1845 : & gmet(3,2)*amet(2,2)*aa(ii,2)+gmet(3,3)*amet(2,3)*aa(ii,3)+&
1846 : & (gmet(3,2)*amet(2,1)+gmet(3,1)*amet(2,2))*aa(ii,6)+&
1847 : & (gmet(3,3)*amet(2,1)+gmet(3,1)*amet(2,3))*aa(ii,5)+&
1848 251536704 : & (gmet(3,3)*amet(2,2)+gmet(3,2)*amet(2,3))*aa(ii,4) )
1849 : bb(ii,4)=bb(ii,4)+3.0d0*(amet(3,1)*gmet(2,1)*aa(ii,1)+&
1850 : & amet(3,2)*gmet(2,2)*aa(ii,2)+amet(3,3)*gmet(2,3)*aa(ii,3)+&
1851 : & (amet(3,2)*gmet(2,1)+amet(3,1)*gmet(2,2))*aa(ii,6)+&
1852 : & (amet(3,3)*gmet(2,1)+amet(3,1)*gmet(2,3))*aa(ii,5)+&
1853 251536704 : & (amet(3,3)*gmet(2,2)+amet(3,2)*gmet(2,3))*aa(ii,4) )
1854 : ! (3,1)->5
1855 : bb(ii,5)=3.0d0*(gmet(3,1)*amet(1,1)*aa(ii,1)+&
1856 : & gmet(3,2)*amet(1,2)*aa(ii,2)+gmet(3,3)*amet(1,3)*aa(ii,3)+&
1857 : & (gmet(3,2)*amet(1,1)+gmet(3,1)*amet(1,2))*aa(ii,6)+&
1858 : & (gmet(3,3)*amet(1,1)+gmet(3,1)*amet(1,3))*aa(ii,5)+&
1859 251536704 : & (gmet(3,3)*amet(1,2)+gmet(3,2)*amet(1,3))*aa(ii,4) )
1860 : bb(ii,5)=bb(ii,5)+3.0d0*(amet(3,1)*gmet(1,1)*aa(ii,1)+&
1861 : & amet(3,2)*gmet(1,2)*aa(ii,2)+amet(3,3)*gmet(1,3)*aa(ii,3)+&
1862 : & (amet(3,2)*gmet(1,1)+amet(3,1)*gmet(1,2))*aa(ii,6)+&
1863 : & (amet(3,3)*gmet(1,1)+amet(3,1)*gmet(1,3))*aa(ii,5)+&
1864 251536704 : & (amet(3,3)*gmet(1,2)+amet(3,2)*gmet(1,3))*aa(ii,4) )
1865 : ! (2,1)->6
1866 : bb(ii,6)=3.0d0*(gmet(2,1)*amet(1,1)*aa(ii,1)+&
1867 : & gmet(2,2)*amet(1,2)*aa(ii,2)+gmet(2,3)*amet(1,3)*aa(ii,3)+&
1868 : & (gmet(2,2)*amet(1,1)+gmet(2,1)*amet(1,2))*aa(ii,6)+&
1869 : & (gmet(2,3)*amet(1,1)+gmet(2,1)*amet(1,3))*aa(ii,5)+&
1870 251536704 : & (gmet(2,3)*amet(1,2)+gmet(2,2)*amet(1,3))*aa(ii,4) )
1871 : bb(ii,6)=bb(ii,6)+3.0d0*(amet(2,1)*gmet(1,1)*aa(ii,1)+&
1872 : & amet(2,2)*gmet(1,2)*aa(ii,2)+amet(2,3)*gmet(1,3)*aa(ii,3)+&
1873 : & (amet(2,2)*gmet(1,1)+amet(2,1)*gmet(1,2))*aa(ii,6)+&
1874 : & (amet(2,3)*gmet(1,1)+amet(2,1)*gmet(1,3))*aa(ii,5)+&
1875 251536704 : & (amet(2,3)*gmet(1,2)+amet(2,2)*gmet(1,3))*aa(ii,4) )
1876 : ! Include appropriate weights for multiplicity
1877 : bb(ii,4)=bb(ii,4)
1878 : bb(ii,5)=bb(ii,5)
1879 125768352 : bb(ii,6)=bb(ii,6)
1880 : end do
1881 :
1882 1130816 : else if (rank==3) then
1883 : ! Apply rank2 expression, re and im
1884 3392448 : do ii=1,2
1885 : ! Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j)
1886 9046528 : do jj=1,3
1887 : tmpiii= gmet(1,1)*amet(jj,1)*aa(ii,1)+&
1888 : & gmet(2,2)*amet(jj,2)*aa(ii,7)+&
1889 6784896 : & gmet(3,3)*amet(jj,3)*aa(ii,10)
1890 : tmpijk= (gmet(1,2)*amet(jj,3)+&
1891 : & gmet(3,1)*amet(jj,2)+&
1892 : & gmet(2,3)*amet(jj,1)+&
1893 : & gmet(3,2)*amet(jj,1)+&
1894 : & gmet(1,3)*amet(jj,2)+&
1895 6784896 : & gmet(2,1)*amet(jj,3)) *aa(ii,4)
1896 : vector(jj)=tmpiii + tmpijk +&
1897 : & (gmet(1,2)*amet(jj,1)+&
1898 : & gmet(2,1)*amet(jj,1)+&
1899 : & gmet(1,1)*amet(jj,2)) *aa(ii,6)+&
1900 : & (gmet(1,2)*amet(jj,2)+&
1901 : & gmet(2,1)*amet(jj,2)+&
1902 : & gmet(2,2)*amet(jj,1)) *aa(ii,2)+&
1903 : & (gmet(1,3)*amet(jj,1)+&
1904 : & gmet(3,1)*amet(jj,1)+&
1905 : & gmet(1,1)*amet(jj,3)) *aa(ii,5)+&
1906 : & (gmet(1,3)*amet(jj,3)+&
1907 : & gmet(3,1)*amet(jj,3)+&
1908 : & gmet(3,3)*amet(jj,1)) *aa(ii,3)+&
1909 : & (gmet(2,3)*amet(jj,2)+&
1910 : & gmet(3,2)*amet(jj,2)+&
1911 : & gmet(2,2)*amet(jj,3)) *aa(ii,9)+&
1912 : & (gmet(2,3)*amet(jj,3)+&
1913 : & gmet(3,2)*amet(jj,3)+&
1914 9046528 : & gmet(3,3)*amet(jj,2)) *aa(ii,8)
1915 : end do
1916 : ! Write out components of contraction
1917 : ! (111)->1
1918 2261632 : bb(ii,1) =7.5d0*con_metso(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1))
1919 : ! (221)->2
1920 : bb(ii,2) =7.5d0*con_metso(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+&
1921 2261632 : & gmet(1,2)*vector(2)+gmet(2,2)*vector(1))
1922 : ! (331)->3
1923 : bb(ii,3) =7.5d0*con_metso(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+&
1924 2261632 : & gmet(1,3)*vector(3)+gmet(3,3)*vector(1))
1925 : ! (321)->4
1926 : bb(ii,4) =7.5d0*con_metso(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+&
1927 2261632 : & gmet(1,2)*vector(3)+gmet(3,2)*vector(1))
1928 : ! (311)->5
1929 : bb(ii,5) =7.5d0*con_metso(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+&
1930 2261632 : & gmet(1,1)*vector(3)+gmet(3,1)*vector(1))
1931 : ! (211)->6
1932 : bb(ii,6) =7.5d0*con_metso(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+&
1933 2261632 : & gmet(1,1)*vector(2)+gmet(2,1)*vector(1))
1934 : ! (222)->7
1935 2261632 : bb(ii,7) =7.5d0*con_metso(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2))
1936 :
1937 : ! (332)->8
1938 : bb(ii,8) =7.5d0*con_metso(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+&
1939 2261632 : & gmet(2,3)*vector(3)+gmet(3,3)*vector(2))
1940 : ! (322)->9
1941 : bb(ii,9) =7.5d0*con_metso(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+&
1942 2261632 : & gmet(2,2)*vector(3)+gmet(3,2)*vector(2))
1943 : ! (333)->10
1944 2261632 : bb(ii,10)=7.5d0*con_metso(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3))
1945 : ! Include appropriate weights for multiplicity
1946 2261632 : bb(ii,2)=3.d0*bb(ii,2)
1947 2261632 : bb(ii,3)=3.d0*bb(ii,3)
1948 2261632 : bb(ii,4)=6.d0*bb(ii,4)
1949 2261632 : bb(ii,5)=3.d0*bb(ii,5)
1950 2261632 : bb(ii,6)=3.d0*bb(ii,6)
1951 2261632 : bb(ii,8)=3.d0*bb(ii,8)
1952 3392448 : bb(ii,9)=3.d0*bb(ii,9)
1953 : end do
1954 :
1955 : else
1956 : write(message, '(a,i0,a,a,a)' )&
1957 0 : & 'Input rank=',rank,' not allowed.',ch10,&
1958 0 : & 'Possible values are 0,1,2,3 only.'
1959 0 : ABI_BUG(message)
1960 : end if
1961 :
1962 : contains
1963 :
1964 : !This function defines the l=3 contraction in
1965 : !terms of the free indices of the contracted tensor (re and im)
1966 :
1967 67848960 : function cona_metso(ii,i1,i2,i3)
1968 :
1969 : real(dp) :: cona_metso
1970 : integer,intent(in) :: ii,i1,i2,i3
1971 : real(dp) :: coniii, conijk
1972 :
1973 : coniii=gmet(i1,1)*gmet(i2,1)*amet(i3,1)*aa(ii,1)+&
1974 : & gmet(i1,2)*gmet(i2,2)*amet(i3,2)*aa(ii,7)+&
1975 67848960 : & gmet(i1,3)*gmet(i2,3)*amet(i3,3)*aa(ii,10)
1976 : conijk=aa(ii,4)*&
1977 : & (gmet(i1,1)*gmet(i2,2)*amet(i3,3)+&
1978 : & gmet(i1,2)*gmet(i2,3)*amet(i3,1)+&
1979 : & gmet(i1,3)*gmet(i2,1)*amet(i3,2)+&
1980 : & gmet(i1,3)*gmet(i2,2)*amet(i3,1)+&
1981 : & gmet(i1,1)*gmet(i2,3)*amet(i3,2)+&
1982 67848960 : & gmet(i1,2)*gmet(i2,1)*amet(i3,3))
1983 : cona_metso=coniii+conijk+&
1984 : & (gmet(i1,1)*gmet(i2,2)*amet(i3,1)+&
1985 : & gmet(i1,2)*gmet(i2,1)*amet(i3,1)+&
1986 : & gmet(i1,1)*gmet(i2,1)*amet(i3,2))*aa(ii,6)+&
1987 : & (gmet(i1,1)*gmet(i2,2)*amet(i3,2)+&
1988 : & gmet(i1,2)*gmet(i2,1)*amet(i3,2)+&
1989 : & gmet(i1,2)*gmet(i2,2)*amet(i3,1))*aa(ii,2)+&
1990 : & (gmet(i1,1)*gmet(i2,3)*amet(i3,1)+&
1991 : & gmet(i1,3)*gmet(i2,1)*amet(i3,1)+&
1992 : & gmet(i1,1)*gmet(i2,1)*amet(i3,3))*aa(ii,5)+&
1993 : & (gmet(i1,1)*gmet(i2,3)*amet(i3,3)+&
1994 : & gmet(i1,3)*gmet(i2,1)*amet(i3,3)+&
1995 : & gmet(i1,3)*gmet(i2,3)*amet(i3,1))*aa(ii,3)+&
1996 : & (gmet(i1,2)*gmet(i2,2)*amet(i3,3)+&
1997 : & gmet(i1,2)*gmet(i2,3)*amet(i3,2)+&
1998 : & gmet(i1,3)*gmet(i2,2)*amet(i3,2))*aa(ii,9)+&
1999 : & (gmet(i1,2)*gmet(i2,3)*amet(i3,3)+&
2000 : & gmet(i1,3)*gmet(i2,2)*amet(i3,3)+&
2001 67848960 : & gmet(i1,3)*gmet(i2,3)*amet(i3,2))*aa(ii,8)
2002 67848960 : end function cona_metso
2003 :
2004 :
2005 22616320 : function con_metso(ii,i1,i2,i3)
2006 :
2007 : real(dp) :: con_metso
2008 : integer,intent(in) :: ii,i1,i2,i3
2009 :
2010 22616320 : con_metso=(cona_metso(ii,i3,i1,i2)+cona_metso(ii,i2,i3,i1)+cona_metso(ii,i1,i2,i3))/3.d0
2011 :
2012 22616320 : end function con_metso
2013 :
2014 : end subroutine metcon_so
2015 : !!***
2016 :
2017 : !!****f* m_contract/metric_so
2018 : !! NAME
2019 : !! metric_so
2020 : !!
2021 : !! FUNCTION
2022 : !! Computes Pauli matrices and antisymmetric tensor needed for
2023 : !! spin-orbit.
2024 : !!
2025 : !! INPUTS
2026 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
2027 : !! (bohr**-1)
2028 : !! soc_weight(3)=prefactors for the spin-orbit components (dimensionless):
2029 : !! (soc_weight(1),soc_weight(2),soc_weight(3))
2030 : !! scale the \sigma_x, \sigma_y, \sigma_z contributions, respectively.
2031 : !! alpha=Euler angle for rotation around z-axis
2032 : !! beta =Euler angle for rotation around y-axis
2033 : !!
2034 : !! OUTPUT
2035 : !! amet(2,3,3,2,2)=the antisymmetric tensor A(Re/Im,y,y'',s,s'')
2036 : !! pauli(2,2,2,3)=Pauli matrixes
2037 : !!
2038 : !! NOTES
2039 : !! (the preprocessing based on CPP does not allow isolated quotes,
2040 : !! so that they have been doubled in the following latex equations)
2041 : !!{{\ \begin{eqnarray}
2042 : !! $2 Pauli(s,s'',1) & = & 0 & 1 \nonumber
2043 : !! & & 1 & 0 \nonumber
2044 : !! 2 Pauli(s,s'',2) & = & 0 & -i \nonumber
2045 : !! & & i & 0 \nonumber
2046 : !! 2 Pauli(s,s'',3) & = & 1 & 0 \nonumber
2047 : !! & & 0 & -1
2048 : !! \end{eqnarray} }}
2049 : !!{{\ \begin{eqnarray}
2050 : !! $Amet(y,y'',s,s'') & = & -i Pauli(s,s'',n) E(n,m,m'') Gprimd(m,y) Gprimd(m'',y'')
2051 : !! \end{eqnarray} }}
2052 : !!
2053 : !! E(n,m,m''): full antisymetric tensor
2054 : !! s,s'': spin indices (1..2)
2055 : !! y,y'',m,m'',n: metric indices (1..3)
2056 : !! a,b: strain indices (1..3)
2057 : !! Amet and Pauli are complex
2058 : !!
2059 : !! SOURCE
2060 :
2061 2339994 : subroutine metric_so(amet,soc_weight,gprimd,pauli,alpha,beta)
2062 :
2063 : !Arguments ------------------------------------
2064 : !scalars
2065 : real(dp),intent(in) :: alpha, beta
2066 : !arrays
2067 : real(dp),intent(in) :: gprimd(3,3),soc_weight(3)
2068 : real(dp),intent(out) :: amet(2,3,3,2,2),pauli(2,2,2,3)
2069 :
2070 : !Local variables-------------------------------
2071 : !scalars
2072 : integer :: iy1,iy2,m1,m2,n
2073 : !arrays
2074 : real(dp) :: buffer1(3,3,2,2), cb2, sb2 !,buffer2(3,3,3,3,2,2)
2075 : complex(dp) :: S(2,2,3), Srot(2,2), U(2,2), ep, em
2076 : ! **********************************************************************
2077 :
2078 : !Fill in Pauli matrices and make them spin matrices:
2079 : !Pauli(Re/Im,up,down,coord):
2080 2339994 : pauli(:,:,:,:)=0.d0
2081 2339994 : pauli(1,1,2,1)= 1.d0;pauli(1,2,1,1)= 1.d0
2082 2339994 : pauli(2,1,2,2)=-1.d0;pauli(2,2,1,2)= 1.d0
2083 2339994 : pauli(1,1,1,3)= 1.d0;pauli(1,2,2,3)=-1.d0
2084 107639724 : pauli(:,:,:,:)= 0.5d0*pauli(:,:,:,:)
2085 :
2086 2339994 : if (.not.(abs(alpha) < tol8 .and. abs(beta) < tol8)) then
2087 280258 : S(:,:,:) = cmplx(pauli(1,:,:,:), pauli(2,:,:,:), kind=dp)
2088 :
2089 12739 : cb2 = cos(half*beta); sb2 = sin(half*beta)
2090 12739 : em = exp(-j_dpc*half*alpha); ep = conjg(em)
2091 :
2092 12739 : U(1,1) = cb2 * em; U(1,2) = -sb2 * em
2093 12739 : U(2,1) = sb2 * ep; U(2,2) = cb2 * ep
2094 :
2095 50956 : do n = 1, 3
2096 1108293 : Srot(:,:) = matmul(conjg(transpose(U)), matmul(S(:,:,n), U))
2097 280258 : S(:,:,n) = Srot(:,:)
2098 : end do
2099 :
2100 12739 : pauli(:,:,:,:) = zero
2101 280258 : pauli(1,:,:,:) = real(S(:,:,:), kind=dp)
2102 2607513 : pauli(2,:,:,:) = aimag(S(:,:,:))
2103 : end if
2104 :
2105 : !Construct the antisymmetric tensor:
2106 2339994 : amet(:,:,:,:,:)=0.d0
2107 9359976 : do iy2=1,3
2108 30419922 : do iy1=1,3
2109 91259766 : do n=1,3
2110 63179838 : m1=mod(n ,3)+1 ! n,m1,m2 is an even permutation
2111 63179838 : m2=mod(m1,3)+1
2112 : amet(1:2,iy1,iy2,1:2,1:2) = amet(:,iy1,iy2,:,:) &
2113 : & + soc_weight(n)*pauli(:,:,:,n) &
2114 : & *(gprimd(m1,iy1)*gprimd(m2,iy2) &
2115 968757516 : & -gprimd(m2,iy1)*gprimd(m1,iy2))
2116 : end do
2117 : end do
2118 : end do
2119 : !amet= -i amet
2120 128699670 : buffer1(:,:,:,:)=-amet(1,:,:,:,:)
2121 128699670 : amet(1,:,:,:,:) = amet(2,:,:,:,:)
2122 128699670 : amet(2,:,:,:,:) =buffer1(:,:,:,:)
2123 :
2124 : !DEBUG
2125 : !!Eventually construct the gradients of Amet wrt strains:
2126 : !! DAmet(y,y'',a,b,s,s'')= d[Amet(y,y'',s,s'')]/d[strain(a,b)]
2127 : !! -i Pauli(s,s'',n)*
2128 : !! ( E(n,a,m'')*Gprimd(b,y )*Gprimd(m'',y'')
2129 : !! +E(n,m,a )*Gprimd(b,y'')*Gprimd(m ,y ) )
2130 : !damet(:,:,:,:,:,:,:)=0.d0
2131 : !do ib=1,3
2132 : !do ia=1,3
2133 : !m1=mod(ia,3)+1 ! ia,m1,m2 is an even permutation
2134 : !m2=mod(m1,3)+1
2135 : !do iy2=1,3
2136 : !do iy1=1,3
2137 : !damet(:,iy1,iy2,ia,ib,:,:) = damet(:,iy1,iy2,ia,ib,:,:) &
2138 : !& + (pauli(:,:,:,m2)*gprimd(m1,iy2) &
2139 : !-pauli(:,:,:,m1)*gprimd(m2,iy2))*gprimd(ib,iy1) &
2140 : !& + (pauli(:,:,:,m1)*gprimd(m2,iy1) &
2141 : !-pauli(:,:,:,m2)*gprimd(m1,iy1))*gprimd(ib,iy2)
2142 : !end do
2143 : !end do
2144 : !end do
2145 : !end do
2146 : !! damet= i damet
2147 : !buffer2(:,:,:,:,:,:)= damet(1,:,:,:,:,:,:)
2148 : !damet(1,:,:,:,:,:,:)= -damet(2,:,:,:,:,:,:)
2149 : !damet(2,:,:,:,:,:,:)=buffer2(:,:,:,:,:,:)
2150 : !! Symetrize damet(:,:,:,a,b,:,:)
2151 : !damet(:,:,:,1,2,:,:)=0.5d0*(damet(:,:,:,1,2,:,:)+damet(:,:,:,2,1,:,:))
2152 : !damet(:,:,:,1,3,:,:)=0.5d0*(damet(:,:,:,1,3,:,:)+damet(:,:,:,3,1,:,:))
2153 : !damet(:,:,:,2,3,:,:)=0.5d0*(damet(:,:,:,2,3,:,:)+damet(:,:,:,3,2,:,:))
2154 : !damet(:,:,:,2,1,:,:)=damet(:,:,:,1,2,:,:)
2155 : !damet(:,:,:,3,1,:,:)=damet(:,:,:,1,3,:,:)
2156 : !damet(:,:,:,3,2,:,:)=damet(:,:,:,2,3,:,:)
2157 : !ENDDEBUG
2158 :
2159 2339994 : end subroutine metric_so
2160 : !!***
2161 :
2162 : end module m_contract
2163 : !!***
|
