TABLE OF CONTENTS


ABINIT/m_contract [ Modules ]

[ Top ] [ Modules ]

NAME

  m_contract

FUNCTION

 Low-level procedeures used in nonlop_pl to contract tensors

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (DCA, XG, MT, GZ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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

m_contract/cont13 [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont13

FUNCTION

 Contract rank1 tensor with rank3 symmetric tensor to
 produce symmetric rank2 tensor

INPUTS

  rank1(2,3)=rank 1 complex tensor (vector of length 3)
  rank3(2,10)=rank 3 complex tensor (symmetric storage)

OUTPUT

  rank2(6)=rank 2 real tensor (symmetric storage)

NOTES

 Tensors are in "symmetric" storage mode.
 For rank1 this is 1, 2, 3;
 for rank2 this is 11, 22, 33, 32, 31, 21;
 for rank3 this is 111, 221, 331, 321, 311, 211, 222, 332, 322, 333.
 rank1 and rank3 are complex; rank2 is real.
 Want $2 Re[contraction]$.
 $rank2(a,b)=2 Re[rank1(i)^"*" rank3(a,b,i)]$.
 In typical usage the input rank1 tensor is actually
 $rank1(i)=gmet(i,j) gxa(j)$

SOURCE

 77 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 & (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 & (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 & (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 & (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 & (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 & (rank1(re,3)*rank3(re,4)+rank1(im,3)*rank3(im,4)))
126 
127 end subroutine cont13

m_contract/cont22 [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont22

FUNCTION

 Contract symmetric rank 2 tensor gxa with itself using gmet to
 produce symmetric rank 2 tensor.

INPUTS

  gxa(2,6)=rank 2 complex tensor
  gmet(3,3)=real symmetric metric tensor (full storage)

OUTPUT

  rank2(6)=rank 2 real tensor (symmetric storage)

NOTES

 Symmetric gxa is stored as 11 22 33 32 31 21;
 gmet(3,3) is symmetric but stored fully (9 elements);
 output rank2 is stored as 11 22 33 32 31 21.
 Want $2 Re[contraction]$.
 $rank2(a,b)=2 Re[gxa(i,a)^"*" gmet(i,j) gxa(j,b)]$.

SOURCE

155 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 & 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 & 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 & 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 & 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 & 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 & gmet(1,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6)))
234 
235 end subroutine cont22

m_contract/cont22cso [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont22cso

FUNCTION

 Contract symmetric rank 2 tensor gxa1 with symmetric rank 2 tensor
 gxa2 using metric tensor gmet to produce rank 2 complex tensor.

INPUTS

  gxa1(2,10)=rank 2 complex symmetric tensor
  gxa2(2,10)=rank 2 complex symmetric tensor
  gmet(3,3)=usual metric tensor (symmetric, real)

OUTPUT

  rank2c(2,6)=rank 2 complex tensor (pseudo-symmetric storage)

NOTES

 This contraction is used for spin-orbit correction in non-local
 contribution to stresses.

 Symmetric gxa1, gxa2 are stored as 11 22 33 32 31 21;
 gmet(3,3) is symmetric but stored fully (9 elements);
 Output rank2c is not symmetric but since
      $rank2c_{gxa1,gxa2}(a,b)=conjg(rank2c_{gxa2,gxa1}(b,a))$
       it is stored as 11 22 33 32 31 21.

 rank2c(1,1), rank2c(2,2), rank3c(3,3) are not needed;
 They are not calculated.

 rank2c(a,b)=3 conjg(gxa1(i,a)) gmet(i,j) gxa2(j,b)

SOURCE

273 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  rank2c(:,:)=0.d0
290 
291 !First compute r2(i,j) = gmet(i,k) gxa2(j,k)
292  r2(:,1,1)=gmet(1,1)*gxa2(:,1)+gmet(1,2)*gxa2(:,6)+gmet(1,3)*gxa2(:,5)
293  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  r2(:,2,1)=gmet(2,1)*gxa2(:,1)+gmet(2,2)*gxa2(:,6)+gmet(2,3)*gxa2(:,5)
296  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  r2(:,3,1)=gmet(3,1)*gxa2(:,1)+gmet(3,2)*gxa2(:,6)+gmet(3,3)*gxa2(:,5)
299  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 & +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 & +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 & +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 & +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 & +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 & +gxa1(re,4)*r2(im,3,1)-gxa1(im,4)*r2(re,3,1))
340 
341 end subroutine cont22cso

m_contract/cont22so [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont22so

FUNCTION

 Contract symmetric rank 2 tensor gxa1 with symmetric rank 2 tensor
 gxa2 using antisymmetric tensor amet to produce rank 2 real tensor.

INPUTS

  gxa1(2,6)=rank 2 complex symmetric tensor
  gxa2(2,6)=rank 2 complex symmetric tensor
  amet(2,3,3)=antisymmetric complex tensor used for spin-orbit

OUTPUT

  rank2(6)=rank 2 real tensor (pseudo-symmetric storage)

NOTES

 This contraction is used for spin-orbit correction in non-local
 contribution to stresses.

 Symmetric gxa1, gxa2 are stored as 11 22 33 32 31 21;
 amet(3,3) is antisymmetric but stored fully (9 elements);
 Output rank2 is not symmetric but since
      $rank2_{gxa1,gxa2}(a,b)=conjg(rank2_{gxa2,gxa1}(b,a))$
       it is stored as 11 22 33 32 31 21.
 Want 2*Re[contraction].

 rank2(a,b)=2 Re[conjg(gxa1(i,a)) amet(i,j) gxa2(j,b)]

 Note that, since amet is antisymmetric, amet(i,i)=0

SOURCE

380 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 & 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 & 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 & 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 & 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 & 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 & amet(2,1,2)*(gxa1(re,6)*gxa2(im,6)-gxa1(im,6)*gxa2(re,6)))
483 
484 end subroutine cont22so

m_contract/cont24 [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont24

FUNCTION

 Contract symmetric rank2 tensor gxa with rank4 symmetric tensor to
 produce symmetric rank2 tensor.

INPUTS

  gxa(2,6)=rank 2 symmetric complex tensor in order 11 22 33 32 31 21
  rank4(2,15)=rank 4 complex tensor (symmetric storage)

OUTPUT

  rank2(6)=rank 2 real tensor (symmetric storage) 11 22 33 32 31 21.

NOTES

 Tensors are in "symmetric" storage mode.
 for gxa and rank2 this is 11, 22, 33, 32, 31, 21;
 for the rank 4 tensor rank4 this is
 1111 2211 3311 3211 3111 2111 2221 3321 3221 3331 2222 3322 3222 3332 3333.
 gxa and rank4 are complex; rank2 is real.
 Want $2 Re[contraction]$.
 $rank2(a,b)=2 Re[gxa(i,j)^"*" rank4(a,b,i,j)]$.

 Note that the input gxa is typically the result of
 $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)$
 where the subroutine "metcon" already includes weights in the
 definition of gxa for off-diagonal elements (weight of 2 for
 symmetry).
 Components 4, 5, and 6 of gxa have already been multiplied by 2
 so the expressions below do not carry the 2.

SOURCE

522 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 & 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 & 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 & 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 & 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 & 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 & gxa(re,6)*rank4(re, 2)+gxa(im,6)*rank4(im, 2))
590 
591 end subroutine cont24

m_contract/cont3 [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont3

FUNCTION

 Compute several specialized contractions needed for the
 l=3 part of the stress tensor.

INPUTS

  gxa(2,10)=complex symmetric rank 3 tensor
  gmet(3,3)=usual metric tensor, a symmetric matrix stored in
            full storage mode (bohr^-2)

OUTPUT

  rank2(6)=2*Re[contraction] given by
   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)]
   where r3(a,i,j)=gmet(j,k) gxa(a,i,k) and r1(a)=gmet(i,j) gxa(i,j,a).
  rank2 is stored in the compressed form 11 22 33 32 31 21.

NOTES

 Input gxa is a completely symmetric rank 3 tensor (complex)
 in compressed storage: 111 221 331 321 311 211 222 332 322 333.
 The output tensor is completely symmetric rank 2, real, and is given by
  $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)]$
  where $r3(a,i,j)=gmet(j,k) gxa(a,i,k)$ and $r1(a)=gmet(i,j) gxa(i,j,a)$.
  rank2 is stored in the compressed form 11 22 33 32 31 21.

SOURCE

623 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  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 &   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 &   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 &   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  do ii=1,2
665    r3(ii, 1)=gmet(1,1)*gxa(ii,1)+gmet(2,1)*gxa(ii,6)+&
666 &   gmet(3,1)*gxa(ii,5)
667    r3(ii, 2)=gmet(1,1)*gxa(ii,2)+gmet(2,1)*gxa(ii,7)+&
668 &   gmet(3,1)*gxa(ii,9)
669    r3(ii, 3)=gmet(1,1)*gxa(ii,3)+gmet(2,1)*gxa(ii,8)+&
670 &   gmet(3,1)*gxa(ii,10)
671    r3(ii, 4)=gmet(1,1)*gxa(ii,4)+gmet(2,1)*gxa(ii,9)+&
672 &   gmet(3,1)*gxa(ii,8)
673    r3(ii, 5)=gmet(1,1)*gxa(ii,5)+gmet(2,1)*gxa(ii,4)+&
674 &   gmet(3,1)*gxa(ii,3)
675    r3(ii, 6)=gmet(1,1)*gxa(ii,6)+gmet(2,1)*gxa(ii,2)+&
676 &   gmet(3,1)*gxa(ii,4)
677    r3(ii, 7)=gmet(2,1)*gxa(ii,1)+gmet(2,2)*gxa(ii,6)+&
678 &   gmet(3,2)*gxa(ii,5)
679    r3(ii, 8)=gmet(2,1)*gxa(ii,2)+gmet(2,2)*gxa(ii,7)+&
680 &   gmet(3,2)*gxa(ii,9)
681    r3(ii, 9)=gmet(2,1)*gxa(ii,3)+gmet(2,2)*gxa(ii,8)+&
682 &   gmet(3,2)*gxa(ii,10)
683    r3(ii,10)=gmet(2,1)*gxa(ii,4)+gmet(2,2)*gxa(ii,9)+&
684 &   gmet(3,2)*gxa(ii,8)
685    r3(ii,11)=gmet(2,1)*gxa(ii,5)+gmet(2,2)*gxa(ii,4)+&
686 &   gmet(3,2)*gxa(ii,3)
687    r3(ii,12)=gmet(2,1)*gxa(ii,6)+gmet(2,2)*gxa(ii,2)+&
688 &   gmet(3,2)*gxa(ii,4)
689    r3(ii,13)=gmet(3,1)*gxa(ii,1)+gmet(3,2)*gxa(ii,6)+&
690 &   gmet(3,3)*gxa(ii,5)
691    r3(ii,14)=gmet(3,1)*gxa(ii,2)+gmet(3,2)*gxa(ii,7)+&
692 &   gmet(3,3)*gxa(ii,9)
693    r3(ii,15)=gmet(3,1)*gxa(ii,3)+gmet(3,2)*gxa(ii,8)+&
694 &   gmet(3,3)*gxa(ii,10)
695    r3(ii,16)=gmet(3,1)*gxa(ii,4)+gmet(3,2)*gxa(ii,9)+&
696 &   gmet(3,3)*gxa(ii,8)
697    r3(ii,17)=gmet(3,1)*gxa(ii,5)+gmet(3,2)*gxa(ii,4)+&
698 &   gmet(3,3)*gxa(ii,3)
699    r3(ii,18)=gmet(3,1)*gxa(ii,6)+gmet(3,2)*gxa(ii,2)+&
700 &   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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & -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 & -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 & -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 & -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 & -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 & -3.d0*(r1(re,2)*r1(re,1)+r1(im,2)*r1(im,1))
806 
807 end subroutine cont3

m_contract/cont33cso [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont33cso

FUNCTION

 Contract symmetric rank 3 tensor gxa1 with symmetric rank 3 tensor
 gxa2 using metric tensor gmet to produce rank 2 complex tensor.

INPUTS

  gxa1(2,10)=rank 3 complex symmetric tensor
  gxa2(2,10)=rank 3 complex symmetric tensor
  gmet(3,3)=usual metric tensor (symmetric, real)

OUTPUT

  rank2c(2,6)=rank 2 complex tensor (pseudo-symmetric storage)

NOTES

 This contraction is used for spin-orbit correction in non-local
 contribution to stresses.

 Symmetric gxa1, gxa2 are stored as
               111 221 331 321 311 211 222 332 322 333;
 gmet(3,3) is symmetric but stored fully (9 elements);
 Output rank2c is not symmetric but since
      $rank2c_{gxa1,gxa2}(a,b)=conjg(rank2c_{gxa2,gxa1}(b,a))$
       it is stored as 11 22 33 32 31 21.

 rank2c(1,1), rank2c(2,2), rank3c(3,3) are not needed;
 They are not calculated.

 rank2c(a,b)=7.5 conjg(gxa1(a,i,j))*r_3(i,j,b) - 1.5 r_{11}(a)*r_{12}(b)
   where:
 r_3(i,j,b) & = & gxa2(b,l,m) gmet(i,l) gmet(j,m) \nonumber
 r_{11}(a)  & = & conjg(gxa1(a,l,m)) gmet(l,m)    \nonumber
 r_{12}(b)  & = & gxa2(b,l,m) gmet(l,m)
 \end{eqnarray} }}

SOURCE

 852 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  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  r3a(:,1,1)=gxa2(:,1 )*gmet(1,1)+gxa2(:,6 )*gmet(1,2)+gxa2(:,5 )*gmet(1,3)
 875  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  r3a(:,4,1)=gxa2(:,4 )*gmet(1,1)+gxa2(:,9 )*gmet(1,2)+gxa2(:,8 )*gmet(1,3)
 878  r3a(:,5,1)=gxa2(:,5 )*gmet(1,1)+gxa2(:,4 )*gmet(1,2)+gxa2(:,3 )*gmet(1,3)
 879  r3a(:,6,1)=gxa2(:,6 )*gmet(1,1)+gxa2(:,2 )*gmet(1,2)+gxa2(:,4 )*gmet(1,3)
 880  r3a(:,1,2)=gxa2(:,1 )*gmet(2,1)+gxa2(:,6 )*gmet(2,2)+gxa2(:,5 )*gmet(2,3)
 881  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  r3a(:,4,2)=gxa2(:,4 )*gmet(2,1)+gxa2(:,9 )*gmet(2,2)+gxa2(:,8 )*gmet(2,3)
 884  r3a(:,5,2)=gxa2(:,5 )*gmet(2,1)+gxa2(:,4 )*gmet(2,2)+gxa2(:,3 )*gmet(2,3)
 885  r3a(:,6,2)=gxa2(:,6 )*gmet(2,1)+gxa2(:,2 )*gmet(2,2)+gxa2(:,4 )*gmet(2,3)
 886  r3a(:,1,3)=gxa2(:,1 )*gmet(3,1)+gxa2(:,6 )*gmet(3,2)+gxa2(:,5 )*gmet(3,3)
 887  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  r3a(:,4,3)=gxa2(:,4 )*gmet(3,1)+gxa2(:,9 )*gmet(3,2)+gxa2(:,8 )*gmet(3,3)
 890  r3a(:,5,3)=gxa2(:,5 )*gmet(3,1)+gxa2(:,4 )*gmet(3,2)+gxa2(:,3 )*gmet(3,3)
 891  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  r3(:,1,1)=r3a(:,1,1)*gmet(1,1)+r3a(:,6,1)*gmet(1,2)+r3a(:,5,1)*gmet(1,3)
 896  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  r3(:,2,1)=r3a(:,1,2)*gmet(2,1)+r3a(:,6,2)*gmet(2,2)+r3a(:,5,2)*gmet(2,3)
 899  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  r3(:,3,1)=r3a(:,1,3)*gmet(3,1)+r3a(:,6,3)*gmet(3,2)+r3a(:,5,3)*gmet(3,3)
 902  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  r3(:,4,1)=r3a(:,1,2)*gmet(3,1)+r3a(:,6,2)*gmet(3,2)+r3a(:,5,2)*gmet(3,3)
 905  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  r3(:,5,1)=r3a(:,1,1)*gmet(3,1)+r3a(:,6,1)*gmet(3,2)+r3a(:,5,1)*gmet(3,3)
 908  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  r3(:,6,1)=r3a(:,1,1)*gmet(2,1)+r3a(:,6,1)*gmet(2,2)+r3a(:,5,1)*gmet(2,3)
 911  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 & +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 & +2.d0*(gxa1(:,8)*gmet(3,2)+gxa1(:,3)*gmet(3,1)+gxa1(:,4 )*gmet(2,1))
 921  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 & +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 & +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 & -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 & -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 & -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 & -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 & -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 & -1.5d0*(r11(re,2)*r12(im,1)+r11(im,2)*r12(re,1))
1016 
1017 end subroutine cont33cso

m_contract/cont33so [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont33so

FUNCTION

 Contract symmetric rank 3 tensor gxa1 with symmetric rank 3 tensor
 gxa2 using metric tensor gmet and antisymmetric tensor amet to
 produce rank 2 real tensor.

INPUTS

  gxa1(2,10)=rank 3 complex symmetric tensor
  gxa2(2,10)=rank 3 complex symmetric tensor
  gmet(3,3)=usual metric tensor (symmetric, real)
  amet(2,3,3)=antisymmetric complex tensor used for spin-orbit

OUTPUT

  rank2(6)=rank 2 real tensor (pseudo-symmetric storage)

NOTES

 This contraction is used for spin-orbit correction in non-local
 contribution to stresses.

 Symmetric gxa1, gxa2 are stored as
               111 221 331 321 311 211 222 332 322 333;
 gmet(3,3) is symmetric but stored fully (9 elements);
 amet(3,3) is antisymmetric but stored fully (9 elements);
 Output rank2 is not symmetric but since
      $rank2_{gxa1,gxa2}(a,b)=conjg(rank2_{gxa2,gxa1}(b,a))$
       it is stored as 11 22 33 32 31 21.
 Want 2*Re[contraction].

 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)]
   where:
   r_1(i)        & = & gxa2(i,j,k) gmet(j,k) \nonumber
   r_{3A}(i,j,k) & = & conjg(gxa1(p,i,j)) amet(p,k) \nonumber
   r_{3G}(i,j,k) & = & gxa2(p,i,j) gmet(p,k)
 \end{eqnarray} }}

SOURCE

1064 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  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 &   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 &   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 &   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  do ii=1,2
1104    r3G(ii, 1)=gmet(1,1)*gxa2(ii,1)+gmet(2,1)*gxa2(ii,6)+gmet(3,1)*gxa2(ii,5)
1105    r3G(ii, 2)=gmet(1,1)*gxa2(ii,2)+gmet(2,1)*gxa2(ii,7)+gmet(3,1)*gxa2(ii,9)
1106    r3G(ii, 3)=gmet(1,1)*gxa2(ii,3)+gmet(2,1)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,10)
1107    r3G(ii, 4)=gmet(1,1)*gxa2(ii,4)+gmet(2,1)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,8)
1108    r3G(ii, 5)=gmet(1,1)*gxa2(ii,5)+gmet(2,1)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,3)
1109    r3G(ii, 6)=gmet(1,1)*gxa2(ii,6)+gmet(2,1)*gxa2(ii,2)+gmet(3,1)*gxa2(ii,4)
1110    r3G(ii, 7)=gmet(2,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,5)
1111    r3G(ii, 8)=gmet(2,1)*gxa2(ii,2)+gmet(2,2)*gxa2(ii,7)+gmet(3,2)*gxa2(ii,9)
1112    r3G(ii, 9)=gmet(2,1)*gxa2(ii,3)+gmet(2,2)*gxa2(ii,8)+gmet(3,2)*gxa2(ii,10)
1113    r3G(ii,10)=gmet(2,1)*gxa2(ii,4)+gmet(2,2)*gxa2(ii,9)+gmet(3,2)*gxa2(ii,8)
1114    r3G(ii,11)=gmet(2,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,3)
1115    r3G(ii,12)=gmet(2,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,4)
1116    r3G(ii,13)=gmet(3,1)*gxa2(ii,1)+gmet(3,2)*gxa2(ii,6)+gmet(3,3)*gxa2(ii,5)
1117    r3G(ii,14)=gmet(3,1)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,7)+gmet(3,3)*gxa2(ii,9)
1118    r3G(ii,15)=gmet(3,1)*gxa2(ii,3)+gmet(3,2)*gxa2(ii,8)+gmet(3,3)*gxa2(ii,10)
1119    r3G(ii,16)=gmet(3,1)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,9)+gmet(3,3)*gxa2(ii,8)
1120    r3G(ii,17)=gmet(3,1)*gxa2(ii,5)+gmet(3,2)*gxa2(ii,4)+gmet(3,3)*gxa2(ii,3)
1121    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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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 & 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  rank2(:)=15.d0*s33(:)-3.d0*s31(:)
1287 
1288 end subroutine cont33so

m_contract/cont35 [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 cont35

FUNCTION

 Contract symmetric rank3 tensor gxa with rank5 symmetric tensor to
 produce symmetric rank2 tensor.

INPUTS

  gxa(2,10)=rank 3 symmetric complex tensor in order
  rank5(2,21)=rank 5 complex tensor (symmetric storage)

OUTPUT

  rank2(6)=rank 2 real tensor (symmetric storage) 11 22 33 32 31 21.

NOTES

 Tensors are in "symmetric" storage mode.
 For rank 3 tensor gxa this is
     111 221 331 321 311 211 222 332 322 333;
 For rank 5 tensor rank5 this is
     11111 22111 33111 32111 31111 21111 22211 33211 32211 33311
     22221 33221 32221 33321 33331 22222 33222 32222 33322 33332 33333;
 For rank 2 tensor rank2 this is 11, 22, 33, 32, 31, 21;
 gxa and rank5 are complex; rank2 is real.
 Want $2 Re[contraction]$.
 $rank2(a,b)=2 Re[gxa(i,j,k)^"*" rank5(a,b,i,j,k)]$.

 Note that the input gxa is typically the result of
 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)
 where the subroutine "metcon" already includes weights in the definition
 of gxa for off-diagonal elements.

SOURCE

1327 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 & 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 & 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 & 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 & 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 & 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 & gxa(re, 4)*rank5(re, 9)+gxa(im, 4)*rank5(im, 9))
1420 
1421 end subroutine cont35

m_contract/metcon [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 metcon

FUNCTION

 Carries out specialized metric tensor contractions needed for
 l=0,1,2,3 nonlocal Kleinman-Bylander pseudopotential operation.
 Full advantage is taken of the full permutational symmetry of these
 tensors.

INPUTS

  rank=0,1,2, or 3 = rank of input tensor aa
  gmet(3,3)=metric tensor (array is symmetric but stored as 3x3)
  aa(2,(rank+1)*(rank+2)/2)=unique elements of complex input tensor

OUTPUT

  bb(2,(rank+1)*(rank+2)/2)=unique elements of complex output tensor

NOTES

 All tensors are stored in a compressed storage mode defined below;
 input and output conform to this scheme.
 When tensor elements occur repeatedly due to symmetry, the
 WEIGHT IS INCLUDED in the output tensor element to simplify later
 contractions with other tensors of the same rank and form, i.e. the
 next contraction is then simply a dot product over the unique elements.

 Definitions of the contractions:

 rank=0:  bb = aa  (simply copy the scalar value, no use of gmet)

 rank=1:  bb(i)= $gmet(i,l) aa(l)$  (3 elements in, 3 elements out)

 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)$
          (6 elements in, 6 elements out)

 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)$
          (10 elements in, 10 elements out)
         In this rank 3 case, the second term is NOT symmetric in all
         permutations of i,j,k, but the final tensor b(ijk) may be
         symmetrized over all permutations because it will be
         contracted with a completely symmetric tensor.

 The compressed storage scheme is based on storing a symmetric 3x3 matrix as
     (1 . .)
     (6 2 .)
     (5 4 3)
 which leads to the following mappings for all ranks
 where the compressed storage index is to the right of the arrow:
 rank=0 1->1 (only a scalar)
 rank=1 1->1 2->2 3->3 (standard vector, no compression)
 rank=2 11->1 22->2 33->3 32->4 31->5 21->6
  weights   1     1     1     2     2     2
 rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
  weights    1      3      3      6      3      3      1      3      3       1

SOURCE

1481 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  if (rank==0) then
1537 !  Simply copy scalar, re and im
1538    bb(1,1)=aa(1,1)
1539    bb(2,1)=aa(2,1)
1540 
1541  else if (rank==1) then
1542 !  Apply gmet to input vector, re and im
1543    do ii=1,2
1544      do jj=1,3
1545        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  else if (rank==2) then
1551 !  Apply rank2 expression, re and im
1552    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 &     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 &     - 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 &     - 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 &     - 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 &     - 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 &     - 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 &     - gmet(2,1)*scalar)
1600 !    Include appropriate weights for multiplicity
1601      bb(ii,4)=2.d0*bb(ii,4)
1602      bb(ii,5)=2.d0*bb(ii,5)
1603      bb(ii,6)=2.d0*bb(ii,6)
1604    end do
1605 
1606  else if (rank==3) then
1607 !  Apply rank2 expression, re and im
1608    do ii=1,2
1609 !    Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j)
1610      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 &       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 &       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 &       gmet(3,3)*gmet(jj,2)) *aa(ii,8)
1639      end do
1640 !    Write out components of contraction
1641 !    (111)->1
1642      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 &     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 &     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 &     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 &     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 &     gmet(1,1)*vector(2)+gmet(2,1)*vector(1))
1658 !    (222)->7
1659      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 &     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 &     gmet(2,2)*vector(3)+gmet(3,2)*vector(2))
1667 !    (333)->10
1668      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      bb(ii,2)=3.d0*bb(ii,2)
1671      bb(ii,3)=3.d0*bb(ii,3)
1672      bb(ii,4)=6.d0*bb(ii,4)
1673      bb(ii,5)=3.d0*bb(ii,5)
1674      bb(ii,6)=3.d0*bb(ii,6)
1675      bb(ii,8)=3.d0*bb(ii,8)
1676      bb(ii,9)=3.d0*bb(ii,9)
1677    end do
1678 
1679  else
1680    write(message, '(a,i0,a,a,a)' )&
1681 &   'Input rank=',rank,' not allowed.',ch10,&
1682 &   'Possible values are 0,1,2,3 only.'
1683    ABI_BUG(message)
1684  end if
1685 
1686  contains
1687 
1688    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 &   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 &   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 &   gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8)
1723 
1724  end function con_met
1725 
1726 end subroutine metcon

m_contract/metcon_so [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 metcon_so

FUNCTION

 Carries out specialized metric tensor contractions needed for
 l=0,1,2,3 nonlocal Kleinman-Bylander pseudopotential operation
 in the spin-orbit case.
 Full advantage is taken of the full permutational symmetry of these tensors.

INPUTS

  rank=0,1,2, or 3 = rank of input tensor aa
  gmet(3,3)=metric tensor (array is symmetric but stored as 3x3)
  amet(3,3)=real or imaginary part of one spin matrix element of the
            "spin metric" tensor
  aa(2,(rank+1)*(rank+2)/2)=unique elements of complex input tensor

OUTPUT

  bb(2,(rank+1)*(rank+2)/2)=unique elements of complex output tensor

NOTES

 All tensors are stored in a compressed storage mode defined below;
 input and output conform to this scheme.
 When tensor elements occur repeatedly due to symmetry, the
 WEIGHT IS INCLUDED in the output tensor element to simplify later
 contractions with other tensors of the same rank and form, i.e. the
 next contraction is then simply a dot product over the unique elements.

 Definitions of the contractions:

 rank=0:  bb=0

 rank=1:  bb(i)= $amet(i,l) aa(l)$  (3 elements in, 3 elements out)

 rank=2:  bb(i,j)= $[3 gmet(i,l) amet(j,m)] aa(l,m)$
          (6 elements in, 6 elements out)

 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)$
          (10 elements in, 10 elements out)
          In this rank 3 case, the second term is NOT symmetric in all
          permutations of i,j,k, but the final tensor b(ijk) may be
          symmetrized over all permutations because it will be
          contracted with a completely symmetric tensor.

 The compressed storage scheme is based on storing
 a symmetric 3x3 matrix as
     (1 . .)
     (6 2 .)
     (5 4 3)
 which leads to the following mappings for all ranks
 where the compressed storage index is to the right of the arrow:
 rank=0 1->1 (only a scalar)
 rank=1 1->1 2->2 3->3 (standard vector, no compression)
 rank=2 11->1 22->2 33->3 32->4 31->5 21->6
  weights   1     1     1     2     2     2
 rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
  weights    1      3      3      6      3      3      1      3      3       1

SOURCE

1789 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  if (rank==0) then
1809 !  Simply copy scalar, re and im
1810    bb(1,1)=0.d0
1811    bb(2,1)=0.d0
1812 !
1813  else if (rank==1) then
1814 !  Apply gmet to input vector, re and im
1815    do ii=1,2
1816      do jj=1,3
1817        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  else if (rank==2) then
1822 !  Apply rank2 expression, re and im
1823    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 &     (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 &     (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 &     (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 &     (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 &     (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 &     (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 &     (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 &     (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 &     (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      bb(ii,6)=bb(ii,6)
1880    end do
1881 
1882  else if (rank==3) then
1883 !  Apply rank2 expression, re and im
1884    do ii=1,2
1885 !    Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j)
1886      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 &       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 &       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 &       gmet(3,3)*amet(jj,2)) *aa(ii,8)
1915      end do
1916 !    Write out components of contraction
1917 !    (111)->1
1918      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 &     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 &     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 &     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 &     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 &     gmet(1,1)*vector(2)+gmet(2,1)*vector(1))
1934 !    (222)->7
1935      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 &     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 &     gmet(2,2)*vector(3)+gmet(3,2)*vector(2))
1943 !    (333)->10
1944      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      bb(ii,2)=3.d0*bb(ii,2)
1947      bb(ii,3)=3.d0*bb(ii,3)
1948      bb(ii,4)=6.d0*bb(ii,4)
1949      bb(ii,5)=3.d0*bb(ii,5)
1950      bb(ii,6)=3.d0*bb(ii,6)
1951      bb(ii,8)=3.d0*bb(ii,8)
1952      bb(ii,9)=3.d0*bb(ii,9)
1953    end do
1954 
1955  else
1956    write(message, '(a,i0,a,a,a)' )&
1957 &   'Input rank=',rank,' not allowed.',ch10,&
1958 &   'Possible values are 0,1,2,3 only.'
1959    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    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 &   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 &   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 &   gmet(i1,3)*gmet(i2,3)*amet(i3,2))*aa(ii,8)
2002  end function cona_metso
2003 
2004 
2005    function con_metso(ii,i1,i2,i3)
2006 
2007    real(dp) :: con_metso
2008    integer,intent(in) :: ii,i1,i2,i3
2009 
2010    con_metso=(cona_metso(ii,i3,i1,i2)+cona_metso(ii,i2,i3,i1)+cona_metso(ii,i1,i2,i3))/3.d0
2011 
2012  end function con_metso
2013 
2014 end subroutine metcon_so

m_contract/metric_so [ Functions ]

[ Top ] [ m_contract ] [ Functions ]

NAME

 metric_so

FUNCTION

 Computes Pauli matrices and antisymmetric tensor needed for
 spin-orbit.

INPUTS

  gprimd(3,3)=dimensional primitive translations for reciprocal space
              (bohr**-1)

OUTPUT

  amet(2,3,3,2,2)=the antisymmetric tensor A(Re/Im,y,y'',s,s'')
  pauli(2,2,2,3)=Pauli matrixes

NOTES

 (the preprocessing based on CPP does not allow isolated quotes,
  so that they have been doubled in the following latex equations)
 $2 Pauli(s,s'',1) & = & 0 & 1 \nonumber
                   &   & 1 & 0 \nonumber
  2 Pauli(s,s'',2) & = & 0 & -i \nonumber
                   &   & i &  0 \nonumber
  2 Pauli(s,s'',3) & = & 1 &  0 \nonumber
                   &   & 0 & -1
 \end{eqnarray} }}
    $Amet(y,y'',s,s'') & = & -i Pauli(s,s'',n) E(n,m,m'') Gprimd(m,y) Gprimd(m'',y'')
 \end{eqnarray} }}

 E(n,m,m''):   full antisymetric tensor
 s,s'':        spin indices (1..2)
 y,y'',m,m'',n: metric indices (1..3)
 a,b:         strain indices (1..3)
 Amet and Pauli are complex

SOURCE

2056 subroutine metric_so(amet,gprimd,pauli)
2057 
2058 !Arguments ------------------------------------
2059 !arrays
2060  real(dp),intent(in) :: gprimd(3,3)
2061  real(dp),intent(out) :: amet(2,3,3,2,2),pauli(2,2,2,3)
2062 
2063 !Local variables-------------------------------
2064 !scalars
2065  integer :: iy1,iy2,m1,m2,n
2066 !arrays
2067  real(dp) :: buffer1(3,3,2,2) !,buffer2(3,3,3,3,2,2)
2068 
2069 ! **********************************************************************
2070 
2071 !Fill in Pauli matrices and make them spin matrices:
2072 !Pauli(Re/Im,up,down,coord):
2073  pauli(:,:,:,:)=0.d0
2074  pauli(1,1,2,1)= 1.d0;pauli(1,2,1,1)= 1.d0
2075  pauli(2,1,2,2)=-1.d0;pauli(2,2,1,2)= 1.d0
2076  pauli(1,1,1,3)= 1.d0;pauli(1,2,2,3)=-1.d0
2077  pauli(:,:,:,:)= 0.5d0*pauli(:,:,:,:)
2078 
2079 !Construct the antisymmetric tensor:
2080  amet(:,:,:,:,:)=0.d0
2081  do iy2=1,3
2082    do iy1=1,3
2083      do n=1,3
2084        m1=mod(n ,3)+1    !  n,m1,m2 is an even permutation
2085        m2=mod(m1,3)+1
2086        amet(1:2,iy1,iy2,1:2,1:2) = amet(:,iy1,iy2,:,:) &
2087 &       + pauli(:,:,:,n) &
2088 &       *(gprimd(m1,iy1)*gprimd(m2,iy2) &
2089 &       -gprimd(m2,iy1)*gprimd(m1,iy2))
2090      end do
2091    end do
2092  end do
2093 !amet= -i amet
2094  buffer1(:,:,:,:)=-amet(1,:,:,:,:)
2095  amet(1,:,:,:,:) = amet(2,:,:,:,:)
2096  amet(2,:,:,:,:) =buffer1(:,:,:,:)
2097 
2098 !DEBUG
2099 !!Eventually construct the gradients of Amet wrt strains:
2100 !! DAmet(y,y'',a,b,s,s'')= d[Amet(y,y'',s,s'')]/d[strain(a,b)]
2101 !!                        -i Pauli(s,s'',n)*
2102 !!                          ( E(n,a,m'')*Gprimd(b,y )*Gprimd(m'',y'')
2103 !!                           +E(n,m,a )*Gprimd(b,y'')*Gprimd(m ,y ) )
2104 !damet(:,:,:,:,:,:,:)=0.d0
2105 !do ib=1,3
2106 !do ia=1,3
2107 !m1=mod(ia,3)+1    !  ia,m1,m2 is an even permutation
2108 !m2=mod(m1,3)+1
2109 !do iy2=1,3
2110 !do iy1=1,3
2111 !damet(:,iy1,iy2,ia,ib,:,:) = damet(:,iy1,iy2,ia,ib,:,:) &
2112 !&        + (pauli(:,:,:,m2)*gprimd(m1,iy2) &
2113 !-pauli(:,:,:,m1)*gprimd(m2,iy2))*gprimd(ib,iy1) &
2114 !&        + (pauli(:,:,:,m1)*gprimd(m2,iy1) &
2115 !-pauli(:,:,:,m2)*gprimd(m1,iy1))*gprimd(ib,iy2)
2116 !end do
2117 !end do
2118 !end do
2119 !end do
2120 !! damet= i damet
2121 !buffer2(:,:,:,:,:,:)= damet(1,:,:,:,:,:,:)
2122 !damet(1,:,:,:,:,:,:)= -damet(2,:,:,:,:,:,:)
2123 !damet(2,:,:,:,:,:,:)=buffer2(:,:,:,:,:,:)
2124 !! Symetrize damet(:,:,:,a,b,:,:)
2125 !damet(:,:,:,1,2,:,:)=0.5d0*(damet(:,:,:,1,2,:,:)+damet(:,:,:,2,1,:,:))
2126 !damet(:,:,:,1,3,:,:)=0.5d0*(damet(:,:,:,1,3,:,:)+damet(:,:,:,3,1,:,:))
2127 !damet(:,:,:,2,3,:,:)=0.5d0*(damet(:,:,:,2,3,:,:)+damet(:,:,:,3,2,:,:))
2128 !damet(:,:,:,2,1,:,:)=damet(:,:,:,1,2,:,:)
2129 !damet(:,:,:,3,1,:,:)=damet(:,:,:,1,3,:,:)
2130 !damet(:,:,:,3,2,:,:)=damet(:,:,:,2,3,:,:)
2131 !ENDDEBUG
2132 
2133 end subroutine metric_so