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-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_contract
28 
29  use defs_basis
30  use m_errors
31 
32  implicit none
33 
34  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)$

PARENTS

      nonlop_pl

CHILDREN

SOURCE

 87 subroutine cont13(rank1,rank3,rank2)
 88 
 89 
 90 !This section has been created automatically by the script Abilint (TD).
 91 !Do not modify the following lines by hand.
 92 #undef ABI_FUNC
 93 #define ABI_FUNC 'cont13'
 94 !End of the abilint section
 95 
 96  implicit none
 97 
 98 !Arguments ------------------------------------
 99 !arrays
100  real(dp),intent(in) :: rank1(2,3),rank3(2,10)
101  real(dp),intent(out) :: rank2(6)
102 
103 !Local variables-------------------------------
104 !scalars
105  integer,parameter :: im=2,re=1
106 
107 ! *************************************************************************
108 
109 !Simply write out index summations
110 !a=1, b=1 in rank2(a,b) --> maps to index 1
111  rank2(1)=2.0d0*(&
112 & (rank1(re,1)*rank3(re,1)+rank1(im,1)*rank3(im,1))+&
113 & (rank1(re,2)*rank3(re,6)+rank1(im,2)*rank3(im,6))+&
114 & (rank1(re,3)*rank3(re,5)+rank1(im,3)*rank3(im,5)))
115 
116 !a=2, b=2 in rank2(a,b) --> maps to index 2
117  rank2(2)=2.0d0*(&
118 & (rank1(re,1)*rank3(re,2)+rank1(im,1)*rank3(im,2))+&
119 & (rank1(re,2)*rank3(re,7)+rank1(im,2)*rank3(im,7))+&
120 & (rank1(re,3)*rank3(re,9)+rank1(im,3)*rank3(im,9)))
121 
122 !a=3, b=3 in rank2(a,b) --> maps to index 3
123  rank2(3)=2.0d0*(&
124 & (rank1(re,1)*rank3(re,3)+rank1(im,1)*rank3(im,3))+&
125 & (rank1(re,2)*rank3(re,8)+rank1(im,2)*rank3(im,8))+&
126 & (rank1(re,3)*rank3(re,10)+rank1(im,3)*rank3(im,10)))
127 
128 !a=3, b=2 in rank2(a,b) --> maps to index 4
129  rank2(4)=2.0d0*(&
130 & (rank1(re,1)*rank3(re,4)+rank1(im,1)*rank3(im,4))+&
131 & (rank1(re,2)*rank3(re,9)+rank1(im,2)*rank3(im,9))+&
132 & (rank1(re,3)*rank3(re,8)+rank1(im,3)*rank3(im,8)))
133 
134 !a=3, b=1 in rank2(a,b) --> maps to index 5
135  rank2(5)=2.0d0*(&
136 & (rank1(re,1)*rank3(re,5)+rank1(im,1)*rank3(im,5))+&
137 & (rank1(re,2)*rank3(re,4)+rank1(im,2)*rank3(im,4))+&
138 & (rank1(re,3)*rank3(re,3)+rank1(im,3)*rank3(im,3)))
139 
140 !a=2, b=1 in rank2(a,b) --> maps to index 6
141  rank2(6)=2.0d0*(&
142 & (rank1(re,1)*rank3(re,6)+rank1(im,1)*rank3(im,6))+&
143 & (rank1(re,2)*rank3(re,2)+rank1(im,2)*rank3(im,2))+&
144 & (rank1(re,3)*rank3(re,4)+rank1(im,3)*rank3(im,4)))
145 
146 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)]$.

PARENTS

      nonlop_pl

CHILDREN

SOURCE

179 subroutine cont22(gxa,gmet,rank2)
180 
181 
182 !This section has been created automatically by the script Abilint (TD).
183 !Do not modify the following lines by hand.
184 #undef ABI_FUNC
185 #define ABI_FUNC 'cont22'
186 !End of the abilint section
187 
188  implicit none
189 
190 !Arguments ------------------------------------
191 !arrays
192  real(dp),intent(in) :: gmet(3,3),gxa(2,6)
193  real(dp),intent(out) :: rank2(6)
194 
195 !Local variables-------------------------------
196 !scalars
197  integer,parameter :: im=2,re=1
198 
199 ! *************************************************************************
200 
201 !Simply write out index summations
202 !a=1, b=1 in rank2(a,b) --> maps to index 1
203  rank2(1)=2.0d0*(&
204 & gmet(1,1)*(gxa(re,1)*gxa(re,1)+gxa(im,1)*gxa(im,1))+&
205 & gmet(2,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6))+&
206 & gmet(3,3)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+&
207 & 2.0d0*(&
208 & gmet(3,2)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6))+&
209 & gmet(3,1)*(gxa(re,5)*gxa(re,1)+gxa(im,5)*gxa(im,1))+&
210 & gmet(2,1)*(gxa(re,6)*gxa(re,1)+gxa(im,6)*gxa(im,1))))
211 
212 !a=2, b=2 in rank2(a,b) --> maps to index 2
213  rank2(2)=2.0d0*(&
214 & gmet(1,1)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6))+&
215 & gmet(2,2)*(gxa(re,2)*gxa(re,2)+gxa(im,2)*gxa(im,2))+&
216 & gmet(3,3)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+&
217 & 2.0d0*(&
218 & gmet(3,2)*(gxa(re,4)*gxa(re,2)+gxa(im,4)*gxa(im,2))+&
219 & gmet(3,1)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
220 & gmet(2,1)*(gxa(re,2)*gxa(re,6)+gxa(im,2)*gxa(im,6))))
221 
222 !a=3, b=3 in rank2(a,b) --> maps to index 3
223  rank2(3)=2.0d0*(&
224 & gmet(1,1)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+&
225 & gmet(2,2)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+&
226 & gmet(3,3)*(gxa(re,3)*gxa(re,3)+gxa(im,3)*gxa(im,3))+&
227 & 2.0d0*(&
228 & gmet(3,2)*(gxa(re,4)*gxa(re,3)+gxa(im,4)*gxa(im,3))+&
229 & gmet(3,1)*(gxa(re,5)*gxa(re,3)+gxa(im,5)*gxa(im,3))+&
230 & gmet(2,1)*(gxa(re,5)*gxa(re,4)+gxa(im,5)*gxa(im,4))))
231 
232 !a=3, b=2 in rank2(a,b) --> maps to index 4
233  rank2(4)=2.0d0*(&
234 & gmet(1,1)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6))+&
235 & gmet(2,2)*(gxa(re,4)*gxa(re,2)+gxa(im,4)*gxa(im,2))+&
236 & gmet(3,3)*(gxa(re,3)*gxa(re,4)+gxa(im,3)*gxa(im,4))+&
237 & gmet(3,2)*(gxa(re,3)*gxa(re,2)+gxa(im,3)*gxa(im,2))+&
238 & gmet(3,1)*(gxa(re,3)*gxa(re,6)+gxa(im,3)*gxa(im,6))+&
239 & gmet(2,1)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
240 & gmet(2,3)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+&
241 & gmet(1,3)*(gxa(re,5)*gxa(re,4)+gxa(im,5)*gxa(im,4))+&
242 & gmet(1,2)*(gxa(re,5)*gxa(re,2)+gxa(im,5)*gxa(im,2)))
243 
244 !a=3, b=1 in rank2(a,b) --> maps to index 5
245  rank2(5)=2.0d0*(&
246 & gmet(1,1)*(gxa(re,5)*gxa(re,1)+gxa(im,5)*gxa(im,1))+&
247 & gmet(2,2)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
248 & gmet(3,3)*(gxa(re,3)*gxa(re,5)+gxa(im,3)*gxa(im,5))+&
249 & gmet(3,2)*(gxa(re,3)*gxa(re,6)+gxa(im,3)*gxa(im,6))+&
250 & gmet(3,1)*(gxa(re,3)*gxa(re,1)+gxa(im,3)*gxa(im,1))+&
251 & gmet(2,1)*(gxa(re,4)*gxa(re,1)+gxa(im,4)*gxa(im,1))+&
252 & gmet(2,3)*(gxa(re,4)*gxa(re,5)+gxa(im,4)*gxa(im,5))+&
253 & gmet(1,3)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+&
254 & gmet(1,2)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6)))
255 
256 !a=2, b=1 in rank2(a,b) --> maps to index 6
257  rank2(6)=2.0d0*(&
258 & gmet(1,1)*(gxa(re,6)*gxa(re,1)+gxa(im,6)*gxa(im,1))+&
259 & gmet(2,2)*(gxa(re,2)*gxa(re,6)+gxa(im,2)*gxa(im,6))+&
260 & gmet(3,3)*(gxa(re,4)*gxa(re,5)+gxa(im,4)*gxa(im,5))+&
261 & gmet(3,2)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+&
262 & gmet(3,1)*(gxa(re,4)*gxa(re,1)+gxa(im,4)*gxa(im,1))+&
263 & gmet(2,1)*(gxa(re,2)*gxa(re,1)+gxa(im,2)*gxa(im,1))+&
264 & gmet(2,3)*(gxa(re,2)*gxa(re,5)+gxa(im,2)*gxa(im,5))+&
265 & gmet(1,3)*(gxa(re,6)*gxa(re,5)+gxa(im,6)*gxa(im,5))+&
266 & gmet(1,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6)))
267 
268 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)

PARENTS

      nonlop_pl

CHILDREN

SOURCE

311 subroutine cont22cso(gxa1,gxa2,gmet,rank2c)
312 
313 
314 !This section has been created automatically by the script Abilint (TD).
315 !Do not modify the following lines by hand.
316 #undef ABI_FUNC
317 #define ABI_FUNC 'cont22cso'
318 !End of the abilint section
319 
320  implicit none
321 
322 !Arguments ------------------------------------
323 !arrays
324  real(dp),intent(in) :: gmet(3,3),gxa1(2,6),gxa2(2,6)
325  real(dp),intent(out) :: rank2c(2,6)
326 
327 !Local variables-------------------------------
328 !scalars
329  integer,parameter :: im=2,re=1
330 !arrays
331  real(dp) :: r2(2,3,3)
332 
333 ! *************************************************************************
334 
335 !Initialize output tensor
336  rank2c(:,:)=0.d0
337 
338 !First compute r2(i,j) = gmet(i,k) gxa2(j,k)
339  r2(:,1,1)=gmet(1,1)*gxa2(:,1)+gmet(1,2)*gxa2(:,6)+gmet(1,3)*gxa2(:,5)
340  r2(:,1,2)=gmet(1,1)*gxa2(:,6)+gmet(1,2)*gxa2(:,2)+gmet(1,3)*gxa2(:,4)
341  r2(:,1,3)=gmet(1,1)*gxa2(:,5)+gmet(1,2)*gxa2(:,4)+gmet(1,3)*gxa2(:,3)
342  r2(:,2,1)=gmet(2,1)*gxa2(:,1)+gmet(2,2)*gxa2(:,6)+gmet(2,3)*gxa2(:,5)
343  r2(:,2,2)=gmet(2,1)*gxa2(:,6)+gmet(2,2)*gxa2(:,2)+gmet(2,3)*gxa2(:,4)
344  r2(:,2,3)=gmet(2,1)*gxa2(:,5)+gmet(2,2)*gxa2(:,4)+gmet(2,3)*gxa2(:,3)
345  r2(:,3,1)=gmet(3,1)*gxa2(:,1)+gmet(3,2)*gxa2(:,6)+gmet(3,3)*gxa2(:,5)
346  r2(:,3,2)=gmet(3,1)*gxa2(:,6)+gmet(3,2)*gxa2(:,2)+gmet(3,3)*gxa2(:,4)
347  r2(:,3,3)=gmet(3,1)*gxa2(:,5)+gmet(3,2)*gxa2(:,4)+gmet(3,3)*gxa2(:,3)
348 
349 !Then compute rank2c(a,b) = 3 conjg(gxa1(a,i)) r2(i,b)
350 !stored as 11 22 33 32 31 21
351 !rank2c(re,1) = 3.d0*(gxa1(re,1)*r2(re,1,1)+gxa1(im,1)*r2(im,1,1)&
352 !&                     +gxa1(re,6)*r2(re,2,1)+gxa1(im,6)*r2(im,2,1)&
353 !&                     +gxa1(re,5)*r2(re,3,1)+gxa1(im,5)*r2(im,3,1))
354 !rank2c(re,2) = 3.d0*(gxa1(re,6)*r2(re,1,2)+gxa1(im,6)*r2(im,1,2)&
355 !&                     +gxa1(re,2)*r2(re,2,2)+gxa1(im,2)*r2(im,2,2)&
356 !&                     +gxa1(re,4)*r2(re,3,2)+gxa1(im,4)*r2(im,3,2))
357 !rank2c(re,3) = 3.d0*(gxa1(re,5)*r2(re,1,3)+gxa1(im,5)*r2(im,1,3)&
358 !&                     +gxa1(re,4)*r2(re,2,3)+gxa1(im,4)*r2(im,2,3)&
359 !&                     +gxa1(re,3)*r2(re,3,3)+gxa1(im,3)*r2(im,3,3))
360  rank2c(re,4) = 3.d0*(gxa1(re,5)*r2(re,1,2)+gxa1(im,5)*r2(im,1,2)&
361 & +gxa1(re,4)*r2(re,2,2)+gxa1(im,4)*r2(im,2,2)&
362 & +gxa1(re,3)*r2(re,3,2)+gxa1(im,3)*r2(im,3,2))
363  rank2c(re,5) = 3.d0*(gxa1(re,5)*r2(re,1,1)+gxa1(im,5)*r2(im,1,1)&
364 & +gxa1(re,4)*r2(re,2,1)+gxa1(im,4)*r2(im,2,1)&
365 & +gxa1(re,3)*r2(re,3,1)+gxa1(im,3)*r2(im,3,1))
366  rank2c(re,6) = 3.d0*(gxa1(re,6)*r2(re,1,1)+gxa1(im,6)*r2(im,1,1)&
367 & +gxa1(re,2)*r2(re,2,1)+gxa1(im,2)*r2(im,2,1)&
368 & +gxa1(re,4)*r2(re,3,1)+gxa1(im,4)*r2(im,3,1))
369 !rank2c(im,1) = 3.d0*(gxa1(re,1)*r2(im,1,1)-gxa1(im,1)*r2(re,1,1)&
370 !&                     +gxa1(re,6)*r2(im,2,1)-gxa1(im,6)*r2(re,2,1)&
371 !&                     +gxa1(re,5)*r2(im,3,1)-gxa1(im,5)*r2(re,3,1))
372 !rank2c(im,2) = 3.d0*(gxa1(re,6)*r2(im,1,2)-gxa1(im,6)*r2(re,1,2)&
373 !&                     +gxa1(re,2)*r2(im,2,2)-gxa1(im,2)*r2(re,2,2)&
374 !&                     +gxa1(re,4)*r2(im,3,2)-gxa1(im,4)*r2(re,3,2))
375 !rank2c(im,3) = 3.d0*(gxa1(re,5)*r2(im,1,3)-gxa1(im,5)*r2(re,1,3)&
376 !&                     +gxa1(re,4)*r2(im,2,3)-gxa1(im,4)*r2(re,2,3)&
377 !&                     +gxa1(re,3)*r2(im,3,3)-gxa1(im,3)*r2(re,3,3))
378  rank2c(im,4) = 3.d0*(gxa1(re,5)*r2(im,1,2)-gxa1(im,5)*r2(re,1,2)&
379 & +gxa1(re,4)*r2(im,2,2)-gxa1(im,4)*r2(re,2,2)&
380 & +gxa1(re,3)*r2(im,3,2)-gxa1(im,3)*r2(re,3,2))
381  rank2c(im,5) = 3.d0*(gxa1(re,5)*r2(im,1,1)-gxa1(im,5)*r2(re,1,1)&
382 & +gxa1(re,4)*r2(im,2,1)-gxa1(im,4)*r2(re,2,1)&
383 & +gxa1(re,3)*r2(im,3,1)-gxa1(im,3)*r2(re,3,1))
384  rank2c(im,6) = 3.d0*(gxa1(re,6)*r2(im,1,1)-gxa1(im,6)*r2(re,1,1)&
385 & +gxa1(re,2)*r2(im,2,1)-gxa1(im,2)*r2(re,2,1)&
386 & +gxa1(re,4)*r2(im,3,1)-gxa1(im,4)*r2(re,3,1))
387 
388 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

PARENTS

      nonlop_pl

CHILDREN

SOURCE

432 subroutine cont22so(gxa1,gxa2,amet,rank2)
433 
434 
435 !This section has been created automatically by the script Abilint (TD).
436 !Do not modify the following lines by hand.
437 #undef ABI_FUNC
438 #define ABI_FUNC 'cont22so'
439 !End of the abilint section
440 
441  implicit none
442 
443 !Arguments ------------------------------------
444 !arrays
445  real(dp),intent(in) :: amet(2,3,3),gxa1(2,6),gxa2(2,6)
446  real(dp),intent(out) :: rank2(6)
447 
448 !Local variables-------------------------------
449 !scalars
450  integer,parameter :: im=2,re=1
451 
452 ! *************************************************************************
453 
454 !Simply write out index summations
455 !a=1, b=1 in rank2(a,b) --> maps to index 1
456  rank2(1)=2.0d0*(&
457 & amet(1,3,2)*(gxa1(re,5)*gxa2(re,6)+gxa1(im,5)*gxa2(im,6))-&
458 & amet(2,3,2)*(gxa1(re,5)*gxa2(im,6)-gxa1(im,5)*gxa2(re,6))+&
459 & amet(1,3,1)*(gxa1(re,5)*gxa2(re,1)+gxa1(im,5)*gxa2(im,1))-&
460 & amet(2,3,1)*(gxa1(re,5)*gxa2(im,1)-gxa1(im,5)*gxa2(re,1))+&
461 & amet(1,2,1)*(gxa1(re,6)*gxa2(re,1)+gxa1(im,6)*gxa2(im,1))-&
462 & amet(2,2,1)*(gxa1(re,6)*gxa2(im,1)-gxa1(im,6)*gxa2(re,1))+&
463 & amet(1,2,3)*(gxa1(re,6)*gxa2(re,5)+gxa1(im,6)*gxa2(im,5))-&
464 & amet(2,2,3)*(gxa1(re,6)*gxa2(im,5)-gxa1(im,6)*gxa2(re,5))+&
465 & amet(1,1,3)*(gxa1(re,1)*gxa2(re,5)+gxa1(im,1)*gxa2(im,5))-&
466 & amet(2,1,3)*(gxa1(re,1)*gxa2(im,5)-gxa1(im,1)*gxa2(re,5))+&
467 & amet(1,1,2)*(gxa1(re,1)*gxa2(re,6)+gxa1(im,1)*gxa2(im,6))-&
468 & amet(2,1,2)*(gxa1(re,1)*gxa2(im,6)-gxa1(im,1)*gxa2(re,6)))
469 
470 !a=2, b=2 in rank2(a,b) --> maps to index 2
471  rank2(2)=2.0d0*(&
472 & amet(1,3,2)*(gxa1(re,4)*gxa2(re,2)+gxa1(im,4)*gxa2(im,2))-&
473 & amet(2,3,2)*(gxa1(re,4)*gxa2(im,2)-gxa1(im,4)*gxa2(re,2))+&
474 & amet(1,3,1)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-&
475 & amet(2,3,1)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+&
476 & amet(1,2,1)*(gxa1(re,2)*gxa2(re,6)+gxa1(im,2)*gxa2(im,6))-&
477 & amet(2,2,1)*(gxa1(re,2)*gxa2(im,6)-gxa1(im,2)*gxa2(re,6))+&
478 & amet(1,2,3)*(gxa1(re,2)*gxa2(re,4)+gxa1(im,2)*gxa2(im,4))-&
479 & amet(2,2,3)*(gxa1(re,2)*gxa2(im,4)-gxa1(im,2)*gxa2(re,4))+&
480 & amet(1,1,3)*(gxa1(re,6)*gxa2(re,4)+gxa1(im,6)*gxa2(im,4))-&
481 & amet(2,1,3)*(gxa1(re,6)*gxa2(im,4)-gxa1(im,6)*gxa2(re,4))+&
482 & amet(1,1,2)*(gxa1(re,6)*gxa2(re,2)+gxa1(im,6)*gxa2(im,2))-&
483 & amet(2,1,2)*(gxa1(re,6)*gxa2(im,2)-gxa1(im,6)*gxa2(re,2)))
484 
485 !a=3, b=3 in rank2(a,b) --> maps to index 3
486  rank2(3)=2.0d0*(&
487 & amet(1,3,2)*(gxa1(re,3)*gxa2(re,4)+gxa1(im,3)*gxa2(im,4))-&
488 & amet(2,3,2)*(gxa1(re,3)*gxa2(im,4)-gxa1(im,3)*gxa2(re,4))+&
489 & amet(1,3,1)*(gxa1(re,3)*gxa2(re,5)+gxa1(im,3)*gxa2(im,5))-&
490 & amet(2,3,1)*(gxa1(re,3)*gxa2(im,5)-gxa1(im,3)*gxa2(re,5))+&
491 & amet(1,2,1)*(gxa1(re,4)*gxa2(re,5)+gxa1(im,4)*gxa2(im,5))-&
492 & amet(2,2,1)*(gxa1(re,4)*gxa2(im,5)-gxa1(im,4)*gxa2(re,5))+&
493 & amet(1,2,3)*(gxa1(re,4)*gxa2(re,3)+gxa1(im,4)*gxa2(im,3))-&
494 & amet(2,2,3)*(gxa1(re,4)*gxa2(im,3)-gxa1(im,4)*gxa2(re,3))+&
495 & amet(1,1,3)*(gxa1(re,5)*gxa2(re,3)+gxa1(im,5)*gxa2(im,3))-&
496 & amet(2,1,3)*(gxa1(re,5)*gxa2(im,3)-gxa1(im,5)*gxa2(re,3))+&
497 & amet(1,1,2)*(gxa1(re,5)*gxa2(re,4)+gxa1(im,5)*gxa2(im,4))-&
498 & amet(2,1,2)*(gxa1(re,5)*gxa2(im,4)-gxa1(im,5)*gxa2(re,4)))
499 
500 !a=3, b=2 in rank2(a,b) --> maps to index 4
501  rank2(4)=2.0d0*(&
502 & amet(1,3,2)*(gxa1(re,3)*gxa2(re,2)+gxa1(im,3)*gxa2(im,2))-&
503 & amet(2,3,2)*(gxa1(re,3)*gxa2(im,2)-gxa1(im,3)*gxa2(re,2))+&
504 & amet(1,3,1)*(gxa1(re,3)*gxa2(re,6)+gxa1(im,3)*gxa2(im,6))-&
505 & amet(2,3,1)*(gxa1(re,3)*gxa2(im,6)-gxa1(im,3)*gxa2(re,6))+&
506 & amet(1,2,1)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-&
507 & amet(2,2,1)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+&
508 & amet(1,2,3)*(gxa1(re,4)*gxa2(re,4)+gxa1(im,4)*gxa2(im,4))-&
509 & amet(2,2,3)*(gxa1(re,4)*gxa2(im,4)-gxa1(im,4)*gxa2(re,4))+&
510 & amet(1,1,3)*(gxa1(re,5)*gxa2(re,4)+gxa1(im,5)*gxa2(im,4))-&
511 & amet(2,1,3)*(gxa1(re,5)*gxa2(im,4)-gxa1(im,5)*gxa2(re,4))+&
512 & amet(1,1,2)*(gxa1(re,5)*gxa2(re,2)+gxa1(im,5)*gxa2(im,2))-&
513 & amet(2,1,2)*(gxa1(re,5)*gxa2(im,2)-gxa1(im,5)*gxa2(re,2)))
514 
515 !a=3, b=1 in rank2(a,b) --> maps to index 5
516  rank2(5)=2.0d0*(&
517 & amet(1,3,2)*(gxa1(re,3)*gxa2(re,6)+gxa1(im,3)*gxa2(im,6))-&
518 & amet(2,3,2)*(gxa1(re,3)*gxa2(im,6)-gxa1(im,3)*gxa2(re,6))+&
519 & amet(1,3,1)*(gxa1(re,3)*gxa2(re,1)+gxa1(im,3)*gxa2(im,1))-&
520 & amet(2,3,1)*(gxa1(re,3)*gxa2(im,1)-gxa1(im,3)*gxa2(re,1))+&
521 & amet(1,2,1)*(gxa1(re,4)*gxa2(re,1)+gxa1(im,4)*gxa2(im,1))-&
522 & amet(2,2,1)*(gxa1(re,4)*gxa2(im,1)-gxa1(im,4)*gxa2(re,1))+&
523 & amet(1,2,3)*(gxa1(re,4)*gxa2(re,5)+gxa1(im,4)*gxa2(im,5))-&
524 & amet(2,2,3)*(gxa1(re,4)*gxa2(im,5)-gxa1(im,4)*gxa2(re,5))+&
525 & amet(1,1,3)*(gxa1(re,5)*gxa2(re,5)+gxa1(im,5)*gxa2(im,5))-&
526 & amet(2,1,3)*(gxa1(re,5)*gxa2(im,5)-gxa1(im,5)*gxa2(re,5))+&
527 & amet(1,1,2)*(gxa1(re,5)*gxa2(re,6)+gxa1(im,5)*gxa2(im,6))-&
528 & amet(2,1,2)*(gxa1(re,5)*gxa2(im,6)-gxa1(im,5)*gxa2(re,6)))
529 
530 !a=2, b=1 in rank2(a,b) --> maps to index 6
531  rank2(6)=2.0d0*(&
532 & amet(1,3,2)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-&
533 & amet(2,3,2)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+&
534 & amet(1,3,1)*(gxa1(re,4)*gxa2(re,1)+gxa1(im,4)*gxa2(im,1))-&
535 & amet(2,3,1)*(gxa1(re,4)*gxa2(im,1)-gxa1(im,4)*gxa2(re,1))+&
536 & amet(1,2,1)*(gxa1(re,2)*gxa2(re,1)+gxa1(im,2)*gxa2(im,1))-&
537 & amet(2,2,1)*(gxa1(re,2)*gxa2(im,1)-gxa1(im,2)*gxa2(re,1))+&
538 & amet(1,2,3)*(gxa1(re,2)*gxa2(re,5)+gxa1(im,2)*gxa2(im,5))-&
539 & amet(2,2,3)*(gxa1(re,2)*gxa2(im,5)-gxa1(im,2)*gxa2(re,5))+&
540 & amet(1,1,3)*(gxa1(re,6)*gxa2(re,5)+gxa1(im,6)*gxa2(im,5))-&
541 & amet(2,1,3)*(gxa1(re,6)*gxa2(im,5)-gxa1(im,6)*gxa2(re,5))+&
542 & amet(1,1,2)*(gxa1(re,6)*gxa2(re,6)+gxa1(im,6)*gxa2(im,6))-&
543 & amet(2,1,2)*(gxa1(re,6)*gxa2(im,6)-gxa1(im,6)*gxa2(re,6)))
544 
545 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.

PARENTS

      nonlop_pl

CHILDREN

SOURCE

588 subroutine cont24(gxa,rank4,rank2)
589 
590 
591 !This section has been created automatically by the script Abilint (TD).
592 !Do not modify the following lines by hand.
593 #undef ABI_FUNC
594 #define ABI_FUNC 'cont24'
595 !End of the abilint section
596 
597  implicit none
598 
599 !Arguments ------------------------------------
600 !arrays
601  real(dp),intent(in) :: gxa(2,6),rank4(2,15)
602  real(dp),intent(out) :: rank2(6)
603 
604 !Local variables-------------------------------
605 !scalars
606  integer,parameter :: im=2,re=1
607 
608 ! *************************************************************************
609 
610 !Simply write out index summations
611 
612 !a=1, b=1 in rank2(a,b) --> maps to index 1
613  rank2(1)=2.0d0*(&
614 & gxa(re,1)*rank4(re, 1)+gxa(im,1)*rank4(im, 1)+&
615 & gxa(re,2)*rank4(re, 2)+gxa(im,2)*rank4(im, 2)+&
616 & gxa(re,3)*rank4(re, 3)+gxa(im,3)*rank4(im, 3)+&
617 & gxa(re,4)*rank4(re, 4)+gxa(im,4)*rank4(im, 4)+&
618 & gxa(re,5)*rank4(re, 5)+gxa(im,5)*rank4(im, 5)+&
619 & gxa(re,6)*rank4(re, 6)+gxa(im,6)*rank4(im, 6))
620 
621 !a=2, b=2 in rank2(a,b) --> maps to index 2
622  rank2(2)=2.0d0*(&
623 & gxa(re,1)*rank4(re, 2)+gxa(im,1)*rank4(im, 2)+&
624 & gxa(re,2)*rank4(re,11)+gxa(im,2)*rank4(im,11)+&
625 & gxa(re,3)*rank4(re,12)+gxa(im,3)*rank4(im,12)+&
626 & gxa(re,4)*rank4(re,13)+gxa(im,4)*rank4(im,13)+&
627 & gxa(re,5)*rank4(re, 9)+gxa(im,5)*rank4(im, 9)+&
628 & gxa(re,6)*rank4(re, 7)+gxa(im,6)*rank4(im, 7))
629 
630 !a=3, b=3 in rank2(a,b) --> maps to index 3
631  rank2(3)=2.0d0*(&
632 & gxa(re,1)*rank4(re, 3)+gxa(im,1)*rank4(im, 3)+&
633 & gxa(re,2)*rank4(re,12)+gxa(im,2)*rank4(im,12)+&
634 & gxa(re,3)*rank4(re,15)+gxa(im,3)*rank4(im,15)+&
635 & gxa(re,4)*rank4(re,14)+gxa(im,4)*rank4(im,14)+&
636 & gxa(re,5)*rank4(re,10)+gxa(im,5)*rank4(im,10)+&
637 & gxa(re,6)*rank4(re, 8)+gxa(im,6)*rank4(im, 8))
638 
639 !a=3, b=2 in rank2(a,b) --> maps to index 4
640  rank2(4)=2.0d0*(&
641 & gxa(re,1)*rank4(re, 4)+gxa(im,1)*rank4(im, 4)+&
642 & gxa(re,2)*rank4(re,13)+gxa(im,2)*rank4(im,13)+&
643 & gxa(re,3)*rank4(re,14)+gxa(im,3)*rank4(im,14)+&
644 & gxa(re,4)*rank4(re,12)+gxa(im,4)*rank4(im,12)+&
645 & gxa(re,5)*rank4(re, 8)+gxa(im,5)*rank4(im, 8)+&
646 & gxa(re,6)*rank4(re, 9)+gxa(im,6)*rank4(im, 9))
647 
648 !a=3, b=1 in rank2(a,b) --> maps to index 5
649  rank2(5)=2.0d0*(&
650 & gxa(re,1)*rank4(re, 5)+gxa(im,1)*rank4(im, 5)+&
651 & gxa(re,2)*rank4(re, 9)+gxa(im,2)*rank4(im, 9)+&
652 & gxa(re,3)*rank4(re,10)+gxa(im,3)*rank4(im,10)+&
653 & gxa(re,4)*rank4(re, 8)+gxa(im,4)*rank4(im, 8)+&
654 & gxa(re,5)*rank4(re, 3)+gxa(im,5)*rank4(im, 3)+&
655 & gxa(re,6)*rank4(re, 4)+gxa(im,6)*rank4(im, 4))
656 
657 !a=2, b=1 in rank2(a,b) --> maps to index 6
658  rank2(6)=2.0d0*(&
659 & gxa(re,1)*rank4(re, 6)+gxa(im,1)*rank4(im, 6)+&
660 & gxa(re,2)*rank4(re, 7)+gxa(im,2)*rank4(im, 7)+&
661 & gxa(re,3)*rank4(re, 8)+gxa(im,3)*rank4(im, 8)+&
662 & gxa(re,4)*rank4(re, 9)+gxa(im,4)*rank4(im, 9)+&
663 & gxa(re,5)*rank4(re, 4)+gxa(im,5)*rank4(im, 4)+&
664 & gxa(re,6)*rank4(re, 2)+gxa(im,6)*rank4(im, 2))
665 
666 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.

PARENTS

      nonlop_pl

CHILDREN

SOURCE

703 subroutine cont3(gxa,gmet,rank2)
704 
705 
706 !This section has been created automatically by the script Abilint (TD).
707 !Do not modify the following lines by hand.
708 #undef ABI_FUNC
709 #define ABI_FUNC 'cont3'
710 !End of the abilint section
711 
712  implicit none
713 
714 !Arguments ------------------------------------
715 !arrays
716  real(dp),intent(in) :: gmet(3,3),gxa(2,10)
717  real(dp),intent(out) :: rank2(6)
718 
719 !Local variables-------------------------------
720 !scalars
721  integer,parameter :: im=2,re=1
722  integer :: ii
723 !arrays
724  real(dp) :: r1(2,3),r3(2,18),s13(6),s33(6)
725 
726 ! *************************************************************************
727 
728 !Compute r1(a) = gmet(i,j) gxa(i,j,a)
729 
730 !Write out components for 3 distinct terms, Re and Im
731  do ii=1,2
732    r1(ii,1)=gmet(1,1)*gxa(ii,1)+gmet(2,2)*gxa(ii,2)+&
733 &   gmet(3,3)*gxa(ii,3)+2.d0*(&
734 &   gmet(3,2)*gxa(ii,4)+gmet(3,1)*gxa(ii,5)+&
735 &   gmet(2,1)*gxa(ii,6))
736    r1(ii,2)=gmet(1,1)*gxa(ii,6)+gmet(2,2)*gxa(ii,7)+&
737 &   gmet(3,3)*gxa(ii,8)+2.d0*(&
738 &   gmet(3,2)*gxa(ii,9)+gmet(3,1)*gxa(ii,4)+&
739 &   gmet(2,1)*gxa(ii,2))
740    r1(ii,3)=gmet(1,1)*gxa(ii,5)+gmet(2,2)*gxa(ii,9)+&
741 &   gmet(3,3)*gxa(ii,10)+2.d0*(&
742 &   gmet(3,2)*gxa(ii,8)+gmet(3,1)*gxa(ii,3)+&
743 &   gmet(2,1)*gxa(ii,4))
744  end do
745 
746 !Compute r3(a,b,k)=gmet(k,n) gxa(a,b,n)
747 
748 !Write out components for 18 distinct terms, Re and Im
749 !(symmetric in first two indices, not in all permutations)
750 !store as 111 221 331 321 311 211
751 !112 222 332 322 312 212
752 !113 223 333 323 313 213
753  do ii=1,2
754    r3(ii, 1)=gmet(1,1)*gxa(ii,1)+gmet(2,1)*gxa(ii,6)+&
755 &   gmet(3,1)*gxa(ii,5)
756    r3(ii, 2)=gmet(1,1)*gxa(ii,2)+gmet(2,1)*gxa(ii,7)+&
757 &   gmet(3,1)*gxa(ii,9)
758    r3(ii, 3)=gmet(1,1)*gxa(ii,3)+gmet(2,1)*gxa(ii,8)+&
759 &   gmet(3,1)*gxa(ii,10)
760    r3(ii, 4)=gmet(1,1)*gxa(ii,4)+gmet(2,1)*gxa(ii,9)+&
761 &   gmet(3,1)*gxa(ii,8)
762    r3(ii, 5)=gmet(1,1)*gxa(ii,5)+gmet(2,1)*gxa(ii,4)+&
763 &   gmet(3,1)*gxa(ii,3)
764    r3(ii, 6)=gmet(1,1)*gxa(ii,6)+gmet(2,1)*gxa(ii,2)+&
765 &   gmet(3,1)*gxa(ii,4)
766    r3(ii, 7)=gmet(2,1)*gxa(ii,1)+gmet(2,2)*gxa(ii,6)+&
767 &   gmet(3,2)*gxa(ii,5)
768    r3(ii, 8)=gmet(2,1)*gxa(ii,2)+gmet(2,2)*gxa(ii,7)+&
769 &   gmet(3,2)*gxa(ii,9)
770    r3(ii, 9)=gmet(2,1)*gxa(ii,3)+gmet(2,2)*gxa(ii,8)+&
771 &   gmet(3,2)*gxa(ii,10)
772    r3(ii,10)=gmet(2,1)*gxa(ii,4)+gmet(2,2)*gxa(ii,9)+&
773 &   gmet(3,2)*gxa(ii,8)
774    r3(ii,11)=gmet(2,1)*gxa(ii,5)+gmet(2,2)*gxa(ii,4)+&
775 &   gmet(3,2)*gxa(ii,3)
776    r3(ii,12)=gmet(2,1)*gxa(ii,6)+gmet(2,2)*gxa(ii,2)+&
777 &   gmet(3,2)*gxa(ii,4)
778    r3(ii,13)=gmet(3,1)*gxa(ii,1)+gmet(3,2)*gxa(ii,6)+&
779 &   gmet(3,3)*gxa(ii,5)
780    r3(ii,14)=gmet(3,1)*gxa(ii,2)+gmet(3,2)*gxa(ii,7)+&
781 &   gmet(3,3)*gxa(ii,9)
782    r3(ii,15)=gmet(3,1)*gxa(ii,3)+gmet(3,2)*gxa(ii,8)+&
783 &   gmet(3,3)*gxa(ii,10)
784    r3(ii,16)=gmet(3,1)*gxa(ii,4)+gmet(3,2)*gxa(ii,9)+&
785 &   gmet(3,3)*gxa(ii,8)
786    r3(ii,17)=gmet(3,1)*gxa(ii,5)+gmet(3,2)*gxa(ii,4)+&
787 &   gmet(3,3)*gxa(ii,3)
788    r3(ii,18)=gmet(3,1)*gxa(ii,6)+gmet(3,2)*gxa(ii,2)+&
789 &   gmet(3,3)*gxa(ii,4)
790 
791  end do
792 
793 !Now need
794 !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)].
795 
796 !Write out s33(a,b)=2*Re[r3(a,i,j)*r3(b,j,i)]
797 
798  s33(1)=2.d0*(r3(re, 1)*r3(re, 1)+r3(im, 1)*r3(im, 1)+&
799 & r3(re,12)*r3(re,12)+r3(im,12)*r3(im,12)+&
800 & r3(re,17)*r3(re,17)+r3(im,17)*r3(im,17)+&
801 & r3(re,11)*r3(re,18)+r3(im,11)*r3(im,18)+&
802 & r3(re,18)*r3(re,11)+r3(im,18)*r3(im,11)+&
803 & r3(re, 5)*r3(re,13)+r3(im, 5)*r3(im,13)+&
804 & r3(re,13)*r3(re, 5)+r3(im,13)*r3(im, 5)+&
805 & r3(re, 6)*r3(re, 7)+r3(im, 6)*r3(im, 7)+&
806 & r3(re, 7)*r3(re, 6)+r3(im, 7)*r3(im, 6))
807 
808  s33(2)=2.d0*(r3(re, 6)*r3(re, 6)+r3(im, 6)*r3(im, 6)+&
809 & r3(re, 8)*r3(re, 8)+r3(im, 8)*r3(im, 8)+&
810 & r3(re,16)*r3(re,16)+r3(im,16)*r3(im,16)+&
811 & r3(re,10)*r3(re,14)+r3(im,10)*r3(im,14)+&
812 & r3(re,14)*r3(re,10)+r3(im,14)*r3(im,10)+&
813 & r3(re, 4)*r3(re,18)+r3(im, 4)*r3(im,18)+&
814 & r3(re,18)*r3(re, 4)+r3(im,18)*r3(im, 4)+&
815 & r3(re, 2)*r3(re,12)+r3(im, 2)*r3(im,12)+&
816 & r3(re,12)*r3(re, 2)+r3(im,12)*r3(im, 2))
817 
818  s33(3)=2.d0*(r3(re, 5)*r3(re, 5)+r3(im, 5)*r3(im, 5)+&
819 & r3(re,10)*r3(re,10)+r3(im,10)*r3(im,10)+&
820 & r3(re,15)*r3(re,15)+r3(im,15)*r3(im,15)+&
821 & r3(re, 9)*r3(re,16)+r3(im, 9)*r3(im,16)+&
822 & r3(re,16)*r3(re, 9)+r3(im,16)*r3(im, 9)+&
823 & r3(re, 3)*r3(re,17)+r3(im, 3)*r3(im,17)+&
824 & r3(re,17)*r3(re, 3)+r3(im,17)*r3(im, 3)+&
825 & r3(re, 4)*r3(re,11)+r3(im, 4)*r3(im,11)+&
826 & r3(re,11)*r3(re, 4)+r3(im,11)*r3(im, 4))
827 
828  s33(4)=2.d0*(r3(re, 5)*r3(re, 6)+r3(im, 5)*r3(im, 6)+&
829 & r3(re,10)*r3(re, 8)+r3(im,10)*r3(im, 8)+&
830 & r3(re,15)*r3(re,16)+r3(im,15)*r3(im,16)+&
831 & r3(re, 9)*r3(re,14)+r3(im, 9)*r3(im,14)+&
832 & r3(re,16)*r3(re,10)+r3(im,16)*r3(im,10)+&
833 & r3(re, 3)*r3(re,18)+r3(im, 3)*r3(im,18)+&
834 & r3(re,17)*r3(re, 4)+r3(im,17)*r3(im, 4)+&
835 & r3(re, 4)*r3(re,12)+r3(im, 4)*r3(im,12)+&
836 & r3(re,11)*r3(re, 2)+r3(im,11)*r3(im, 2))
837 
838  s33(5)=2.d0*(r3(re, 5)*r3(re, 1)+r3(im, 5)*r3(im, 1)+&
839 & r3(re,10)*r3(re,12)+r3(im,10)*r3(im,12)+&
840 & r3(re,15)*r3(re,17)+r3(im,15)*r3(im,17)+&
841 & r3(re, 9)*r3(re,18)+r3(im, 9)*r3(im,18)+&
842 & r3(re,16)*r3(re,11)+r3(im,16)*r3(im,11)+&
843 & r3(re, 3)*r3(re,13)+r3(im, 3)*r3(im,13)+&
844 & r3(re,17)*r3(re, 5)+r3(im,17)*r3(im, 5)+&
845 & r3(re, 4)*r3(re, 7)+r3(im, 4)*r3(im, 7)+&
846 & r3(re,11)*r3(re, 6)+r3(im,11)*r3(im, 6))
847 
848  s33(6)=2.d0*(r3(re, 6)*r3(re, 1)+r3(im, 6)*r3(im, 1)+&
849 & r3(re, 8)*r3(re,12)+r3(im, 8)*r3(im,12)+&
850 & r3(re,16)*r3(re,17)+r3(im,16)*r3(im,17)+&
851 & r3(re,10)*r3(re,18)+r3(im,10)*r3(im,18)+&
852 & r3(re,14)*r3(re,11)+r3(im,14)*r3(im,11)+&
853 & r3(re, 4)*r3(re,13)+r3(im, 4)*r3(im,13)+&
854 & r3(re,18)*r3(re, 5)+r3(im,18)*r3(im, 5)+&
855 & r3(re, 2)*r3(re, 7)+r3(im, 2)*r3(im, 7)+&
856 & r3(re,12)*r3(re, 6)+r3(im,12)*r3(im, 6))
857 
858 
859 !Write out s13(a,b)=2*Re[r1(i)*r3(a,b,i)]
860 
861  s13(1)=2.d0*(r1(re,1)*r3(re, 1)+r1(im,1)*r3(im, 1)+&
862 & r1(re,2)*r3(re, 7)+r1(im,2)*r3(im, 7)+&
863 & r1(re,3)*r3(re,13)+r1(im,3)*r3(im,13))
864  s13(2)=2.d0*(r1(re,1)*r3(re, 2)+r1(im,1)*r3(im, 2)+&
865 & r1(re,2)*r3(re, 8)+r1(im,2)*r3(im, 8)+&
866 & r1(re,3)*r3(re,14)+r1(im,3)*r3(im,14))
867  s13(3)=2.d0*(r1(re,1)*r3(re, 3)+r1(im,1)*r3(im, 3)+&
868 & r1(re,2)*r3(re, 9)+r1(im,2)*r3(im, 9)+&
869 & r1(re,3)*r3(re,15)+r1(im,3)*r3(im,15))
870  s13(4)=2.d0*(r1(re,1)*r3(re, 4)+r1(im,1)*r3(im, 4)+&
871 & r1(re,2)*r3(re,10)+r1(im,2)*r3(im,10)+&
872 & r1(re,3)*r3(re,16)+r1(im,3)*r3(im,16))
873  s13(5)=2.d0*(r1(re,1)*r3(re, 5)+r1(im,1)*r3(im, 5)+&
874 & r1(re,2)*r3(re,11)+r1(im,2)*r3(im,11)+&
875 & r1(re,3)*r3(re,17)+r1(im,3)*r3(im,17))
876  s13(6)=2.d0*(r1(re,1)*r3(re, 6)+r1(im,1)*r3(im, 6)+&
877 & r1(re,2)*r3(re,12)+r1(im,2)*r3(im,12)+&
878 & r1(re,3)*r3(re,18)+r1(im,3)*r3(im,18))
879 
880 !Finally, write out the six terms as final answer
881 !rank2(a,b)=(15/2)*s33(a,b)-3*s13(a,b)-(3/2)*2*Re[r1(a)*r1(b)]
882 
883  rank2(1)=7.5d0*s33(1)-3.d0*s13(1)&
884 & -3.d0*(r1(re,1)*r1(re,1)+r1(im,1)*r1(im,1))
885  rank2(2)=7.5d0*s33(2)-3.d0*s13(2)&
886 & -3.d0*(r1(re,2)*r1(re,2)+r1(im,2)*r1(im,2))
887  rank2(3)=7.5d0*s33(3)-3.d0*s13(3)&
888 & -3.d0*(r1(re,3)*r1(re,3)+r1(im,3)*r1(im,3))
889  rank2(4)=7.5d0*s33(4)-3.d0*s13(4)&
890 & -3.d0*(r1(re,3)*r1(re,2)+r1(im,3)*r1(im,2))
891  rank2(5)=7.5d0*s33(5)-3.d0*s13(5)&
892 & -3.d0*(r1(re,3)*r1(re,1)+r1(im,3)*r1(im,1))
893  rank2(6)=7.5d0*s33(6)-3.d0*s13(6)&
894 & -3.d0*(r1(re,2)*r1(re,1)+r1(im,2)*r1(im,1))
895 
896 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} }}

PARENTS

      nonlop_pl

CHILDREN

SOURCE

 946 subroutine cont33cso(gxa1,gxa2,gmet,rank2c)
 947 
 948 
 949 !This section has been created automatically by the script Abilint (TD).
 950 !Do not modify the following lines by hand.
 951 #undef ABI_FUNC
 952 #define ABI_FUNC 'cont33cso'
 953 !End of the abilint section
 954 
 955  implicit none
 956 
 957 !Arguments ------------------------------------
 958 !arrays
 959  real(dp),intent(in) :: gmet(3,3),gxa1(2,10),gxa2(2,10)
 960  real(dp),intent(out) :: rank2c(2,6)
 961 
 962 !Local variables-------------------------------
 963 !scalars
 964  integer,parameter :: im=2,re=1
 965 !arrays
 966  real(dp) :: r11(2,3),r12(2,3),r3(2,6,3),r3a(2,6,3)
 967 
 968 ! *************************************************************************
 969 
 970 !Initialize output tensor
 971  rank2c(:,:)=0.d0
 972 
 973 !Compute r3(i,j,b)=gxa2(b,l,m)*gmet(i,l)*gmet(j,m)
 974 !stored as 11 22 33 32 31 21 for (i,j)
 975 !First compute r3a(b,l,j)=gxa2(b,l,m)*gmet(j,m)
 976 !stored as r3a(re/im,(b,l),j)
 977  r3a(:,1,1)=gxa2(:,1 )*gmet(1,1)+gxa2(:,6 )*gmet(1,2)+gxa2(:,5 )*gmet(1,3)
 978  r3a(:,2,1)=gxa2(:,2 )*gmet(1,1)+gxa2(:,7 )*gmet(1,2)+gxa2(:,9 )*gmet(1,3)
 979  r3a(:,3,1)=gxa2(:,3 )*gmet(1,1)+gxa2(:,8 )*gmet(1,2)+gxa2(:,10)*gmet(1,3)
 980  r3a(:,4,1)=gxa2(:,4 )*gmet(1,1)+gxa2(:,9 )*gmet(1,2)+gxa2(:,8 )*gmet(1,3)
 981  r3a(:,5,1)=gxa2(:,5 )*gmet(1,1)+gxa2(:,4 )*gmet(1,2)+gxa2(:,3 )*gmet(1,3)
 982  r3a(:,6,1)=gxa2(:,6 )*gmet(1,1)+gxa2(:,2 )*gmet(1,2)+gxa2(:,4 )*gmet(1,3)
 983  r3a(:,1,2)=gxa2(:,1 )*gmet(2,1)+gxa2(:,6 )*gmet(2,2)+gxa2(:,5 )*gmet(2,3)
 984  r3a(:,2,2)=gxa2(:,2 )*gmet(2,1)+gxa2(:,7 )*gmet(2,2)+gxa2(:,9 )*gmet(2,3)
 985  r3a(:,3,2)=gxa2(:,3 )*gmet(2,1)+gxa2(:,8 )*gmet(2,2)+gxa2(:,10)*gmet(2,3)
 986  r3a(:,4,2)=gxa2(:,4 )*gmet(2,1)+gxa2(:,9 )*gmet(2,2)+gxa2(:,8 )*gmet(2,3)
 987  r3a(:,5,2)=gxa2(:,5 )*gmet(2,1)+gxa2(:,4 )*gmet(2,2)+gxa2(:,3 )*gmet(2,3)
 988  r3a(:,6,2)=gxa2(:,6 )*gmet(2,1)+gxa2(:,2 )*gmet(2,2)+gxa2(:,4 )*gmet(2,3)
 989  r3a(:,1,3)=gxa2(:,1 )*gmet(3,1)+gxa2(:,6 )*gmet(3,2)+gxa2(:,5 )*gmet(3,3)
 990  r3a(:,2,3)=gxa2(:,2 )*gmet(3,1)+gxa2(:,7 )*gmet(3,2)+gxa2(:,9 )*gmet(3,3)
 991  r3a(:,3,3)=gxa2(:,3 )*gmet(3,1)+gxa2(:,8 )*gmet(3,2)+gxa2(:,10)*gmet(3,3)
 992  r3a(:,4,3)=gxa2(:,4 )*gmet(3,1)+gxa2(:,9 )*gmet(3,2)+gxa2(:,8 )*gmet(3,3)
 993  r3a(:,5,3)=gxa2(:,5 )*gmet(3,1)+gxa2(:,4 )*gmet(3,2)+gxa2(:,3 )*gmet(3,3)
 994  r3a(:,6,3)=gxa2(:,6 )*gmet(3,1)+gxa2(:,2 )*gmet(3,2)+gxa2(:,4 )*gmet(3,3)
 995 
 996 !Then  compute r3(i,j,b)=r3a(b,l,j)*gmet(i,l)
 997 !stored as r3(re/im,(i,j),b)
 998  r3(:,1,1)=r3a(:,1,1)*gmet(1,1)+r3a(:,6,1)*gmet(1,2)+r3a(:,5,1)*gmet(1,3)
 999  r3(:,1,2)=r3a(:,6,1)*gmet(1,1)+r3a(:,2,1)*gmet(1,2)+r3a(:,4,1)*gmet(1,3)
1000  r3(:,1,3)=r3a(:,5,1)*gmet(1,1)+r3a(:,4,1)*gmet(1,2)+r3a(:,3,1)*gmet(1,3)
1001  r3(:,2,1)=r3a(:,1,2)*gmet(2,1)+r3a(:,6,2)*gmet(2,2)+r3a(:,5,2)*gmet(2,3)
1002  r3(:,2,2)=r3a(:,6,2)*gmet(2,1)+r3a(:,2,2)*gmet(2,2)+r3a(:,4,2)*gmet(2,3)
1003  r3(:,2,3)=r3a(:,5,2)*gmet(2,1)+r3a(:,4,2)*gmet(2,2)+r3a(:,3,2)*gmet(2,3)
1004  r3(:,3,1)=r3a(:,1,3)*gmet(3,1)+r3a(:,6,3)*gmet(3,2)+r3a(:,5,3)*gmet(3,3)
1005  r3(:,3,2)=r3a(:,6,3)*gmet(3,1)+r3a(:,2,3)*gmet(3,2)+r3a(:,4,3)*gmet(3,3)
1006  r3(:,3,3)=r3a(:,5,3)*gmet(3,1)+r3a(:,4,3)*gmet(3,2)+r3a(:,3,3)*gmet(3,3)
1007  r3(:,4,1)=r3a(:,1,2)*gmet(3,1)+r3a(:,6,2)*gmet(3,2)+r3a(:,5,2)*gmet(3,3)
1008  r3(:,4,2)=r3a(:,6,2)*gmet(3,1)+r3a(:,2,2)*gmet(3,2)+r3a(:,4,2)*gmet(3,3)
1009  r3(:,4,3)=r3a(:,5,2)*gmet(3,1)+r3a(:,4,2)*gmet(3,2)+r3a(:,3,2)*gmet(3,3)
1010  r3(:,5,1)=r3a(:,1,1)*gmet(3,1)+r3a(:,6,1)*gmet(3,2)+r3a(:,5,1)*gmet(3,3)
1011  r3(:,5,2)=r3a(:,6,1)*gmet(3,1)+r3a(:,2,1)*gmet(3,2)+r3a(:,4,1)*gmet(3,3)
1012  r3(:,5,3)=r3a(:,5,1)*gmet(3,1)+r3a(:,4,1)*gmet(3,2)+r3a(:,3,1)*gmet(3,3)
1013  r3(:,6,1)=r3a(:,1,1)*gmet(2,1)+r3a(:,6,1)*gmet(2,2)+r3a(:,5,1)*gmet(2,3)
1014  r3(:,6,2)=r3a(:,6,1)*gmet(2,1)+r3a(:,2,1)*gmet(2,2)+r3a(:,4,1)*gmet(2,3)
1015  r3(:,6,3)=r3a(:,5,1)*gmet(2,1)+r3a(:,4,1)*gmet(2,2)+r3a(:,3,1)*gmet(2,3)
1016 
1017 !Compute r11(a)=conjg(gxa1(a,l,m))*gmet(l,m)
1018  r11(:,1)=gxa1(:,1)*gmet(1,1)+gxa1(:,2)*gmet(2,2)+gxa1(:,3 )*gmet(3,3)&
1019 & +2.d0*(gxa1(:,4)*gmet(3,2)+gxa1(:,5)*gmet(3,1)+gxa1(:,6 )*gmet(2,1))
1020  r11(:,2)=gxa1(:,6)*gmet(1,1)+gxa1(:,7)*gmet(2,2)+gxa1(:,8 )*gmet(3,3)&
1021 & +2.d0*(gxa1(:,9)*gmet(3,2)+gxa1(:,4)*gmet(3,1)+gxa1(:,2 )*gmet(2,1))
1022  r11(:,3)=gxa1(:,5)*gmet(1,1)+gxa1(:,9)*gmet(2,2)+gxa1(:,10)*gmet(3,3)&
1023 & +2.d0*(gxa1(:,8)*gmet(3,2)+gxa1(:,3)*gmet(3,1)+gxa1(:,4 )*gmet(2,1))
1024  r11(im,1)=-r11(im,1);r11(im,2)=-r11(im,2);r11(im,3)=-r11(im,3)
1025 
1026 !Compute r12(b)=gxa2(b,l,m)*gmet(l,m)
1027  r12(:,1)=gxa2(:,1)*gmet(1,1)+gxa2(:,2)*gmet(2,2)+gxa2(:,3 )*gmet(3,3)&
1028 & +2.d0*(gxa2(:,4)*gmet(3,2)+gxa2(:,5)*gmet(3,1)+gxa2(:,6 )*gmet(2,1))
1029  r12(:,2)=gxa2(:,6)*gmet(1,1)+gxa2(:,7)*gmet(2,2)+gxa2(:,8 )*gmet(3,3)&
1030 & +2.d0*(gxa2(:,9)*gmet(3,2)+gxa2(:,4)*gmet(3,1)+gxa2(:,2 )*gmet(2,1))
1031  r12(:,3)=gxa2(:,5)*gmet(1,1)+gxa2(:,9)*gmet(2,2)+gxa2(:,10)*gmet(3,3)&
1032 & +2.d0*(gxa2(:,8)*gmet(3,2)+gxa2(:,3)*gmet(3,1)+gxa2(:,4 )*gmet(2,1))
1033 
1034 !Finally compute rank2c(a,b)=7.5*conjg(gxa1(a,i,j))*r3(i,j,b) - 1.5*r11(a)*r12(b)
1035 !rank2c(re,1)=7.5d0*(gxa1(re,1 )*r3(re,1,1)+gxa1(im,1 )*r3(im,1,1)&
1036 !&                    +gxa1(re,2 )*r3(re,2,1)+gxa1(im,2 )*r3(im,2,1)&
1037 !&                    +gxa1(re,3 )*r3(re,3,1)+gxa1(im,3 )*r3(im,3,1))&
1038 !&             +15.d0*(gxa1(re,4 )*r3(re,4,1)+gxa1(im,4 )*r3(im,4,1)&
1039 !&                    +gxa1(re,5 )*r3(re,5,1)+gxa1(im,5 )*r3(im,5,1)&
1040 !&                    +gxa1(re,6 )*r3(re,6,1)+gxa1(im,6 )*r3(im,6,1))&
1041 !&             -1.5d0*(r11(re,1)*r12(re,1)-r11(im,1)*r12(im,1))
1042 !rank2c(re,2)=7.5d0*(gxa1(re,6 )*r3(re,1,2)+gxa1(im,6 )*r3(im,1,2)&
1043 !&                    +gxa1(re,7 )*r3(re,2,2)+gxa1(im,7 )*r3(im,2,2)&
1044 !&                    +gxa1(re,8 )*r3(re,3,2)+gxa1(im,8 )*r3(im,3,2))&
1045 !&             +15.d0*(gxa1(re,9 )*r3(re,4,2)+gxa1(im,9 )*r3(im,4,2)&
1046 !&                    +gxa1(re,4 )*r3(re,5,2)+gxa1(im,4 )*r3(im,5,2)&
1047 !&                    +gxa1(re,2 )*r3(re,6,2)+gxa1(im,2 )*r3(im,6,2))&
1048 !&             -1.5d0*(r11(re,2)*r12(re,2)-r11(im,2)*r12(im,2))
1049 !rank2c(re,3)=7.5d0*(gxa1(re,5 )*r3(re,1,3)+gxa1(im,5 )*r3(im,1,3)&
1050 !&                    +gxa1(re,9 )*r3(re,2,3)+gxa1(im,9 )*r3(im,2,3)&
1051 !&                    +gxa1(re,10)*r3(re,3,3)+gxa1(im,10)*r3(im,3,3))&
1052 !&             +15.d0*(gxa1(re,8 )*r3(re,4,3)+gxa1(im,8 )*r3(im,4,3)&
1053 !&                    +gxa1(re,3 )*r3(re,5,3)+gxa1(im,3 )*r3(im,5,3)&
1054 !&                    +gxa1(re,4 )*r3(re,6,3)+gxa1(im,4 )*r3(im,6,3))&
1055 !&             -1.5d0*(r11(re,3)*r12(re,3)-r11(im,3)*r12(im,3))
1056  rank2c(re,4)=7.5d0*(gxa1(re,5 )*r3(re,1,2)+gxa1(im,5 )*r3(im,1,2)&
1057 & +gxa1(re,9 )*r3(re,2,2)+gxa1(im,9 )*r3(im,2,2)&
1058 & +gxa1(re,10)*r3(re,3,2)+gxa1(im,10)*r3(im,3,2))&
1059 & +15.d0*(gxa1(re,8 )*r3(re,4,2)+gxa1(im,8 )*r3(im,4,2)&
1060 & +gxa1(re,3 )*r3(re,5,2)+gxa1(im,3 )*r3(im,5,2)&
1061 & +gxa1(re,4 )*r3(re,6,2)+gxa1(im,4 )*r3(im,6,2))&
1062 & -1.5d0*(r11(re,3)*r12(re,2)-r11(im,3)*r12(im,2))
1063  rank2c(re,5)=7.5d0*(gxa1(re,5 )*r3(re,1,1)+gxa1(im,5 )*r3(im,1,1)&
1064 & +gxa1(re,9 )*r3(re,2,1)+gxa1(im,9 )*r3(im,2,1)&
1065 & +gxa1(re,10)*r3(re,3,1)+gxa1(im,10)*r3(im,3,1))&
1066 & +15.d0*(gxa1(re,8 )*r3(re,4,1)+gxa1(im,8 )*r3(im,4,1)&
1067 & +gxa1(re,3 )*r3(re,5,1)+gxa1(im,3 )*r3(im,5,1)&
1068 & +gxa1(re,4 )*r3(re,6,1)+gxa1(im,4 )*r3(im,6,1))&
1069 & -1.5d0*(r11(re,3)*r12(re,1)-r11(im,3)*r12(im,1))
1070  rank2c(re,6)=7.5d0*(gxa1(re,6 )*r3(re,1,1)+gxa1(im,6 )*r3(im,1,1)&
1071 & +gxa1(re,7 )*r3(re,2,1)+gxa1(im,7 )*r3(im,2,1)&
1072 & +gxa1(re,8 )*r3(re,3,1)+gxa1(im,8 )*r3(im,3,1))&
1073 & +15.d0*(gxa1(re,9 )*r3(re,4,1)+gxa1(im,9 )*r3(im,4,1)&
1074 & +gxa1(re,4 )*r3(re,5,1)+gxa1(im,4 )*r3(im,5,1)&
1075 & +gxa1(re,2 )*r3(re,6,1)+gxa1(im,2 )*r3(im,6,1))&
1076 & -1.5d0*(r11(re,2)*r12(re,1)-r11(im,2)*r12(im,1))
1077 !rank2c(im,1)=7.5d0*(gxa1(re,1 )*r3(im,1,1)-gxa1(im,1 )*r3(re,1,1)&
1078 !&                    +gxa1(re,2 )*r3(im,2,1)-gxa1(im,2 )*r3(re,2,1)&
1079 !&                    +gxa1(re,3 )*r3(im,3,1)-gxa1(im,3 )*r3(re,3,1))&
1080 !&             +15.d0*(gxa1(re,4 )*r3(im,4,1)-gxa1(im,4 )*r3(re,4,1)&
1081 !&                    +gxa1(re,5 )*r3(im,5,1)-gxa1(im,5 )*r3(re,5,1)&
1082 !&                    +gxa1(re,6 )*r3(im,6,1)-gxa1(im,6 )*r3(re,6,1))&
1083 !&             -1.5d0*(r11(re,1)*r12(im,1)+r11(im,1)*r12(re,1))
1084 !rank2c(im,2)=7.5d0*(gxa1(re,6 )*r3(im,1,2)-gxa1(im,6 )*r3(re,1,2)&
1085 !&                    +gxa1(re,7 )*r3(im,2,2)-gxa1(im,7 )*r3(re,2,2)&
1086 !&                    +gxa1(re,8 )*r3(im,3,2)-gxa1(im,8 )*r3(re,3,2))&
1087 !&             +15.d0*(gxa1(re,9 )*r3(im,4,2)-gxa1(im,9 )*r3(re,4,2)&
1088 !&                    +gxa1(re,4 )*r3(im,5,2)-gxa1(im,4 )*r3(re,5,2)&
1089 !&                    +gxa1(re,2 )*r3(im,6,2)-gxa1(im,2 )*r3(re,6,2))&
1090 !&             -1.5d0*(r11(re,2)*r12(im,2)+r11(im,2)*r12(re,2))
1091 !rank2c(im,3)=7.5d0*(gxa1(re,5 )*r3(im,1,3)-gxa1(im,5 )*r3(re,1,3)&
1092 !&                    +gxa1(re,9 )*r3(im,2,3)-gxa1(im,9 )*r3(re,2,3)&
1093 !&                    +gxa1(re,10)*r3(im,3,3)-gxa1(im,10)*r3(re,3,3))&
1094 !&             +15.d0*(gxa1(re,8 )*r3(im,4,3)-gxa1(im,8 )*r3(re,4,3)&
1095 !&                    +gxa1(re,3 )*r3(im,5,3)-gxa1(im,3 )*r3(re,5,3)&
1096 !&                    +gxa1(re,4 )*r3(im,6,3)-gxa1(im,4 )*r3(re,6,3))&
1097 !&             -1.5d0*(r11(re,3)*r12(im,3)+r11(im,3)*r12(re,3))
1098  rank2c(im,4)=7.5d0*(gxa1(re,5 )*r3(im,1,2)-gxa1(im,5 )*r3(re,1,2)&
1099 & +gxa1(re,9 )*r3(im,2,2)-gxa1(im,9 )*r3(re,2,2)&
1100 & +gxa1(re,10)*r3(im,3,2)-gxa1(im,10)*r3(re,3,2))&
1101 & +15.d0*(gxa1(re,8 )*r3(im,4,2)-gxa1(im,8 )*r3(re,4,2)&
1102 & +gxa1(re,3 )*r3(im,5,2)-gxa1(im,3 )*r3(re,5,2)&
1103 & +gxa1(re,4 )*r3(im,6,2)-gxa1(im,4 )*r3(re,6,2))&
1104 & -1.5d0*(r11(re,3)*r12(im,2)+r11(im,3)*r12(re,2))
1105  rank2c(im,5)=7.5d0*(gxa1(re,5 )*r3(im,1,1)-gxa1(im,5 )*r3(re,1,1)&
1106 & +gxa1(re,9 )*r3(im,2,1)-gxa1(im,9 )*r3(re,2,1)&
1107 & +gxa1(re,10)*r3(im,3,1)-gxa1(im,10)*r3(re,3,1))&
1108 & +15.d0*(gxa1(re,8 )*r3(im,4,1)-gxa1(im,8 )*r3(re,4,1)&
1109 & +gxa1(re,3 )*r3(im,5,1)-gxa1(im,3 )*r3(re,5,1)&
1110 & +gxa1(re,4 )*r3(im,6,1)-gxa1(im,4 )*r3(re,6,1))&
1111 & -1.5d0*(r11(re,3)*r12(im,1)+r11(im,3)*r12(re,1))
1112  rank2c(im,6)=7.5d0*(gxa1(re,6 )*r3(im,1,1)-gxa1(im,6 )*r3(re,1,1)&
1113 & +gxa1(re,7 )*r3(im,2,1)-gxa1(im,7 )*r3(re,2,1)&
1114 & +gxa1(re,8 )*r3(im,3,1)-gxa1(im,8 )*r3(re,3,1))&
1115 & +15.d0*(gxa1(re,9 )*r3(im,4,1)-gxa1(im,9 )*r3(re,4,1)&
1116 & +gxa1(re,4 )*r3(im,5,1)-gxa1(im,4 )*r3(re,5,1)&
1117 & +gxa1(re,2 )*r3(im,6,1)-gxa1(im,2 )*r3(re,6,1))&
1118 & -1.5d0*(r11(re,2)*r12(im,1)+r11(im,2)*r12(re,1))
1119 
1120 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} }}

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1172 subroutine cont33so(gxa1,gxa2,gmet,amet,rank2)
1173 
1174 
1175 !This section has been created automatically by the script Abilint (TD).
1176 !Do not modify the following lines by hand.
1177 #undef ABI_FUNC
1178 #define ABI_FUNC 'cont33so'
1179 !End of the abilint section
1180 
1181  implicit none
1182 
1183 !Arguments ------------------------------------
1184 !arrays
1185  real(dp),intent(in) :: amet(2,3,3),gmet(3,3),gxa1(2,10),gxa2(2,10)
1186  real(dp),intent(out) :: rank2(6)
1187 
1188 !Local variables-------------------------------
1189 !scalars
1190  integer,parameter :: im=2,re=1
1191  integer :: ii
1192 !arrays
1193  real(dp) :: r1(2,3),r3A(2,18),r3G(2,18),s31(6),s33(6)
1194 
1195 ! *************************************************************************
1196 
1197 !Compute r1(i)=gxa2(i,j,k)*gmet(j,k)
1198 
1199  do ii=1,2
1200    r1(ii,1)=gmet(1,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,2)+&
1201 &   gmet(3,3)*gxa2(ii,3)+2.d0*(&
1202 &   gmet(3,2)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,5)+&
1203 &   gmet(2,1)*gxa2(ii,6))
1204    r1(ii,2)=gmet(1,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,7)+&
1205 &   gmet(3,3)*gxa2(ii,8)+2.d0*(&
1206 &   gmet(3,2)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,4)+&
1207 &   gmet(2,1)*gxa2(ii,2))
1208    r1(ii,3)=gmet(1,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,9)+&
1209 &   gmet(3,3)*gxa2(ii,10)+2.d0*(&
1210 &   gmet(3,2)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,3)+&
1211 &   gmet(2,1)*gxa2(ii,4))
1212  end do
1213 
1214 !Compute r3G(i,j,k)=gxa2(p,i,j)*gmet(p,k)
1215 !Write out components for 18 distinct terms, Re and Im
1216 !(r3G is symmetric in first two indices, not in all permutations)
1217 !Store as 111 221 331 321 311 211
1218 !112 222 332 322 312 212
1219 !113 223 333 323 313 213
1220  do ii=1,2
1221    r3G(ii, 1)=gmet(1,1)*gxa2(ii,1)+gmet(2,1)*gxa2(ii,6)+gmet(3,1)*gxa2(ii,5)
1222    r3G(ii, 2)=gmet(1,1)*gxa2(ii,2)+gmet(2,1)*gxa2(ii,7)+gmet(3,1)*gxa2(ii,9)
1223    r3G(ii, 3)=gmet(1,1)*gxa2(ii,3)+gmet(2,1)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,10)
1224    r3G(ii, 4)=gmet(1,1)*gxa2(ii,4)+gmet(2,1)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,8)
1225    r3G(ii, 5)=gmet(1,1)*gxa2(ii,5)+gmet(2,1)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,3)
1226    r3G(ii, 6)=gmet(1,1)*gxa2(ii,6)+gmet(2,1)*gxa2(ii,2)+gmet(3,1)*gxa2(ii,4)
1227    r3G(ii, 7)=gmet(2,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,5)
1228    r3G(ii, 8)=gmet(2,1)*gxa2(ii,2)+gmet(2,2)*gxa2(ii,7)+gmet(3,2)*gxa2(ii,9)
1229    r3G(ii, 9)=gmet(2,1)*gxa2(ii,3)+gmet(2,2)*gxa2(ii,8)+gmet(3,2)*gxa2(ii,10)
1230    r3G(ii,10)=gmet(2,1)*gxa2(ii,4)+gmet(2,2)*gxa2(ii,9)+gmet(3,2)*gxa2(ii,8)
1231    r3G(ii,11)=gmet(2,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,3)
1232    r3G(ii,12)=gmet(2,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,4)
1233    r3G(ii,13)=gmet(3,1)*gxa2(ii,1)+gmet(3,2)*gxa2(ii,6)+gmet(3,3)*gxa2(ii,5)
1234    r3G(ii,14)=gmet(3,1)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,7)+gmet(3,3)*gxa2(ii,9)
1235    r3G(ii,15)=gmet(3,1)*gxa2(ii,3)+gmet(3,2)*gxa2(ii,8)+gmet(3,3)*gxa2(ii,10)
1236    r3G(ii,16)=gmet(3,1)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,9)+gmet(3,3)*gxa2(ii,8)
1237    r3G(ii,17)=gmet(3,1)*gxa2(ii,5)+gmet(3,2)*gxa2(ii,4)+gmet(3,3)*gxa2(ii,3)
1238    r3G(ii,18)=gmet(3,1)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,2)+gmet(3,3)*gxa2(ii,4)
1239  end do
1240 
1241 !Compute r3A(i,j,k)=conjg(gxa1(p,i,j))*amet(p,k)
1242 !Write out components for 18 distinct terms, Re and Im
1243 !(r3A is symmetric in first two indices, not in all permutations)
1244 !Store as 111 221 331 321 311 211
1245 !112 222 332 322 312 212
1246 !113 223 333 323 313 213
1247 !Note that, since amet is antisymmetric, amet(i,i)=0
1248 
1249  r3A(re, 1)=amet(re,2,1)*gxa1(re,6 )+amet(im,2,1)*gxa1(im,6 )+&
1250 & amet(re,3,1)*gxa1(re,5 )+amet(im,3,1)*gxa1(im,5 )
1251  r3A(re, 2)=amet(re,2,1)*gxa1(re,7 )+amet(im,2,1)*gxa1(im,7 )+&
1252 & amet(re,3,1)*gxa1(re,9 )+amet(im,3,1)*gxa1(im,9 )
1253  r3A(re, 3)=amet(re,2,1)*gxa1(re,8 )+amet(im,2,1)*gxa1(im,8 )+&
1254 & amet(re,3,1)*gxa1(re,10)+amet(im,3,1)*gxa1(im,10)
1255  r3A(re, 4)=amet(re,2,1)*gxa1(re,9 )+amet(im,2,1)*gxa1(im,9 )+&
1256 & amet(re,3,1)*gxa1(re,8 )+amet(im,3,1)*gxa1(im,8 )
1257  r3A(re, 5)=amet(re,2,1)*gxa1(re,4 )+amet(im,2,1)*gxa1(im,4 )+&
1258 & amet(re,3,1)*gxa1(re,3 )+amet(im,3,1)*gxa1(im,3 )
1259  r3A(re, 6)=amet(re,2,1)*gxa1(re,2 )+amet(im,2,1)*gxa1(im,2 )+&
1260 & amet(re,3,1)*gxa1(re,4 )+amet(im,3,1)*gxa1(im,4 )
1261  r3A(re, 7)=amet(re,1,2)*gxa1(re,1 )+amet(im,1,2)*gxa1(im,1 )+&
1262 & amet(re,3,2)*gxa1(re,5 )+amet(im,3,2)*gxa1(im,5 )
1263  r3A(re, 8)=amet(re,1,2)*gxa1(re,2 )+amet(im,1,2)*gxa1(im,2 )+&
1264 & amet(re,3,2)*gxa1(re,9 )+amet(im,3,2)*gxa1(im,9 )
1265  r3A(re, 9)=amet(re,1,2)*gxa1(re,3 )+amet(im,1,2)*gxa1(im,3 )+&
1266 & amet(re,3,2)*gxa1(re,10)+amet(im,3,2)*gxa1(im,10)
1267  r3A(re,10)=amet(re,1,2)*gxa1(re,4 )+amet(im,1,2)*gxa1(im,4 )+&
1268 & amet(re,3,2)*gxa1(re,8 )+amet(im,3,2)*gxa1(im,8 )
1269  r3A(re,11)=amet(re,1,2)*gxa1(re,5 )+amet(im,1,2)*gxa1(im,5 )+&
1270 & amet(re,3,2)*gxa1(re,3 )+amet(im,3,2)*gxa1(im,3 )
1271  r3A(re,12)=amet(re,1,2)*gxa1(re,6 )+amet(im,1,2)*gxa1(im,6 )+&
1272 & amet(re,3,2)*gxa1(re,4 )+amet(im,3,2)*gxa1(im,4 )
1273  r3A(re,13)=amet(re,1,3)*gxa1(re,1 )+amet(im,1,3)*gxa1(im,1 )+&
1274 & amet(re,2,3)*gxa1(re,6 )+amet(im,2,3)*gxa1(im,6 )
1275  r3A(re,14)=amet(re,1,3)*gxa1(re,2 )+amet(im,1,3)*gxa1(im,2 )+&
1276 & amet(re,2,3)*gxa1(re,7 )+amet(im,2,3)*gxa1(im,7 )
1277  r3A(re,15)=amet(re,1,3)*gxa1(re,3 )+amet(im,1,3)*gxa1(im,3 )+&
1278 & amet(re,2,3)*gxa1(re,8 )+amet(im,2,3)*gxa1(im,8 )
1279  r3A(re,16)=amet(re,1,3)*gxa1(re,4 )+amet(im,1,3)*gxa1(im,4 )+&
1280 & amet(re,2,3)*gxa1(re,9 )+amet(im,2,3)*gxa1(im,9 )
1281  r3A(re,17)=amet(re,1,3)*gxa1(re,5 )+amet(im,1,3)*gxa1(im,5 )+&
1282 & amet(re,2,3)*gxa1(re,4 )+amet(im,2,3)*gxa1(im,4 )
1283  r3A(re,18)=amet(re,1,3)*gxa1(re,6 )+amet(im,1,3)*gxa1(im,6 )+&
1284 & amet(re,2,3)*gxa1(re,2 )+amet(im,2,3)*gxa1(im,2 )
1285 
1286  r3A(im, 1)=amet(im,2,1)*gxa1(re,6 )-amet(re,2,1)*gxa1(im,6 )+&
1287 & amet(im,3,1)*gxa1(re,5 )-amet(re,3,1)*gxa1(im,5 )
1288  r3A(im, 2)=amet(im,2,1)*gxa1(re,7 )-amet(re,2,1)*gxa1(im,7 )+&
1289 & amet(im,3,1)*gxa1(re,9 )-amet(re,3,1)*gxa1(im,9 )
1290  r3A(im, 3)=amet(im,2,1)*gxa1(re,8 )-amet(re,2,1)*gxa1(im,8 )+&
1291 & amet(im,3,1)*gxa1(re,10)-amet(re,3,1)*gxa1(im,10)
1292  r3A(im, 4)=amet(im,2,1)*gxa1(re,9 )-amet(re,2,1)*gxa1(im,9 )+&
1293 & amet(im,3,1)*gxa1(re,8 )-amet(re,3,1)*gxa1(im,8 )
1294  r3A(im, 5)=amet(im,2,1)*gxa1(re,4 )-amet(re,2,1)*gxa1(im,4 )+&
1295 & amet(im,3,1)*gxa1(re,3 )-amet(re,3,1)*gxa1(im,3 )
1296  r3A(im, 6)=amet(im,2,1)*gxa1(re,2 )-amet(re,2,1)*gxa1(im,2 )+&
1297 & amet(im,3,1)*gxa1(re,4 )-amet(re,3,1)*gxa1(im,4 )
1298  r3A(im, 7)=amet(im,1,2)*gxa1(re,1 )-amet(re,1,2)*gxa1(im,1 )+&
1299 & amet(im,3,2)*gxa1(re,5 )-amet(re,3,2)*gxa1(im,5 )
1300  r3A(im, 8)=amet(im,1,2)*gxa1(re,2 )-amet(re,1,2)*gxa1(im,2 )+&
1301 & amet(im,3,2)*gxa1(re,9 )-amet(re,3,2)*gxa1(im,9 )
1302  r3A(im, 9)=amet(im,1,2)*gxa1(re,3 )-amet(re,1,2)*gxa1(im,3 )+&
1303 & amet(im,3,2)*gxa1(re,10)-amet(re,3,2)*gxa1(im,10)
1304  r3A(im,10)=amet(im,1,2)*gxa1(re,4 )-amet(re,1,2)*gxa1(im,4 )+&
1305 & amet(im,3,2)*gxa1(re,8 )-amet(re,3,2)*gxa1(im,8 )
1306  r3A(im,11)=amet(im,1,2)*gxa1(re,5 )-amet(re,1,2)*gxa1(im,5 )+&
1307 & amet(im,3,2)*gxa1(re,3 )-amet(re,3,2)*gxa1(im,3 )
1308  r3A(im,12)=amet(im,1,2)*gxa1(re,6 )-amet(re,1,2)*gxa1(im,6 )+&
1309 & amet(im,3,2)*gxa1(re,4 )-amet(re,3,2)*gxa1(im,4 )
1310  r3A(im,13)=amet(im,1,3)*gxa1(re,1 )-amet(re,1,3)*gxa1(im,1 )+&
1311 & amet(im,2,3)*gxa1(re,6 )-amet(re,2,3)*gxa1(im,6 )
1312  r3A(im,14)=amet(im,1,3)*gxa1(re,2 )-amet(re,1,3)*gxa1(im,2 )+&
1313 & amet(im,2,3)*gxa1(re,7 )-amet(re,2,3)*gxa1(im,7 )
1314  r3A(im,15)=amet(im,1,3)*gxa1(re,3 )-amet(re,1,3)*gxa1(im,3 )+&
1315 & amet(im,2,3)*gxa1(re,8 )-amet(re,2,3)*gxa1(im,8 )
1316  r3A(im,16)=amet(im,1,3)*gxa1(re,4 )-amet(re,1,3)*gxa1(im,4 )+&
1317 & amet(im,2,3)*gxa1(re,9 )-amet(re,2,3)*gxa1(im,9 )
1318  r3A(im,17)=amet(im,1,3)*gxa1(re,5 )-amet(re,1,3)*gxa1(im,5 )+&
1319 & amet(im,2,3)*gxa1(re,4 )-amet(re,2,3)*gxa1(im,4 )
1320  r3A(im,18)=amet(im,1,3)*gxa1(re,6 )-amet(re,1,3)*gxa1(im,6 )+&
1321 & amet(im,2,3)*gxa1(re,2 )-amet(re,2,3)*gxa1(im,2 )
1322 
1323 !Compute s33(a,b)=2*Re[r3A(a,i,j)*r3G(b,j,i)]
1324 
1325  s33(1)=2.d0*(r3A(re, 1)*r3G(re, 1)-r3A(im, 1)*r3G(im, 1)+&
1326 & r3A(re,12)*r3G(re,12)-r3A(im,12)*r3G(im,12)+&
1327 & r3A(re,17)*r3G(re,17)-r3A(im,17)*r3G(im,17)+&
1328 & r3A(re,11)*r3G(re,18)-r3A(im,11)*r3G(im,18)+&
1329 & r3A(re,18)*r3G(re,11)-r3A(im,18)*r3G(im,11)+&
1330 & r3A(re, 5)*r3G(re,13)-r3A(im, 5)*r3G(im,13)+&
1331 & r3A(re,13)*r3G(re, 5)-r3A(im,13)*r3G(im, 5)+&
1332 & r3A(re, 6)*r3G(re, 7)-r3A(im, 6)*r3G(im, 7)+&
1333 & r3A(re, 7)*r3G(re, 6)-r3A(im, 7)*r3G(im, 6))
1334  s33(2)=2.d0*(r3A(re, 6)*r3G(re, 6)-r3A(im, 6)*r3G(im, 6)+&
1335 & r3A(re, 8)*r3G(re, 8)-r3A(im, 8)*r3G(im, 8)+&
1336 & r3A(re,16)*r3G(re,16)-r3A(im,16)*r3G(im,16)+&
1337 & r3A(re,10)*r3G(re,14)-r3A(im,10)*r3G(im,14)+&
1338 & r3A(re,14)*r3G(re,10)-r3A(im,14)*r3G(im,10)+&
1339 & r3A(re, 4)*r3G(re,18)-r3A(im, 4)*r3G(im,18)+&
1340 & r3A(re,18)*r3G(re, 4)-r3A(im,18)*r3G(im, 4)+&
1341 & r3A(re, 2)*r3G(re,12)-r3A(im, 2)*r3G(im,12)+&
1342 & r3A(re,12)*r3G(re, 2)-r3A(im,12)*r3G(im, 2))
1343  s33(3)=2.d0*(r3A(re, 5)*r3G(re, 5)-r3A(im, 5)*r3G(im, 5)+&
1344 & r3A(re,10)*r3G(re,10)-r3A(im,10)*r3G(im,10)+&
1345 & r3A(re,15)*r3G(re,15)-r3A(im,15)*r3G(im,15)+&
1346 & r3A(re, 9)*r3G(re,16)-r3A(im, 9)*r3G(im,16)+&
1347 & r3A(re,16)*r3G(re, 9)-r3A(im,16)*r3G(im, 9)+&
1348 & r3A(re, 3)*r3G(re,17)-r3A(im, 3)*r3G(im,17)+&
1349 & r3A(re,17)*r3G(re, 3)-r3A(im,17)*r3G(im, 3)+&
1350 & r3A(re, 4)*r3G(re,11)-r3A(im, 4)*r3G(im,11)+&
1351 & r3A(re,11)*r3G(re, 4)-r3A(im,11)*r3G(im, 4))
1352  s33(4)=2.d0*(r3A(re, 5)*r3G(re, 6)-r3A(im, 5)*r3G(im, 6)+&
1353 & r3A(re,10)*r3G(re, 8)-r3A(im,10)*r3G(im, 8)+&
1354 & r3A(re,15)*r3G(re,16)-r3A(im,15)*r3G(im,16)+&
1355 & r3A(re, 9)*r3G(re,14)-r3A(im, 9)*r3G(im,14)+&
1356 & r3A(re,16)*r3G(re,10)-r3A(im,16)*r3G(im,10)+&
1357 & r3A(re, 3)*r3G(re,18)-r3A(im, 3)*r3G(im,18)+&
1358 & r3A(re,17)*r3G(re, 4)-r3A(im,17)*r3G(im, 4)+&
1359 & r3A(re, 4)*r3G(re,12)-r3A(im, 4)*r3G(im,12)+&
1360 & r3A(re,11)*r3G(re, 2)-r3A(im,11)*r3G(im, 2))
1361  s33(5)=2.d0*(r3A(re, 5)*r3G(re, 1)-r3A(im, 5)*r3G(im, 1)+&
1362 & r3A(re,10)*r3G(re,12)-r3A(im,10)*r3G(im,12)+&
1363 & r3A(re,15)*r3G(re,17)-r3A(im,15)*r3G(im,17)+&
1364 & r3A(re, 9)*r3G(re,18)-r3A(im, 9)*r3G(im,18)+&
1365 & r3A(re,16)*r3G(re,11)-r3A(im,16)*r3G(im,11)+&
1366 & r3A(re, 3)*r3G(re,13)-r3A(im, 3)*r3G(im,13)+&
1367 & r3A(re,17)*r3G(re, 5)-r3A(im,17)*r3G(im, 5)+&
1368 & r3A(re, 4)*r3G(re, 7)-r3A(im, 4)*r3G(im, 7)+&
1369 & r3A(re,11)*r3G(re, 6)-r3A(im,11)*r3G(im, 6))
1370  s33(6)=2.d0*(r3A(re, 6)*r3G(re, 1)-r3A(im, 6)*r3G(im, 1)+&
1371 & r3A(re, 8)*r3G(re,12)-r3A(im, 8)*r3G(im,12)+&
1372 & r3A(re,16)*r3G(re,17)-r3A(im,16)*r3G(im,17)+&
1373 & r3A(re,10)*r3G(re,18)-r3A(im,10)*r3G(im,18)+&
1374 & r3A(re,14)*r3G(re,11)-r3A(im,14)*r3G(im,11)+&
1375 & r3A(re, 4)*r3G(re,13)-r3A(im, 4)*r3G(im,13)+&
1376 & r3A(re,18)*r3G(re, 5)-r3A(im,18)*r3G(im, 5)+&
1377 & r3A(re, 2)*r3G(re, 7)-r3A(im, 2)*r3G(im, 7)+&
1378 & r3A(re,12)*r3G(re, 6)-r3A(im,12)*r3G(im, 6))
1379 
1380 !Compute s31(a,b)=2*Re[r3A(a,b,i)*r1(i)]
1381 
1382  s31(1)=2.d0*(r1(re,1)*r3A(re, 1)-r1(im,1)*r3A(im, 1)+&
1383 & r1(re,2)*r3A(re, 7)-r1(im,2)*r3A(im, 7)+&
1384 & r1(re,3)*r3A(re,13)-r1(im,3)*r3A(im,13))
1385  s31(2)=2.d0*(r1(re,1)*r3A(re, 2)-r1(im,1)*r3A(im, 2)+&
1386 & r1(re,2)*r3A(re, 8)-r1(im,2)*r3A(im, 8)+&
1387 & r1(re,3)*r3A(re,14)-r1(im,3)*r3A(im,14))
1388  s31(3)=2.d0*(r1(re,1)*r3A(re, 3)-r1(im,1)*r3A(im, 3)+&
1389 & r1(re,2)*r3A(re, 9)-r1(im,2)*r3A(im, 9)+&
1390 & r1(re,3)*r3A(re,15)-r1(im,3)*r3A(im,15))
1391  s31(4)=2.d0*(r1(re,1)*r3A(re, 4)-r1(im,1)*r3A(im, 4)+&
1392 & r1(re,2)*r3A(re,10)-r1(im,2)*r3A(im,10)+&
1393 & r1(re,3)*r3A(re,16)-r1(im,3)*r3A(im,16))
1394  s31(5)=2.d0*(r1(re,1)*r3A(re, 5)-r1(im,1)*r3A(im, 5)+&
1395 & r1(re,2)*r3A(re,11)-r1(im,2)*r3A(im,11)+&
1396 & r1(re,3)*r3A(re,17)-r1(im,3)*r3A(im,17))
1397  s31(6)=2.d0*(r1(re,1)*r3A(re, 6)-r1(im,1)*r3A(im, 6)+&
1398 & r1(re,2)*r3A(re,12)-r1(im,2)*r3A(im,12)+&
1399 & r1(re,3)*r3A(re,18)-r1(im,3)*r3A(im,18))
1400 
1401 !Finally, compute rank2(a,b)=-15*s33(a,b)+3*s31(a,b)
1402 
1403  rank2(:)=15.d0*s33(:)-3.d0*s31(:)
1404 
1405 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.

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1449 subroutine cont35(gxa,rank5,rank2)
1450 
1451 
1452 !This section has been created automatically by the script Abilint (TD).
1453 !Do not modify the following lines by hand.
1454 #undef ABI_FUNC
1455 #define ABI_FUNC 'cont35'
1456 !End of the abilint section
1457 
1458  implicit none
1459 
1460 !Arguments ------------------------------------
1461 !arrays
1462  real(dp),intent(in) :: gxa(2,10),rank5(2,21)
1463  real(dp),intent(out) :: rank2(6)
1464 
1465 !Local variables-------------------------------
1466 !scalars
1467  integer,parameter :: im=2,re=1
1468 
1469 ! *************************************************************************
1470 
1471 !Simply write out index summations
1472 
1473 !a=1, b=1 in rank2(a,b) --> maps to index 1
1474  rank2(1)=2.0d0*(&
1475 & gxa(re, 1)*rank5(re, 1)+gxa(im, 1)*rank5(im, 1)+&
1476 & gxa(re, 7)*rank5(re, 7)+gxa(im, 7)*rank5(im, 7)+&
1477 & gxa(re,10)*rank5(re,10)+gxa(im,10)*rank5(im,10)+&
1478 & gxa(re, 2)*rank5(re, 2)+gxa(im, 2)*rank5(im, 2)+&
1479 & gxa(re, 3)*rank5(re, 3)+gxa(im, 3)*rank5(im, 3)+&
1480 & gxa(re, 5)*rank5(re, 5)+gxa(im, 5)*rank5(im, 5)+&
1481 & gxa(re, 6)*rank5(re, 6)+gxa(im, 6)*rank5(im, 6)+&
1482 & gxa(re, 8)*rank5(re, 8)+gxa(im, 8)*rank5(im, 8)+&
1483 & gxa(re, 9)*rank5(re, 9)+gxa(im, 9)*rank5(im, 9)+&
1484 & gxa(re, 4)*rank5(re, 4)+gxa(im, 4)*rank5(im, 4))
1485 
1486 
1487 !a=2, b=2 in rank2(a,b) --> maps to index 2
1488  rank2(2)=2.0d0*(&
1489 & gxa(re, 1)*rank5(re, 2)+gxa(im, 1)*rank5(im, 2)+&
1490 & gxa(re, 7)*rank5(re,16)+gxa(im, 7)*rank5(im,16)+&
1491 & gxa(re,10)*rank5(re,19)+gxa(im,10)*rank5(im,19)+&
1492 & gxa(re, 2)*rank5(re,11)+gxa(im, 2)*rank5(im,11)+&
1493 & gxa(re, 3)*rank5(re,12)+gxa(im, 3)*rank5(im,12)+&
1494 & gxa(re, 5)*rank5(re, 9)+gxa(im, 5)*rank5(im, 9)+&
1495 & gxa(re, 6)*rank5(re, 7)+gxa(im, 6)*rank5(im, 7)+&
1496 & gxa(re, 8)*rank5(re,17)+gxa(im, 8)*rank5(im,17)+&
1497 & gxa(re, 9)*rank5(re,18)+gxa(im, 9)*rank5(im,18)+&
1498 & gxa(re, 4)*rank5(re,13)+gxa(im, 4)*rank5(im,13))
1499 
1500 !a=3, b=3 in rank2(a,b) --> maps to index 3
1501  rank2(3)=2.0d0*(&
1502 & gxa(re, 1)*rank5(re, 3)+gxa(im, 1)*rank5(im, 3)+&
1503 & gxa(re, 7)*rank5(re,17)+gxa(im, 7)*rank5(im,17)+&
1504 & gxa(re,10)*rank5(re,21)+gxa(im,10)*rank5(im,21)+&
1505 & gxa(re, 2)*rank5(re,12)+gxa(im, 2)*rank5(im,12)+&
1506 & gxa(re, 3)*rank5(re,15)+gxa(im, 3)*rank5(im,15)+&
1507 & gxa(re, 5)*rank5(re,10)+gxa(im, 5)*rank5(im,10)+&
1508 & gxa(re, 6)*rank5(re, 8)+gxa(im, 6)*rank5(im, 8)+&
1509 & gxa(re, 8)*rank5(re,20)+gxa(im, 8)*rank5(im,20)+&
1510 & gxa(re, 9)*rank5(re,19)+gxa(im, 9)*rank5(im,19)+&
1511 & gxa(re, 4)*rank5(re,14)+gxa(im, 4)*rank5(im,14))
1512 
1513 !a=3, b=2 in rank2(a,b) --> maps to index 4
1514  rank2(4)=2.0d0*(&
1515 & gxa(re, 1)*rank5(re, 4)+gxa(im, 1)*rank5(im, 4)+&
1516 & gxa(re, 7)*rank5(re,18)+gxa(im, 7)*rank5(im,18)+&
1517 & gxa(re,10)*rank5(re,20)+gxa(im,10)*rank5(im,20)+&
1518 & gxa(re, 2)*rank5(re,13)+gxa(im, 2)*rank5(im,13)+&
1519 & gxa(re, 3)*rank5(re,14)+gxa(im, 3)*rank5(im,14)+&
1520 & gxa(re, 5)*rank5(re, 8)+gxa(im, 5)*rank5(im, 8)+&
1521 & gxa(re, 6)*rank5(re, 9)+gxa(im, 6)*rank5(im, 9)+&
1522 & gxa(re, 8)*rank5(re,19)+gxa(im, 8)*rank5(im,19)+&
1523 & gxa(re, 9)*rank5(re,17)+gxa(im, 9)*rank5(im,17)+&
1524 & gxa(re, 4)*rank5(re,12)+gxa(im, 4)*rank5(im,12))
1525 
1526 !a=3, b=1 in rank2(a,b) --> maps to index 5
1527  rank2(5)=2.0d0*(&
1528 & gxa(re, 1)*rank5(re, 5)+gxa(im, 1)*rank5(im, 5)+&
1529 & gxa(re, 7)*rank5(re,13)+gxa(im, 7)*rank5(im,13)+&
1530 & gxa(re,10)*rank5(re,15)+gxa(im,10)*rank5(im,15)+&
1531 & gxa(re, 2)*rank5(re, 9)+gxa(im, 2)*rank5(im, 9)+&
1532 & gxa(re, 3)*rank5(re,10)+gxa(im, 3)*rank5(im,10)+&
1533 & gxa(re, 5)*rank5(re, 3)+gxa(im, 5)*rank5(im, 3)+&
1534 & gxa(re, 6)*rank5(re, 4)+gxa(im, 6)*rank5(im, 4)+&
1535 & gxa(re, 8)*rank5(re,14)+gxa(im, 8)*rank5(im,14)+&
1536 & gxa(re, 9)*rank5(re,12)+gxa(im, 9)*rank5(im,12)+&
1537 & gxa(re, 4)*rank5(re, 8)+gxa(im, 4)*rank5(im, 8))
1538 
1539 !a=2, b=1 in rank2(a,b) --> maps to index 6
1540  rank2(6)=2.0d0*(&
1541 & gxa(re, 1)*rank5(re, 6)+gxa(im, 1)*rank5(im, 6)+&
1542 & gxa(re, 7)*rank5(re,11)+gxa(im, 7)*rank5(im,11)+&
1543 & gxa(re,10)*rank5(re,14)+gxa(im,10)*rank5(im,14)+&
1544 & gxa(re, 2)*rank5(re, 7)+gxa(im, 2)*rank5(im, 7)+&
1545 & gxa(re, 3)*rank5(re, 8)+gxa(im, 3)*rank5(im, 8)+&
1546 & gxa(re, 5)*rank5(re, 4)+gxa(im, 5)*rank5(im, 4)+&
1547 & gxa(re, 6)*rank5(re, 2)+gxa(im, 6)*rank5(im, 2)+&
1548 & gxa(re, 8)*rank5(re,12)+gxa(im, 8)*rank5(im,12)+&
1549 & gxa(re, 9)*rank5(re,13)+gxa(im, 9)*rank5(im,13)+&
1550 & gxa(re, 4)*rank5(re, 9)+gxa(im, 4)*rank5(im, 9))
1551 
1552 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

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1617 subroutine metcon(rank,gmet,aa,bb)
1618 
1619 
1620 !This section has been created automatically by the script Abilint (TD).
1621 !Do not modify the following lines by hand.
1622 #undef ABI_FUNC
1623 #define ABI_FUNC 'metcon'
1624 !End of the abilint section
1625 
1626  implicit none
1627 
1628 !Arguments ------------------------------------
1629 !scalars
1630  integer,intent(in) :: rank
1631 !arrays
1632  real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),gmet(3,3)
1633  real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2)
1634 
1635 !Local variables-------------------------------
1636 !scalars
1637  integer :: ii,jj
1638  real(dp) :: scalar,tmpiii,tmpijk
1639  character(len=500) :: message
1640 !arrays
1641  real(dp) :: vector(3)
1642 
1643 ! *************************************************************************
1644 
1645 !This statement function defines the l=3 contraction in
1646 !terms of the free indices of the contracted tensor (re and im)
1647 ! coniii(ii,i1,i2,i3)=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+&
1648 !& gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+&
1649 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10)
1650 ! conijk(ii,i1,i2,i3)=aa(ii,4)*&
1651 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+&
1652 !& gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+&
1653 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+&
1654 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+&
1655 !& gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+&
1656 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,3))
1657 ! con(ii,i1,i2,i3)=coniii(ii,i1,i2,i3)+conijk(ii,i1,i2,i3)+&
1658 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+&
1659 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+&
1660 !& gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+&
1661 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+&
1662 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+&
1663 !& gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+&
1664 !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+&
1665 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+&
1666 !& gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+&
1667 !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+&
1668 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+&
1669 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+&
1670 !& (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+&
1671 !& gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+&
1672 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+&
1673 !& (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+&
1674 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+&
1675 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8)
1676 
1677 !DEBUG
1678 !write(std_out,*)' metcon : enter '
1679 !stop
1680 !ENDDEBUG
1681  if (rank==0) then
1682 !  Simply copy scalar, re and im
1683    bb(1,1)=aa(1,1)
1684    bb(2,1)=aa(2,1)
1685 
1686  else if (rank==1) then
1687 !  Apply gmet to input vector, re and im
1688    do ii=1,2
1689      do jj=1,3
1690        bb(ii,jj)=gmet(jj,1)*aa(ii,1)+gmet(jj,2)*aa(ii,2)+gmet(jj,3)*aa(ii,3)
1691      end do
1692    end do
1693 
1694 
1695  else if (rank==2) then
1696 !  Apply rank2 expression, re and im
1697    do ii=1,2
1698 !    Carry out g(c,d)*aa(c,d) contraction to get scalar
1699      scalar=gmet(1,1)*aa(ii,1)+gmet(2,2)*aa(ii,2)+&
1700 &     gmet(3,3)*aa(ii,3)+2.0d0*(gmet(2,1)*aa(ii,6)+&
1701 &     gmet(3,1)*aa(ii,5)+gmet(3,2)*aa(ii,4) )
1702 !    Write out components of contraction
1703 !    (1,1)->1
1704      bb(ii,1)=0.5d0*(3.0d0*(gmet(1,1)*gmet(1,1)*aa(ii,1)+&
1705 &     gmet(1,2)*gmet(1,2)*aa(ii,2)+gmet(1,3)*gmet(1,3)*aa(ii,3)+&
1706 &     (gmet(1,2)*gmet(1,1)+gmet(1,1)*gmet(1,2))*aa(ii,6)+&
1707 &     (gmet(1,3)*gmet(1,1)+gmet(1,1)*gmet(1,3))*aa(ii,5)+&
1708 &     (gmet(1,3)*gmet(1,2)+gmet(1,2)*gmet(1,3))*aa(ii,4) ) &
1709 &     - gmet(1,1)*scalar)
1710 !    (2,2)->2
1711      bb(ii,2)=0.5d0*(3.0d0*(gmet(2,1)*gmet(2,1)*aa(ii,1)+&
1712 &     gmet(2,2)*gmet(2,2)*aa(ii,2)+gmet(2,3)*gmet(2,3)*aa(ii,3)+&
1713 &     (gmet(2,2)*gmet(2,1)+gmet(2,1)*gmet(2,2))*aa(ii,6)+&
1714 &     (gmet(2,3)*gmet(2,1)+gmet(2,1)*gmet(2,3))*aa(ii,5)+&
1715 &     (gmet(2,3)*gmet(2,2)+gmet(2,2)*gmet(2,3))*aa(ii,4) )&
1716 &     - gmet(2,2)*scalar)
1717 !    (3,3)->3
1718      bb(ii,3)=0.5d0*(3.0d0*(gmet(3,1)*gmet(3,1)*aa(ii,1)+&
1719 &     gmet(3,2)*gmet(3,2)*aa(ii,2)+gmet(3,3)*gmet(3,3)*aa(ii,3)+&
1720 &     (gmet(3,2)*gmet(3,1)+gmet(3,1)*gmet(3,2))*aa(ii,6)+&
1721 &     (gmet(3,3)*gmet(3,1)+gmet(3,1)*gmet(3,3))*aa(ii,5)+&
1722 &     (gmet(3,3)*gmet(3,2)+gmet(3,2)*gmet(3,3))*aa(ii,4) )&
1723 &     - gmet(3,3)*scalar)
1724 !    (3,2)->4
1725      bb(ii,4)=0.5d0*(3.0d0*(gmet(3,1)*gmet(2,1)*aa(ii,1)+&
1726 &     gmet(3,2)*gmet(2,2)*aa(ii,2)+gmet(3,3)*gmet(2,3)*aa(ii,3)+&
1727 &     (gmet(3,2)*gmet(2,1)+gmet(3,1)*gmet(2,2))*aa(ii,6)+&
1728 &     (gmet(3,3)*gmet(2,1)+gmet(3,1)*gmet(2,3))*aa(ii,5)+&
1729 &     (gmet(3,3)*gmet(2,2)+gmet(3,2)*gmet(2,3))*aa(ii,4) )&
1730 &     - gmet(3,2)*scalar)
1731 !    (3,1)->5
1732      bb(ii,5)=0.5d0*(3.0d0*(gmet(3,1)*gmet(1,1)*aa(ii,1)+&
1733 &     gmet(3,2)*gmet(1,2)*aa(ii,2)+gmet(3,3)*gmet(1,3)*aa(ii,3)+&
1734 &     (gmet(3,2)*gmet(1,1)+gmet(3,1)*gmet(1,2))*aa(ii,6)+&
1735 &     (gmet(3,3)*gmet(1,1)+gmet(3,1)*gmet(1,3))*aa(ii,5)+&
1736 &     (gmet(3,3)*gmet(1,2)+gmet(3,2)*gmet(1,3))*aa(ii,4) )&
1737 &     - gmet(3,1)*scalar)
1738 !    (2,1)->6
1739      bb(ii,6)=0.5d0*(3.0d0*(gmet(2,1)*gmet(1,1)*aa(ii,1)+&
1740 &     gmet(2,2)*gmet(1,2)*aa(ii,2)+gmet(2,3)*gmet(1,3)*aa(ii,3)+&
1741 &     (gmet(2,2)*gmet(1,1)+gmet(2,1)*gmet(1,2))*aa(ii,6)+&
1742 &     (gmet(2,3)*gmet(1,1)+gmet(2,1)*gmet(1,3))*aa(ii,5)+&
1743 &     (gmet(2,3)*gmet(1,2)+gmet(2,2)*gmet(1,3))*aa(ii,4) ) &
1744 &     - gmet(2,1)*scalar)
1745 !    Include appropriate weights for multiplicity
1746      bb(ii,4)=2.d0*bb(ii,4)
1747      bb(ii,5)=2.d0*bb(ii,5)
1748      bb(ii,6)=2.d0*bb(ii,6)
1749    end do
1750 
1751  else if (rank==3) then
1752 !  Apply rank2 expression, re and im
1753    do ii=1,2
1754 !    Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j)
1755      do jj=1,3
1756        tmpiii=   gmet(1,1)*gmet(jj,1)*aa(ii,1)+&
1757 &       gmet(2,2)*gmet(jj,2)*aa(ii,7)+&
1758 &       gmet(3,3)*gmet(jj,3)*aa(ii,10)
1759        tmpijk=  (gmet(1,2)*gmet(jj,3)+&
1760 &       gmet(3,1)*gmet(jj,2)+&
1761 &       gmet(2,3)*gmet(jj,1)+&
1762 &       gmet(3,2)*gmet(jj,1)+&
1763 &       gmet(1,3)*gmet(jj,2)+&
1764 &       gmet(2,1)*gmet(jj,3)) *aa(ii,4)
1765        vector(jj)=tmpiii + tmpijk +&
1766 &       (gmet(1,2)*gmet(jj,1)+&
1767 &       gmet(2,1)*gmet(jj,1)+&
1768 &       gmet(1,1)*gmet(jj,2)) *aa(ii,6)+&
1769 &       (gmet(1,2)*gmet(jj,2)+&
1770 &       gmet(2,1)*gmet(jj,2)+&
1771 &       gmet(2,2)*gmet(jj,1)) *aa(ii,2)+&
1772 &       (gmet(1,3)*gmet(jj,1)+&
1773 &       gmet(3,1)*gmet(jj,1)+&
1774 &       gmet(1,1)*gmet(jj,3)) *aa(ii,5)+&
1775 &       (gmet(1,3)*gmet(jj,3)+&
1776 &       gmet(3,1)*gmet(jj,3)+&
1777 &       gmet(3,3)*gmet(jj,1)) *aa(ii,3)+&
1778 &       (gmet(2,3)*gmet(jj,2)+&
1779 &       gmet(3,2)*gmet(jj,2)+&
1780 &       gmet(2,2)*gmet(jj,3)) *aa(ii,9)+&
1781 &       (gmet(2,3)*gmet(jj,3)+&
1782 &       gmet(3,2)*gmet(jj,3)+&
1783 &       gmet(3,3)*gmet(jj,2)) *aa(ii,8)
1784      end do
1785 !    Write out components of contraction
1786 !    (111)->1
1787      bb(ii,1) =2.5d0*con_met(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1))
1788 !    (221)->2
1789      bb(ii,2) =2.5d0*con_met(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+&
1790 &     gmet(1,2)*vector(2)+gmet(2,2)*vector(1))
1791 !    (331)->3
1792      bb(ii,3) =2.5d0*con_met(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+&
1793 &     gmet(1,3)*vector(3)+gmet(3,3)*vector(1))
1794 !    (321)->4
1795      bb(ii,4) =2.5d0*con_met(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+&
1796 &     gmet(1,2)*vector(3)+gmet(3,2)*vector(1))
1797 !    (311)->5
1798      bb(ii,5) =2.5d0*con_met(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+&
1799 &     gmet(1,1)*vector(3)+gmet(3,1)*vector(1))
1800 !    (211)->6
1801      bb(ii,6) =2.5d0*con_met(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+&
1802 &     gmet(1,1)*vector(2)+gmet(2,1)*vector(1))
1803 !    (222)->7
1804      bb(ii,7) =2.5d0*con_met(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2))
1805 
1806 !    (332)->8
1807      bb(ii,8) =2.5d0*con_met(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+&
1808 &     gmet(2,3)*vector(3)+gmet(3,3)*vector(2))
1809 !    (322)->9
1810      bb(ii,9) =2.5d0*con_met(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+&
1811 &     gmet(2,2)*vector(3)+gmet(3,2)*vector(2))
1812 !    (333)->10
1813      bb(ii,10)=2.5d0*con_met(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3))
1814 !    Include appropriate weights for multiplicity
1815      bb(ii,2)=3.d0*bb(ii,2)
1816      bb(ii,3)=3.d0*bb(ii,3)
1817      bb(ii,4)=6.d0*bb(ii,4)
1818      bb(ii,5)=3.d0*bb(ii,5)
1819      bb(ii,6)=3.d0*bb(ii,6)
1820      bb(ii,8)=3.d0*bb(ii,8)
1821      bb(ii,9)=3.d0*bb(ii,9)
1822    end do
1823 
1824  else
1825    write(message, '(a,i0,a,a,a)' )&
1826 &   'Input rank=',rank,' not allowed.',ch10,&
1827 &   'Possible values are 0,1,2,3 only.'
1828    MSG_BUG(message)
1829  end if
1830 
1831  contains
1832 
1833    function con_met(ii,i1,i2,i3)
1834 
1835 
1836 !This section has been created automatically by the script Abilint (TD).
1837 !Do not modify the following lines by hand.
1838 #undef ABI_FUNC
1839 #define ABI_FUNC 'con_met'
1840 !End of the abilint section
1841 
1842    real(dp) :: con_met
1843    integer :: ii,i1,i2,i3
1844    real(dp)::coniii,conijk
1845 
1846    coniii=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+&
1847 &   gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+&
1848 &   gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10)
1849    conijk=aa(ii,4)*&
1850 &   (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+&
1851 &   gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+&
1852 &   gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+&
1853 &   gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+&
1854 &   gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+&
1855 &   gmet(i1,2)*gmet(i2,1)*gmet(i3,3))
1856    con_met=coniii+conijk+&
1857 &   (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+&
1858 &   gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+&
1859 &   gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+&
1860 &   (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+&
1861 &   gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+&
1862 &   gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+&
1863 &   (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+&
1864 &   gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+&
1865 &   gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+&
1866 &   (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+&
1867 &   gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+&
1868 &   gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+&
1869 &   (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+&
1870 &   gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+&
1871 &   gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+&
1872 &   (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+&
1873 &   gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+&
1874 &   gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8)
1875 
1876  end function con_met
1877 
1878 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

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1946 subroutine metcon_so(rank,gmet,amet,aa,bb)
1947 
1948 
1949 !This section has been created automatically by the script Abilint (TD).
1950 !Do not modify the following lines by hand.
1951 #undef ABI_FUNC
1952 #define ABI_FUNC 'metcon_so'
1953 !End of the abilint section
1954 
1955  implicit none
1956 
1957 !Arguments ------------------------------------
1958 !scalars
1959  integer,intent(in) :: rank
1960 !arrays
1961  real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),amet(3,3),gmet(3,3)
1962  real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2)
1963 
1964 !Local variables-------------------------------
1965 !scalars
1966  integer :: ii,jj
1967  real(dp) :: tmpiii,tmpijk
1968  character(len=500) :: message
1969 !arrays
1970  real(dp) :: vector(3)
1971 
1972 ! *************************************************************************
1973 
1974  if (rank==0) then
1975 !  Simply copy scalar, re and im
1976    bb(1,1)=0.d0
1977    bb(2,1)=0.d0
1978 !
1979  else if (rank==1) then
1980 !  Apply gmet to input vector, re and im
1981    do ii=1,2
1982      do jj=1,3
1983        bb(ii,jj)=amet(jj,1)*aa(ii,1)+amet(jj,2)*aa(ii,2)+amet(jj,3)*aa(ii,3)
1984      end do
1985    end do
1986 
1987  else if (rank==2) then
1988 !  Apply rank2 expression, re and im
1989    do ii=1,2
1990 !    Write out components of contraction
1991 !    (1,1)->1
1992      bb(ii,1)=3.0d0*(gmet(1,1)*amet(1,1)*aa(ii,1)+&
1993 &     gmet(1,2)*amet(1,2)*aa(ii,2)+gmet(1,3)*amet(1,3)*aa(ii,3)+&
1994 &     (gmet(1,2)*amet(1,1)+gmet(1,1)*amet(1,2))*aa(ii,6)+&
1995 &     (gmet(1,3)*amet(1,1)+gmet(1,1)*amet(1,3))*aa(ii,5)+&
1996 &     (gmet(1,3)*amet(1,2)+gmet(1,2)*amet(1,3))*aa(ii,4))
1997 !    (2,2)->2
1998      bb(ii,2)=3.0d0*(gmet(2,1)*amet(2,1)*aa(ii,1)+&
1999 &     gmet(2,2)*amet(2,2)*aa(ii,2)+gmet(2,3)*amet(2,3)*aa(ii,3)+&
2000 &     (gmet(2,2)*amet(2,1)+gmet(2,1)*amet(2,2))*aa(ii,6)+&
2001 &     (gmet(2,3)*amet(2,1)+gmet(2,1)*amet(2,3))*aa(ii,5)+&
2002 &     (gmet(2,3)*amet(2,2)+gmet(2,2)*amet(2,3))*aa(ii,4) )
2003 !    (3,3)->3
2004      bb(ii,3)=3.0d0*(gmet(3,1)*amet(3,1)*aa(ii,1)+&
2005 &     gmet(3,2)*amet(3,2)*aa(ii,2)+gmet(3,3)*amet(3,3)*aa(ii,3)+&
2006 &     (gmet(3,2)*amet(3,1)+gmet(3,1)*amet(3,2))*aa(ii,6)+&
2007 &     (gmet(3,3)*amet(3,1)+gmet(3,1)*amet(3,3))*aa(ii,5)+&
2008 &     (gmet(3,3)*amet(3,2)+gmet(3,2)*amet(3,3))*aa(ii,4) )
2009 !    (3,2)->4
2010      bb(ii,4)=3.0d0*(gmet(3,1)*amet(2,1)*aa(ii,1)+&
2011 &     gmet(3,2)*amet(2,2)*aa(ii,2)+gmet(3,3)*amet(2,3)*aa(ii,3)+&
2012 &     (gmet(3,2)*amet(2,1)+gmet(3,1)*amet(2,2))*aa(ii,6)+&
2013 &     (gmet(3,3)*amet(2,1)+gmet(3,1)*amet(2,3))*aa(ii,5)+&
2014 &     (gmet(3,3)*amet(2,2)+gmet(3,2)*amet(2,3))*aa(ii,4) )
2015      bb(ii,4)=bb(ii,4)+3.0d0*(amet(3,1)*gmet(2,1)*aa(ii,1)+&
2016 &     amet(3,2)*gmet(2,2)*aa(ii,2)+amet(3,3)*gmet(2,3)*aa(ii,3)+&
2017 &     (amet(3,2)*gmet(2,1)+amet(3,1)*gmet(2,2))*aa(ii,6)+&
2018 &     (amet(3,3)*gmet(2,1)+amet(3,1)*gmet(2,3))*aa(ii,5)+&
2019 &     (amet(3,3)*gmet(2,2)+amet(3,2)*gmet(2,3))*aa(ii,4) )
2020 !    (3,1)->5
2021      bb(ii,5)=3.0d0*(gmet(3,1)*amet(1,1)*aa(ii,1)+&
2022 &     gmet(3,2)*amet(1,2)*aa(ii,2)+gmet(3,3)*amet(1,3)*aa(ii,3)+&
2023 &     (gmet(3,2)*amet(1,1)+gmet(3,1)*amet(1,2))*aa(ii,6)+&
2024 &     (gmet(3,3)*amet(1,1)+gmet(3,1)*amet(1,3))*aa(ii,5)+&
2025 &     (gmet(3,3)*amet(1,2)+gmet(3,2)*amet(1,3))*aa(ii,4) )
2026      bb(ii,5)=bb(ii,5)+3.0d0*(amet(3,1)*gmet(1,1)*aa(ii,1)+&
2027 &     amet(3,2)*gmet(1,2)*aa(ii,2)+amet(3,3)*gmet(1,3)*aa(ii,3)+&
2028 &     (amet(3,2)*gmet(1,1)+amet(3,1)*gmet(1,2))*aa(ii,6)+&
2029 &     (amet(3,3)*gmet(1,1)+amet(3,1)*gmet(1,3))*aa(ii,5)+&
2030 &     (amet(3,3)*gmet(1,2)+amet(3,2)*gmet(1,3))*aa(ii,4) )
2031 !    (2,1)->6
2032      bb(ii,6)=3.0d0*(gmet(2,1)*amet(1,1)*aa(ii,1)+&
2033 &     gmet(2,2)*amet(1,2)*aa(ii,2)+gmet(2,3)*amet(1,3)*aa(ii,3)+&
2034 &     (gmet(2,2)*amet(1,1)+gmet(2,1)*amet(1,2))*aa(ii,6)+&
2035 &     (gmet(2,3)*amet(1,1)+gmet(2,1)*amet(1,3))*aa(ii,5)+&
2036 &     (gmet(2,3)*amet(1,2)+gmet(2,2)*amet(1,3))*aa(ii,4) )
2037      bb(ii,6)=bb(ii,6)+3.0d0*(amet(2,1)*gmet(1,1)*aa(ii,1)+&
2038 &     amet(2,2)*gmet(1,2)*aa(ii,2)+amet(2,3)*gmet(1,3)*aa(ii,3)+&
2039 &     (amet(2,2)*gmet(1,1)+amet(2,1)*gmet(1,2))*aa(ii,6)+&
2040 &     (amet(2,3)*gmet(1,1)+amet(2,1)*gmet(1,3))*aa(ii,5)+&
2041 &     (amet(2,3)*gmet(1,2)+amet(2,2)*gmet(1,3))*aa(ii,4) )
2042 !    Include appropriate weights for multiplicity
2043      bb(ii,4)=bb(ii,4)
2044      bb(ii,5)=bb(ii,5)
2045      bb(ii,6)=bb(ii,6)
2046    end do
2047 
2048  else if (rank==3) then
2049 !  Apply rank2 expression, re and im
2050    do ii=1,2
2051 !    Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j)
2052      do jj=1,3
2053        tmpiii=   gmet(1,1)*amet(jj,1)*aa(ii,1)+&
2054 &       gmet(2,2)*amet(jj,2)*aa(ii,7)+&
2055 &       gmet(3,3)*amet(jj,3)*aa(ii,10)
2056        tmpijk=  (gmet(1,2)*amet(jj,3)+&
2057 &       gmet(3,1)*amet(jj,2)+&
2058 &       gmet(2,3)*amet(jj,1)+&
2059 &       gmet(3,2)*amet(jj,1)+&
2060 &       gmet(1,3)*amet(jj,2)+&
2061 &       gmet(2,1)*amet(jj,3)) *aa(ii,4)
2062        vector(jj)=tmpiii + tmpijk +&
2063 &       (gmet(1,2)*amet(jj,1)+&
2064 &       gmet(2,1)*amet(jj,1)+&
2065 &       gmet(1,1)*amet(jj,2)) *aa(ii,6)+&
2066 &       (gmet(1,2)*amet(jj,2)+&
2067 &       gmet(2,1)*amet(jj,2)+&
2068 &       gmet(2,2)*amet(jj,1)) *aa(ii,2)+&
2069 &       (gmet(1,3)*amet(jj,1)+&
2070 &       gmet(3,1)*amet(jj,1)+&
2071 &       gmet(1,1)*amet(jj,3)) *aa(ii,5)+&
2072 &       (gmet(1,3)*amet(jj,3)+&
2073 &       gmet(3,1)*amet(jj,3)+&
2074 &       gmet(3,3)*amet(jj,1)) *aa(ii,3)+&
2075 &       (gmet(2,3)*amet(jj,2)+&
2076 &       gmet(3,2)*amet(jj,2)+&
2077 &       gmet(2,2)*amet(jj,3)) *aa(ii,9)+&
2078 &       (gmet(2,3)*amet(jj,3)+&
2079 &       gmet(3,2)*amet(jj,3)+&
2080 &       gmet(3,3)*amet(jj,2)) *aa(ii,8)
2081      end do
2082 !    Write out components of contraction
2083 !    (111)->1
2084      bb(ii,1) =7.5d0*con_metso(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1))
2085 !    (221)->2
2086      bb(ii,2) =7.5d0*con_metso(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+&
2087 &     gmet(1,2)*vector(2)+gmet(2,2)*vector(1))
2088 !    (331)->3
2089      bb(ii,3) =7.5d0*con_metso(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+&
2090 &     gmet(1,3)*vector(3)+gmet(3,3)*vector(1))
2091 !    (321)->4
2092      bb(ii,4) =7.5d0*con_metso(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+&
2093 &     gmet(1,2)*vector(3)+gmet(3,2)*vector(1))
2094 !    (311)->5
2095      bb(ii,5) =7.5d0*con_metso(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+&
2096 &     gmet(1,1)*vector(3)+gmet(3,1)*vector(1))
2097 !    (211)->6
2098      bb(ii,6) =7.5d0*con_metso(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+&
2099 &     gmet(1,1)*vector(2)+gmet(2,1)*vector(1))
2100 !    (222)->7
2101      bb(ii,7) =7.5d0*con_metso(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2))
2102 
2103 !    (332)->8
2104      bb(ii,8) =7.5d0*con_metso(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+&
2105 &     gmet(2,3)*vector(3)+gmet(3,3)*vector(2))
2106 !    (322)->9
2107      bb(ii,9) =7.5d0*con_metso(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+&
2108 &     gmet(2,2)*vector(3)+gmet(3,2)*vector(2))
2109 !    (333)->10
2110      bb(ii,10)=7.5d0*con_metso(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3))
2111 !    Include appropriate weights for multiplicity
2112      bb(ii,2)=3.d0*bb(ii,2)
2113      bb(ii,3)=3.d0*bb(ii,3)
2114      bb(ii,4)=6.d0*bb(ii,4)
2115      bb(ii,5)=3.d0*bb(ii,5)
2116      bb(ii,6)=3.d0*bb(ii,6)
2117      bb(ii,8)=3.d0*bb(ii,8)
2118      bb(ii,9)=3.d0*bb(ii,9)
2119    end do
2120 
2121  else
2122    write(message, '(a,i0,a,a,a)' )&
2123 &   'Input rank=',rank,' not allowed.',ch10,&
2124 &   'Possible values are 0,1,2,3 only.'
2125    MSG_BUG(message)
2126  end if
2127 
2128  contains
2129 
2130 !This function defines the l=3 contraction in
2131 !terms of the free indices of the contracted tensor (re and im)
2132 
2133    function cona_metso(ii,i1,i2,i3)
2134 
2135 
2136 !This section has been created automatically by the script Abilint (TD).
2137 !Do not modify the following lines by hand.
2138 #undef ABI_FUNC
2139 #define ABI_FUNC 'cona_metso'
2140 !End of the abilint section
2141 
2142    real(dp) :: cona_metso
2143    integer,intent(in) :: ii,i1,i2,i3
2144    real(dp) :: coniii, conijk
2145 
2146    coniii=gmet(i1,1)*gmet(i2,1)*amet(i3,1)*aa(ii,1)+&
2147 &   gmet(i1,2)*gmet(i2,2)*amet(i3,2)*aa(ii,7)+&
2148 &   gmet(i1,3)*gmet(i2,3)*amet(i3,3)*aa(ii,10)
2149    conijk=aa(ii,4)*&
2150 &   (gmet(i1,1)*gmet(i2,2)*amet(i3,3)+&
2151 &   gmet(i1,2)*gmet(i2,3)*amet(i3,1)+&
2152 &   gmet(i1,3)*gmet(i2,1)*amet(i3,2)+&
2153 &   gmet(i1,3)*gmet(i2,2)*amet(i3,1)+&
2154 &   gmet(i1,1)*gmet(i2,3)*amet(i3,2)+&
2155 &   gmet(i1,2)*gmet(i2,1)*amet(i3,3))
2156    cona_metso=coniii+conijk+&
2157 &   (gmet(i1,1)*gmet(i2,2)*amet(i3,1)+&
2158 &   gmet(i1,2)*gmet(i2,1)*amet(i3,1)+&
2159 &   gmet(i1,1)*gmet(i2,1)*amet(i3,2))*aa(ii,6)+&
2160 &   (gmet(i1,1)*gmet(i2,2)*amet(i3,2)+&
2161 &   gmet(i1,2)*gmet(i2,1)*amet(i3,2)+&
2162 &   gmet(i1,2)*gmet(i2,2)*amet(i3,1))*aa(ii,2)+&
2163 &   (gmet(i1,1)*gmet(i2,3)*amet(i3,1)+&
2164 &   gmet(i1,3)*gmet(i2,1)*amet(i3,1)+&
2165 &   gmet(i1,1)*gmet(i2,1)*amet(i3,3))*aa(ii,5)+&
2166 &   (gmet(i1,1)*gmet(i2,3)*amet(i3,3)+&
2167 &   gmet(i1,3)*gmet(i2,1)*amet(i3,3)+&
2168 &   gmet(i1,3)*gmet(i2,3)*amet(i3,1))*aa(ii,3)+&
2169 &   (gmet(i1,2)*gmet(i2,2)*amet(i3,3)+&
2170 &   gmet(i1,2)*gmet(i2,3)*amet(i3,2)+&
2171 &   gmet(i1,3)*gmet(i2,2)*amet(i3,2))*aa(ii,9)+&
2172 &   (gmet(i1,2)*gmet(i2,3)*amet(i3,3)+&
2173 &   gmet(i1,3)*gmet(i2,2)*amet(i3,3)+&
2174 &   gmet(i1,3)*gmet(i2,3)*amet(i3,2))*aa(ii,8)
2175  end function cona_metso
2176 
2177 
2178    function con_metso(ii,i1,i2,i3)
2179 
2180 
2181 !This section has been created automatically by the script Abilint (TD).
2182 !Do not modify the following lines by hand.
2183 #undef ABI_FUNC
2184 #define ABI_FUNC 'con_metso'
2185 !End of the abilint section
2186 
2187    real(dp) :: con_metso
2188    integer,intent(in) :: ii,i1,i2,i3
2189 
2190    con_metso=(cona_metso(ii,i3,i1,i2)+cona_metso(ii,i2,i3,i1)+cona_metso(ii,i1,i2,i3))/3.d0
2191 
2192  end function con_metso
2193 
2194 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

PARENTS

      nonlop_pl

CHILDREN

SOURCE

2241 subroutine metric_so(amet,gprimd,pauli)
2242 
2243 
2244 !This section has been created automatically by the script Abilint (TD).
2245 !Do not modify the following lines by hand.
2246 #undef ABI_FUNC
2247 #define ABI_FUNC 'metric_so'
2248 !End of the abilint section
2249 
2250  implicit none
2251 
2252 !Arguments ------------------------------------
2253 !arrays
2254  real(dp),intent(in) :: gprimd(3,3)
2255  real(dp),intent(out) :: amet(2,3,3,2,2),pauli(2,2,2,3)
2256 
2257 !Local variables-------------------------------
2258 !scalars
2259  integer :: iy1,iy2,m1,m2,n
2260 !arrays
2261  real(dp) :: buffer1(3,3,2,2) !,buffer2(3,3,3,3,2,2)
2262 
2263 ! **********************************************************************
2264 
2265 !Fill in Pauli matrices and make them spin matrices:
2266 !Pauli(Re/Im,up,down,coord):
2267  pauli(:,:,:,:)=0.d0
2268  pauli(1,1,2,1)= 1.d0;pauli(1,2,1,1)= 1.d0
2269  pauli(2,1,2,2)=-1.d0;pauli(2,2,1,2)= 1.d0
2270  pauli(1,1,1,3)= 1.d0;pauli(1,2,2,3)=-1.d0
2271  pauli(:,:,:,:)= 0.5d0*pauli(:,:,:,:)
2272 
2273 !Construct the antisymmetric tensor:
2274  amet(:,:,:,:,:)=0.d0
2275  do iy2=1,3
2276    do iy1=1,3
2277      do n=1,3
2278        m1=mod(n ,3)+1    !  n,m1,m2 is an even permutation
2279        m2=mod(m1,3)+1
2280        amet(1:2,iy1,iy2,1:2,1:2) = amet(:,iy1,iy2,:,:) &
2281 &       + pauli(:,:,:,n) &
2282 &       *(gprimd(m1,iy1)*gprimd(m2,iy2) &
2283 &       -gprimd(m2,iy1)*gprimd(m1,iy2))
2284      end do
2285    end do
2286  end do
2287 !amet= -i amet
2288  buffer1(:,:,:,:)=-amet(1,:,:,:,:)
2289  amet(1,:,:,:,:) = amet(2,:,:,:,:)
2290  amet(2,:,:,:,:) =buffer1(:,:,:,:)
2291 
2292 !DEBUG
2293 !!Eventually construct the gradients of Amet wrt strains:
2294 !! DAmet(y,y'',a,b,s,s'')= d[Amet(y,y'',s,s'')]/d[strain(a,b)]
2295 !!                        -i Pauli(s,s'',n)*
2296 !!                          ( E(n,a,m'')*Gprimd(b,y )*Gprimd(m'',y'')
2297 !!                           +E(n,m,a )*Gprimd(b,y'')*Gprimd(m ,y ) )
2298 !damet(:,:,:,:,:,:,:)=0.d0
2299 !do ib=1,3
2300 !do ia=1,3
2301 !m1=mod(ia,3)+1    !  ia,m1,m2 is an even permutation
2302 !m2=mod(m1,3)+1
2303 !do iy2=1,3
2304 !do iy1=1,3
2305 !damet(:,iy1,iy2,ia,ib,:,:) = damet(:,iy1,iy2,ia,ib,:,:) &
2306 !&        + (pauli(:,:,:,m2)*gprimd(m1,iy2) &
2307 !-pauli(:,:,:,m1)*gprimd(m2,iy2))*gprimd(ib,iy1) &
2308 !&        + (pauli(:,:,:,m1)*gprimd(m2,iy1) &
2309 !-pauli(:,:,:,m2)*gprimd(m1,iy1))*gprimd(ib,iy2)
2310 !end do
2311 !end do
2312 !end do
2313 !end do
2314 !! damet= i damet
2315 !buffer2(:,:,:,:,:,:)= damet(1,:,:,:,:,:,:)
2316 !damet(1,:,:,:,:,:,:)= -damet(2,:,:,:,:,:,:)
2317 !damet(2,:,:,:,:,:,:)=buffer2(:,:,:,:,:,:)
2318 !! Symetrize damet(:,:,:,a,b,:,:)
2319 !damet(:,:,:,1,2,:,:)=0.5d0*(damet(:,:,:,1,2,:,:)+damet(:,:,:,2,1,:,:))
2320 !damet(:,:,:,1,3,:,:)=0.5d0*(damet(:,:,:,1,3,:,:)+damet(:,:,:,3,1,:,:))
2321 !damet(:,:,:,2,3,:,:)=0.5d0*(damet(:,:,:,2,3,:,:)+damet(:,:,:,3,2,:,:))
2322 !damet(:,:,:,2,1,:,:)=damet(:,:,:,1,2,:,:)
2323 !damet(:,:,:,3,1,:,:)=damet(:,:,:,1,3,:,:)
2324 !damet(:,:,:,3,2,:,:)=damet(:,:,:,2,3,:,:)
2325 !ENDDEBUG
2326 
2327 end subroutine metric_so