Branch data 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 : 872029 : 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 : 872029 : & (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 : 872029 : & (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 : 872029 : & (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 : 872029 : & (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 : 872029 : & (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 : 872029 : & (rank1(re,3)*rank3(re,4)+rank1(im,3)*rank3(im,4)))
126 : :
127 : 872029 : 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 : 493901 : 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 : 493901 : & 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 : 493901 : & 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 : 493901 : & 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 : 493901 : & 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 : 493901 : & 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 : 493901 : & gmet(1,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6)))
234 : :
235 : 493901 : 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 : 636849 : 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 : 636849 : & 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 : 636849 : & 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 : 636849 : & 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 : 636849 : & 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 : 636849 : & 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 : 636849 : & gxa(re,6)*rank4(re, 2)+gxa(im,6)*rank4(im, 2))
590 : :
591 : 636849 : 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 : 293469401 : 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 [ + + ]: 293469401 : if (rank==0) then
1537 : : ! Simply copy scalar, re and im
1538 : 115515745 : bb(1,1)=aa(1,1)
1539 : 115515745 : bb(2,1)=aa(2,1)
1540 : :
1541 [ + + ]: 177953656 : else if (rank==1) then
1542 : : ! Apply gmet to input vector, re and im
1543 [ + + ]: 327047547 : do ii=1,2
1544 [ + + ]: 981142641 : do jj=1,3
1545 : 872126792 : 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 [ + + ]: 68937807 : else if (rank==2) then
1551 : : ! Apply rank2 expression, re and im
1552 [ + + ]: 204490266 : 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 : 136326844 : & 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 : 136326844 : & - 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 : 136326844 : & - 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 : 136326844 : & - 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 : 136326844 : & - 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 : 136326844 : & - 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 : 136326844 : & - gmet(2,1)*scalar)
1600 : : ! Include appropriate weights for multiplicity
1601 : 136326844 : bb(ii,4)=2.d0*bb(ii,4)
1602 : 136326844 : bb(ii,5)=2.d0*bb(ii,5)
1603 : 204490266 : 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 : 277228900 : 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 [ - + ]: 277228900 : 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 [ + + ]: 277228900 : else if (rank==1) then
1814 : : ! Apply gmet to input vector, re and im
1815 [ + + ]: 450186864 : do ii=1,2
1816 [ + + ]: 1350560592 : do jj=1,3
1817 : 1200498304 : 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 [ + + ]: 127166612 : else if (rank==2) then
1822 : : ! Apply rank2 expression, re and im
1823 [ + + ]: 378107388 : 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 : 252071592 : & (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 : 252071592 : & (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 : 252071592 : & (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 : 252071592 : & (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 : 252071592 : & (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 : 252071592 : & (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 : 252071592 : & (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 : 252071592 : & (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 : 252071592 : & (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 : 126035796 : 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 : 2424708 : 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 : 2424708 : pauli(:,:,:,:)=0.d0
2081 : 2424708 : pauli(1,1,2,1)= 1.d0;pauli(1,2,1,1)= 1.d0
2082 : 2424708 : pauli(2,1,2,2)=-1.d0;pauli(2,2,1,2)= 1.d0
2083 : 2424708 : pauli(1,1,1,3)= 1.d0;pauli(1,2,2,3)=-1.d0
2084 [ + + + + : 111536568 : pauli(:,:,:,:)= 0.5d0*pauli(:,:,:,:)
+ + + + ]
2085 : :
2086 [ + - + + ]: 2424708 : 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 [ + + + + : 2692227 : pauli(2,:,:,:) = aimag(S(:,:,:))
+ + ]
2103 : : end if
2104 : :
2105 : : !Construct the antisymmetric tensor:
2106 : 2424708 : amet(:,:,:,:,:)=0.d0
2107 [ + + ]: 9698832 : do iy2=1,3
2108 [ + + ]: 31521204 : do iy1=1,3
2109 [ + + ]: 94563612 : do n=1,3
2110 : 65467116 : m1=mod(n ,3)+1 ! n,m1,m2 is an even permutation
2111 : 65467116 : 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 [ + + + + : 1003829112 : & -gprimd(m2,iy1)*gprimd(m1,iy2))
+ + ]
2116 : : end do
2117 : : end do
2118 : : end do
2119 : : !amet= -i amet
2120 [ + + + + : 133358940 : buffer1(:,:,:,:)=-amet(1,:,:,:,:)
+ + + + ]
2121 [ + + + + : 133358940 : amet(1,:,:,:,:) = amet(2,:,:,:,:)
+ + + + ]
2122 [ + + + + : 133358940 : 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 : 2424708 : end subroutine metric_so
2160 : : !!***
2161 : :
2162 : : end module m_contract
2163 : : !!***
|
