TABLE OF CONTENTS
- ABINIT/m_contract
- m_contract/cont13
- m_contract/cont22
- m_contract/cont22cso
- m_contract/cont22so
- m_contract/cont24
- m_contract/cont3
- m_contract/cont33cso
- m_contract/cont33so
- m_contract/cont35
- m_contract/metcon
- m_contract/metcon_so
- m_contract/metric_so
ABINIT/m_contract [ Modules ]
NAME
m_contract
FUNCTION
Low-level procedeures used in nonlop_pl to contract tensors
COPYRIGHT
Copyright (C) 1998-2022 ABINIT group (DCA, XG, MT, GZ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_contract 23 24 use defs_basis 25 use m_errors 26 27 implicit none 28 29 private
m_contract/cont13 [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont13
FUNCTION
Contract rank1 tensor with rank3 symmetric tensor to produce symmetric rank2 tensor
INPUTS
rank1(2,3)=rank 1 complex tensor (vector of length 3) rank3(2,10)=rank 3 complex tensor (symmetric storage)
OUTPUT
rank2(6)=rank 2 real tensor (symmetric storage)
NOTES
Tensors are in "symmetric" storage mode. For rank1 this is 1, 2, 3; for rank2 this is 11, 22, 33, 32, 31, 21; for rank3 this is 111, 221, 331, 321, 311, 211, 222, 332, 322, 333. rank1 and rank3 are complex; rank2 is real. Want $2 Re[contraction]$. $rank2(a,b)=2 Re[rank1(i)^"*" rank3(a,b,i)]$. In typical usage the input rank1 tensor is actually $rank1(i)=gmet(i,j) gxa(j)$
SOURCE
77 subroutine cont13(rank1,rank3,rank2) 78 79 !Arguments ------------------------------------ 80 !arrays 81 real(dp),intent(in) :: rank1(2,3),rank3(2,10) 82 real(dp),intent(out) :: rank2(6) 83 84 !Local variables------------------------------- 85 !scalars 86 integer,parameter :: im=2,re=1 87 88 ! ************************************************************************* 89 90 !Simply write out index summations 91 !a=1, b=1 in rank2(a,b) --> maps to index 1 92 rank2(1)=2.0d0*(& 93 & (rank1(re,1)*rank3(re,1)+rank1(im,1)*rank3(im,1))+& 94 & (rank1(re,2)*rank3(re,6)+rank1(im,2)*rank3(im,6))+& 95 & (rank1(re,3)*rank3(re,5)+rank1(im,3)*rank3(im,5))) 96 97 !a=2, b=2 in rank2(a,b) --> maps to index 2 98 rank2(2)=2.0d0*(& 99 & (rank1(re,1)*rank3(re,2)+rank1(im,1)*rank3(im,2))+& 100 & (rank1(re,2)*rank3(re,7)+rank1(im,2)*rank3(im,7))+& 101 & (rank1(re,3)*rank3(re,9)+rank1(im,3)*rank3(im,9))) 102 103 !a=3, b=3 in rank2(a,b) --> maps to index 3 104 rank2(3)=2.0d0*(& 105 & (rank1(re,1)*rank3(re,3)+rank1(im,1)*rank3(im,3))+& 106 & (rank1(re,2)*rank3(re,8)+rank1(im,2)*rank3(im,8))+& 107 & (rank1(re,3)*rank3(re,10)+rank1(im,3)*rank3(im,10))) 108 109 !a=3, b=2 in rank2(a,b) --> maps to index 4 110 rank2(4)=2.0d0*(& 111 & (rank1(re,1)*rank3(re,4)+rank1(im,1)*rank3(im,4))+& 112 & (rank1(re,2)*rank3(re,9)+rank1(im,2)*rank3(im,9))+& 113 & (rank1(re,3)*rank3(re,8)+rank1(im,3)*rank3(im,8))) 114 115 !a=3, b=1 in rank2(a,b) --> maps to index 5 116 rank2(5)=2.0d0*(& 117 & (rank1(re,1)*rank3(re,5)+rank1(im,1)*rank3(im,5))+& 118 & (rank1(re,2)*rank3(re,4)+rank1(im,2)*rank3(im,4))+& 119 & (rank1(re,3)*rank3(re,3)+rank1(im,3)*rank3(im,3))) 120 121 !a=2, b=1 in rank2(a,b) --> maps to index 6 122 rank2(6)=2.0d0*(& 123 & (rank1(re,1)*rank3(re,6)+rank1(im,1)*rank3(im,6))+& 124 & (rank1(re,2)*rank3(re,2)+rank1(im,2)*rank3(im,2))+& 125 & (rank1(re,3)*rank3(re,4)+rank1(im,3)*rank3(im,4))) 126 127 end subroutine cont13
m_contract/cont22 [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont22
FUNCTION
Contract symmetric rank 2 tensor gxa with itself using gmet to produce symmetric rank 2 tensor.
INPUTS
gxa(2,6)=rank 2 complex tensor gmet(3,3)=real symmetric metric tensor (full storage)
OUTPUT
rank2(6)=rank 2 real tensor (symmetric storage)
NOTES
Symmetric gxa is stored as 11 22 33 32 31 21; gmet(3,3) is symmetric but stored fully (9 elements); output rank2 is stored as 11 22 33 32 31 21. Want $2 Re[contraction]$. $rank2(a,b)=2 Re[gxa(i,a)^"*" gmet(i,j) gxa(j,b)]$.
SOURCE
155 subroutine cont22(gxa,gmet,rank2) 156 157 !Arguments ------------------------------------ 158 !arrays 159 real(dp),intent(in) :: gmet(3,3),gxa(2,6) 160 real(dp),intent(out) :: rank2(6) 161 162 !Local variables------------------------------- 163 !scalars 164 integer,parameter :: im=2,re=1 165 166 ! ************************************************************************* 167 168 !Simply write out index summations 169 !a=1, b=1 in rank2(a,b) --> maps to index 1 170 rank2(1)=2.0d0*(& 171 & gmet(1,1)*(gxa(re,1)*gxa(re,1)+gxa(im,1)*gxa(im,1))+& 172 & gmet(2,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6))+& 173 & gmet(3,3)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+& 174 & 2.0d0*(& 175 & gmet(3,2)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6))+& 176 & gmet(3,1)*(gxa(re,5)*gxa(re,1)+gxa(im,5)*gxa(im,1))+& 177 & gmet(2,1)*(gxa(re,6)*gxa(re,1)+gxa(im,6)*gxa(im,1)))) 178 179 !a=2, b=2 in rank2(a,b) --> maps to index 2 180 rank2(2)=2.0d0*(& 181 & gmet(1,1)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6))+& 182 & gmet(2,2)*(gxa(re,2)*gxa(re,2)+gxa(im,2)*gxa(im,2))+& 183 & gmet(3,3)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+& 184 & 2.0d0*(& 185 & gmet(3,2)*(gxa(re,4)*gxa(re,2)+gxa(im,4)*gxa(im,2))+& 186 & gmet(3,1)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+& 187 & gmet(2,1)*(gxa(re,2)*gxa(re,6)+gxa(im,2)*gxa(im,6)))) 188 189 !a=3, b=3 in rank2(a,b) --> maps to index 3 190 rank2(3)=2.0d0*(& 191 & gmet(1,1)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+& 192 & gmet(2,2)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+& 193 & gmet(3,3)*(gxa(re,3)*gxa(re,3)+gxa(im,3)*gxa(im,3))+& 194 & 2.0d0*(& 195 & gmet(3,2)*(gxa(re,4)*gxa(re,3)+gxa(im,4)*gxa(im,3))+& 196 & gmet(3,1)*(gxa(re,5)*gxa(re,3)+gxa(im,5)*gxa(im,3))+& 197 & gmet(2,1)*(gxa(re,5)*gxa(re,4)+gxa(im,5)*gxa(im,4)))) 198 199 !a=3, b=2 in rank2(a,b) --> maps to index 4 200 rank2(4)=2.0d0*(& 201 & gmet(1,1)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6))+& 202 & gmet(2,2)*(gxa(re,4)*gxa(re,2)+gxa(im,4)*gxa(im,2))+& 203 & gmet(3,3)*(gxa(re,3)*gxa(re,4)+gxa(im,3)*gxa(im,4))+& 204 & gmet(3,2)*(gxa(re,3)*gxa(re,2)+gxa(im,3)*gxa(im,2))+& 205 & gmet(3,1)*(gxa(re,3)*gxa(re,6)+gxa(im,3)*gxa(im,6))+& 206 & gmet(2,1)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+& 207 & gmet(2,3)*(gxa(re,4)*gxa(re,4)+gxa(im,4)*gxa(im,4))+& 208 & gmet(1,3)*(gxa(re,5)*gxa(re,4)+gxa(im,5)*gxa(im,4))+& 209 & gmet(1,2)*(gxa(re,5)*gxa(re,2)+gxa(im,5)*gxa(im,2))) 210 211 !a=3, b=1 in rank2(a,b) --> maps to index 5 212 rank2(5)=2.0d0*(& 213 & gmet(1,1)*(gxa(re,5)*gxa(re,1)+gxa(im,5)*gxa(im,1))+& 214 & gmet(2,2)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+& 215 & gmet(3,3)*(gxa(re,3)*gxa(re,5)+gxa(im,3)*gxa(im,5))+& 216 & gmet(3,2)*(gxa(re,3)*gxa(re,6)+gxa(im,3)*gxa(im,6))+& 217 & gmet(3,1)*(gxa(re,3)*gxa(re,1)+gxa(im,3)*gxa(im,1))+& 218 & gmet(2,1)*(gxa(re,4)*gxa(re,1)+gxa(im,4)*gxa(im,1))+& 219 & gmet(2,3)*(gxa(re,4)*gxa(re,5)+gxa(im,4)*gxa(im,5))+& 220 & gmet(1,3)*(gxa(re,5)*gxa(re,5)+gxa(im,5)*gxa(im,5))+& 221 & gmet(1,2)*(gxa(re,5)*gxa(re,6)+gxa(im,5)*gxa(im,6))) 222 223 !a=2, b=1 in rank2(a,b) --> maps to index 6 224 rank2(6)=2.0d0*(& 225 & gmet(1,1)*(gxa(re,6)*gxa(re,1)+gxa(im,6)*gxa(im,1))+& 226 & gmet(2,2)*(gxa(re,2)*gxa(re,6)+gxa(im,2)*gxa(im,6))+& 227 & gmet(3,3)*(gxa(re,4)*gxa(re,5)+gxa(im,4)*gxa(im,5))+& 228 & gmet(3,2)*(gxa(re,4)*gxa(re,6)+gxa(im,4)*gxa(im,6))+& 229 & gmet(3,1)*(gxa(re,4)*gxa(re,1)+gxa(im,4)*gxa(im,1))+& 230 & gmet(2,1)*(gxa(re,2)*gxa(re,1)+gxa(im,2)*gxa(im,1))+& 231 & gmet(2,3)*(gxa(re,2)*gxa(re,5)+gxa(im,2)*gxa(im,5))+& 232 & gmet(1,3)*(gxa(re,6)*gxa(re,5)+gxa(im,6)*gxa(im,5))+& 233 & gmet(1,2)*(gxa(re,6)*gxa(re,6)+gxa(im,6)*gxa(im,6))) 234 235 end subroutine cont22
m_contract/cont22cso [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont22cso
FUNCTION
Contract symmetric rank 2 tensor gxa1 with symmetric rank 2 tensor gxa2 using metric tensor gmet to produce rank 2 complex tensor.
INPUTS
gxa1(2,10)=rank 2 complex symmetric tensor gxa2(2,10)=rank 2 complex symmetric tensor gmet(3,3)=usual metric tensor (symmetric, real)
OUTPUT
rank2c(2,6)=rank 2 complex tensor (pseudo-symmetric storage)
NOTES
This contraction is used for spin-orbit correction in non-local contribution to stresses. Symmetric gxa1, gxa2 are stored as 11 22 33 32 31 21; gmet(3,3) is symmetric but stored fully (9 elements); Output rank2c is not symmetric but since $rank2c_{gxa1,gxa2}(a,b)=conjg(rank2c_{gxa2,gxa1}(b,a))$ it is stored as 11 22 33 32 31 21. rank2c(1,1), rank2c(2,2), rank3c(3,3) are not needed; They are not calculated. rank2c(a,b)=3 conjg(gxa1(i,a)) gmet(i,j) gxa2(j,b)
SOURCE
273 subroutine cont22cso(gxa1,gxa2,gmet,rank2c) 274 275 !Arguments ------------------------------------ 276 !arrays 277 real(dp),intent(in) :: gmet(3,3),gxa1(2,6),gxa2(2,6) 278 real(dp),intent(out) :: rank2c(2,6) 279 280 !Local variables------------------------------- 281 !scalars 282 integer,parameter :: im=2,re=1 283 !arrays 284 real(dp) :: r2(2,3,3) 285 286 ! ************************************************************************* 287 288 !Initialize output tensor 289 rank2c(:,:)=0.d0 290 291 !First compute r2(i,j) = gmet(i,k) gxa2(j,k) 292 r2(:,1,1)=gmet(1,1)*gxa2(:,1)+gmet(1,2)*gxa2(:,6)+gmet(1,3)*gxa2(:,5) 293 r2(:,1,2)=gmet(1,1)*gxa2(:,6)+gmet(1,2)*gxa2(:,2)+gmet(1,3)*gxa2(:,4) 294 r2(:,1,3)=gmet(1,1)*gxa2(:,5)+gmet(1,2)*gxa2(:,4)+gmet(1,3)*gxa2(:,3) 295 r2(:,2,1)=gmet(2,1)*gxa2(:,1)+gmet(2,2)*gxa2(:,6)+gmet(2,3)*gxa2(:,5) 296 r2(:,2,2)=gmet(2,1)*gxa2(:,6)+gmet(2,2)*gxa2(:,2)+gmet(2,3)*gxa2(:,4) 297 r2(:,2,3)=gmet(2,1)*gxa2(:,5)+gmet(2,2)*gxa2(:,4)+gmet(2,3)*gxa2(:,3) 298 r2(:,3,1)=gmet(3,1)*gxa2(:,1)+gmet(3,2)*gxa2(:,6)+gmet(3,3)*gxa2(:,5) 299 r2(:,3,2)=gmet(3,1)*gxa2(:,6)+gmet(3,2)*gxa2(:,2)+gmet(3,3)*gxa2(:,4) 300 r2(:,3,3)=gmet(3,1)*gxa2(:,5)+gmet(3,2)*gxa2(:,4)+gmet(3,3)*gxa2(:,3) 301 302 !Then compute rank2c(a,b) = 3 conjg(gxa1(a,i)) r2(i,b) 303 !stored as 11 22 33 32 31 21 304 !rank2c(re,1) = 3.d0*(gxa1(re,1)*r2(re,1,1)+gxa1(im,1)*r2(im,1,1)& 305 !& +gxa1(re,6)*r2(re,2,1)+gxa1(im,6)*r2(im,2,1)& 306 !& +gxa1(re,5)*r2(re,3,1)+gxa1(im,5)*r2(im,3,1)) 307 !rank2c(re,2) = 3.d0*(gxa1(re,6)*r2(re,1,2)+gxa1(im,6)*r2(im,1,2)& 308 !& +gxa1(re,2)*r2(re,2,2)+gxa1(im,2)*r2(im,2,2)& 309 !& +gxa1(re,4)*r2(re,3,2)+gxa1(im,4)*r2(im,3,2)) 310 !rank2c(re,3) = 3.d0*(gxa1(re,5)*r2(re,1,3)+gxa1(im,5)*r2(im,1,3)& 311 !& +gxa1(re,4)*r2(re,2,3)+gxa1(im,4)*r2(im,2,3)& 312 !& +gxa1(re,3)*r2(re,3,3)+gxa1(im,3)*r2(im,3,3)) 313 rank2c(re,4) = 3.d0*(gxa1(re,5)*r2(re,1,2)+gxa1(im,5)*r2(im,1,2)& 314 & +gxa1(re,4)*r2(re,2,2)+gxa1(im,4)*r2(im,2,2)& 315 & +gxa1(re,3)*r2(re,3,2)+gxa1(im,3)*r2(im,3,2)) 316 rank2c(re,5) = 3.d0*(gxa1(re,5)*r2(re,1,1)+gxa1(im,5)*r2(im,1,1)& 317 & +gxa1(re,4)*r2(re,2,1)+gxa1(im,4)*r2(im,2,1)& 318 & +gxa1(re,3)*r2(re,3,1)+gxa1(im,3)*r2(im,3,1)) 319 rank2c(re,6) = 3.d0*(gxa1(re,6)*r2(re,1,1)+gxa1(im,6)*r2(im,1,1)& 320 & +gxa1(re,2)*r2(re,2,1)+gxa1(im,2)*r2(im,2,1)& 321 & +gxa1(re,4)*r2(re,3,1)+gxa1(im,4)*r2(im,3,1)) 322 !rank2c(im,1) = 3.d0*(gxa1(re,1)*r2(im,1,1)-gxa1(im,1)*r2(re,1,1)& 323 !& +gxa1(re,6)*r2(im,2,1)-gxa1(im,6)*r2(re,2,1)& 324 !& +gxa1(re,5)*r2(im,3,1)-gxa1(im,5)*r2(re,3,1)) 325 !rank2c(im,2) = 3.d0*(gxa1(re,6)*r2(im,1,2)-gxa1(im,6)*r2(re,1,2)& 326 !& +gxa1(re,2)*r2(im,2,2)-gxa1(im,2)*r2(re,2,2)& 327 !& +gxa1(re,4)*r2(im,3,2)-gxa1(im,4)*r2(re,3,2)) 328 !rank2c(im,3) = 3.d0*(gxa1(re,5)*r2(im,1,3)-gxa1(im,5)*r2(re,1,3)& 329 !& +gxa1(re,4)*r2(im,2,3)-gxa1(im,4)*r2(re,2,3)& 330 !& +gxa1(re,3)*r2(im,3,3)-gxa1(im,3)*r2(re,3,3)) 331 rank2c(im,4) = 3.d0*(gxa1(re,5)*r2(im,1,2)-gxa1(im,5)*r2(re,1,2)& 332 & +gxa1(re,4)*r2(im,2,2)-gxa1(im,4)*r2(re,2,2)& 333 & +gxa1(re,3)*r2(im,3,2)-gxa1(im,3)*r2(re,3,2)) 334 rank2c(im,5) = 3.d0*(gxa1(re,5)*r2(im,1,1)-gxa1(im,5)*r2(re,1,1)& 335 & +gxa1(re,4)*r2(im,2,1)-gxa1(im,4)*r2(re,2,1)& 336 & +gxa1(re,3)*r2(im,3,1)-gxa1(im,3)*r2(re,3,1)) 337 rank2c(im,6) = 3.d0*(gxa1(re,6)*r2(im,1,1)-gxa1(im,6)*r2(re,1,1)& 338 & +gxa1(re,2)*r2(im,2,1)-gxa1(im,2)*r2(re,2,1)& 339 & +gxa1(re,4)*r2(im,3,1)-gxa1(im,4)*r2(re,3,1)) 340 341 end subroutine cont22cso
m_contract/cont22so [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont22so
FUNCTION
Contract symmetric rank 2 tensor gxa1 with symmetric rank 2 tensor gxa2 using antisymmetric tensor amet to produce rank 2 real tensor.
INPUTS
gxa1(2,6)=rank 2 complex symmetric tensor gxa2(2,6)=rank 2 complex symmetric tensor amet(2,3,3)=antisymmetric complex tensor used for spin-orbit
OUTPUT
rank2(6)=rank 2 real tensor (pseudo-symmetric storage)
NOTES
This contraction is used for spin-orbit correction in non-local contribution to stresses. Symmetric gxa1, gxa2 are stored as 11 22 33 32 31 21; amet(3,3) is antisymmetric but stored fully (9 elements); Output rank2 is not symmetric but since $rank2_{gxa1,gxa2}(a,b)=conjg(rank2_{gxa2,gxa1}(b,a))$ it is stored as 11 22 33 32 31 21. Want 2*Re[contraction]. rank2(a,b)=2 Re[conjg(gxa1(i,a)) amet(i,j) gxa2(j,b)] Note that, since amet is antisymmetric, amet(i,i)=0
SOURCE
380 subroutine cont22so(gxa1,gxa2,amet,rank2) 381 382 !Arguments ------------------------------------ 383 !arrays 384 real(dp),intent(in) :: amet(2,3,3),gxa1(2,6),gxa2(2,6) 385 real(dp),intent(out) :: rank2(6) 386 387 !Local variables------------------------------- 388 !scalars 389 integer,parameter :: im=2,re=1 390 391 ! ************************************************************************* 392 393 !Simply write out index summations 394 !a=1, b=1 in rank2(a,b) --> maps to index 1 395 rank2(1)=2.0d0*(& 396 & amet(1,3,2)*(gxa1(re,5)*gxa2(re,6)+gxa1(im,5)*gxa2(im,6))-& 397 & amet(2,3,2)*(gxa1(re,5)*gxa2(im,6)-gxa1(im,5)*gxa2(re,6))+& 398 & amet(1,3,1)*(gxa1(re,5)*gxa2(re,1)+gxa1(im,5)*gxa2(im,1))-& 399 & amet(2,3,1)*(gxa1(re,5)*gxa2(im,1)-gxa1(im,5)*gxa2(re,1))+& 400 & amet(1,2,1)*(gxa1(re,6)*gxa2(re,1)+gxa1(im,6)*gxa2(im,1))-& 401 & amet(2,2,1)*(gxa1(re,6)*gxa2(im,1)-gxa1(im,6)*gxa2(re,1))+& 402 & amet(1,2,3)*(gxa1(re,6)*gxa2(re,5)+gxa1(im,6)*gxa2(im,5))-& 403 & amet(2,2,3)*(gxa1(re,6)*gxa2(im,5)-gxa1(im,6)*gxa2(re,5))+& 404 & amet(1,1,3)*(gxa1(re,1)*gxa2(re,5)+gxa1(im,1)*gxa2(im,5))-& 405 & amet(2,1,3)*(gxa1(re,1)*gxa2(im,5)-gxa1(im,1)*gxa2(re,5))+& 406 & amet(1,1,2)*(gxa1(re,1)*gxa2(re,6)+gxa1(im,1)*gxa2(im,6))-& 407 & amet(2,1,2)*(gxa1(re,1)*gxa2(im,6)-gxa1(im,1)*gxa2(re,6))) 408 409 !a=2, b=2 in rank2(a,b) --> maps to index 2 410 rank2(2)=2.0d0*(& 411 & amet(1,3,2)*(gxa1(re,4)*gxa2(re,2)+gxa1(im,4)*gxa2(im,2))-& 412 & amet(2,3,2)*(gxa1(re,4)*gxa2(im,2)-gxa1(im,4)*gxa2(re,2))+& 413 & amet(1,3,1)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-& 414 & amet(2,3,1)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+& 415 & amet(1,2,1)*(gxa1(re,2)*gxa2(re,6)+gxa1(im,2)*gxa2(im,6))-& 416 & amet(2,2,1)*(gxa1(re,2)*gxa2(im,6)-gxa1(im,2)*gxa2(re,6))+& 417 & amet(1,2,3)*(gxa1(re,2)*gxa2(re,4)+gxa1(im,2)*gxa2(im,4))-& 418 & amet(2,2,3)*(gxa1(re,2)*gxa2(im,4)-gxa1(im,2)*gxa2(re,4))+& 419 & amet(1,1,3)*(gxa1(re,6)*gxa2(re,4)+gxa1(im,6)*gxa2(im,4))-& 420 & amet(2,1,3)*(gxa1(re,6)*gxa2(im,4)-gxa1(im,6)*gxa2(re,4))+& 421 & amet(1,1,2)*(gxa1(re,6)*gxa2(re,2)+gxa1(im,6)*gxa2(im,2))-& 422 & amet(2,1,2)*(gxa1(re,6)*gxa2(im,2)-gxa1(im,6)*gxa2(re,2))) 423 424 !a=3, b=3 in rank2(a,b) --> maps to index 3 425 rank2(3)=2.0d0*(& 426 & amet(1,3,2)*(gxa1(re,3)*gxa2(re,4)+gxa1(im,3)*gxa2(im,4))-& 427 & amet(2,3,2)*(gxa1(re,3)*gxa2(im,4)-gxa1(im,3)*gxa2(re,4))+& 428 & amet(1,3,1)*(gxa1(re,3)*gxa2(re,5)+gxa1(im,3)*gxa2(im,5))-& 429 & amet(2,3,1)*(gxa1(re,3)*gxa2(im,5)-gxa1(im,3)*gxa2(re,5))+& 430 & amet(1,2,1)*(gxa1(re,4)*gxa2(re,5)+gxa1(im,4)*gxa2(im,5))-& 431 & amet(2,2,1)*(gxa1(re,4)*gxa2(im,5)-gxa1(im,4)*gxa2(re,5))+& 432 & amet(1,2,3)*(gxa1(re,4)*gxa2(re,3)+gxa1(im,4)*gxa2(im,3))-& 433 & amet(2,2,3)*(gxa1(re,4)*gxa2(im,3)-gxa1(im,4)*gxa2(re,3))+& 434 & amet(1,1,3)*(gxa1(re,5)*gxa2(re,3)+gxa1(im,5)*gxa2(im,3))-& 435 & amet(2,1,3)*(gxa1(re,5)*gxa2(im,3)-gxa1(im,5)*gxa2(re,3))+& 436 & amet(1,1,2)*(gxa1(re,5)*gxa2(re,4)+gxa1(im,5)*gxa2(im,4))-& 437 & amet(2,1,2)*(gxa1(re,5)*gxa2(im,4)-gxa1(im,5)*gxa2(re,4))) 438 439 !a=3, b=2 in rank2(a,b) --> maps to index 4 440 rank2(4)=2.0d0*(& 441 & amet(1,3,2)*(gxa1(re,3)*gxa2(re,2)+gxa1(im,3)*gxa2(im,2))-& 442 & amet(2,3,2)*(gxa1(re,3)*gxa2(im,2)-gxa1(im,3)*gxa2(re,2))+& 443 & amet(1,3,1)*(gxa1(re,3)*gxa2(re,6)+gxa1(im,3)*gxa2(im,6))-& 444 & amet(2,3,1)*(gxa1(re,3)*gxa2(im,6)-gxa1(im,3)*gxa2(re,6))+& 445 & amet(1,2,1)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-& 446 & amet(2,2,1)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+& 447 & amet(1,2,3)*(gxa1(re,4)*gxa2(re,4)+gxa1(im,4)*gxa2(im,4))-& 448 & amet(2,2,3)*(gxa1(re,4)*gxa2(im,4)-gxa1(im,4)*gxa2(re,4))+& 449 & amet(1,1,3)*(gxa1(re,5)*gxa2(re,4)+gxa1(im,5)*gxa2(im,4))-& 450 & amet(2,1,3)*(gxa1(re,5)*gxa2(im,4)-gxa1(im,5)*gxa2(re,4))+& 451 & amet(1,1,2)*(gxa1(re,5)*gxa2(re,2)+gxa1(im,5)*gxa2(im,2))-& 452 & amet(2,1,2)*(gxa1(re,5)*gxa2(im,2)-gxa1(im,5)*gxa2(re,2))) 453 454 !a=3, b=1 in rank2(a,b) --> maps to index 5 455 rank2(5)=2.0d0*(& 456 & amet(1,3,2)*(gxa1(re,3)*gxa2(re,6)+gxa1(im,3)*gxa2(im,6))-& 457 & amet(2,3,2)*(gxa1(re,3)*gxa2(im,6)-gxa1(im,3)*gxa2(re,6))+& 458 & amet(1,3,1)*(gxa1(re,3)*gxa2(re,1)+gxa1(im,3)*gxa2(im,1))-& 459 & amet(2,3,1)*(gxa1(re,3)*gxa2(im,1)-gxa1(im,3)*gxa2(re,1))+& 460 & amet(1,2,1)*(gxa1(re,4)*gxa2(re,1)+gxa1(im,4)*gxa2(im,1))-& 461 & amet(2,2,1)*(gxa1(re,4)*gxa2(im,1)-gxa1(im,4)*gxa2(re,1))+& 462 & amet(1,2,3)*(gxa1(re,4)*gxa2(re,5)+gxa1(im,4)*gxa2(im,5))-& 463 & amet(2,2,3)*(gxa1(re,4)*gxa2(im,5)-gxa1(im,4)*gxa2(re,5))+& 464 & amet(1,1,3)*(gxa1(re,5)*gxa2(re,5)+gxa1(im,5)*gxa2(im,5))-& 465 & amet(2,1,3)*(gxa1(re,5)*gxa2(im,5)-gxa1(im,5)*gxa2(re,5))+& 466 & amet(1,1,2)*(gxa1(re,5)*gxa2(re,6)+gxa1(im,5)*gxa2(im,6))-& 467 & amet(2,1,2)*(gxa1(re,5)*gxa2(im,6)-gxa1(im,5)*gxa2(re,6))) 468 469 !a=2, b=1 in rank2(a,b) --> maps to index 6 470 rank2(6)=2.0d0*(& 471 & amet(1,3,2)*(gxa1(re,4)*gxa2(re,6)+gxa1(im,4)*gxa2(im,6))-& 472 & amet(2,3,2)*(gxa1(re,4)*gxa2(im,6)-gxa1(im,4)*gxa2(re,6))+& 473 & amet(1,3,1)*(gxa1(re,4)*gxa2(re,1)+gxa1(im,4)*gxa2(im,1))-& 474 & amet(2,3,1)*(gxa1(re,4)*gxa2(im,1)-gxa1(im,4)*gxa2(re,1))+& 475 & amet(1,2,1)*(gxa1(re,2)*gxa2(re,1)+gxa1(im,2)*gxa2(im,1))-& 476 & amet(2,2,1)*(gxa1(re,2)*gxa2(im,1)-gxa1(im,2)*gxa2(re,1))+& 477 & amet(1,2,3)*(gxa1(re,2)*gxa2(re,5)+gxa1(im,2)*gxa2(im,5))-& 478 & amet(2,2,3)*(gxa1(re,2)*gxa2(im,5)-gxa1(im,2)*gxa2(re,5))+& 479 & amet(1,1,3)*(gxa1(re,6)*gxa2(re,5)+gxa1(im,6)*gxa2(im,5))-& 480 & amet(2,1,3)*(gxa1(re,6)*gxa2(im,5)-gxa1(im,6)*gxa2(re,5))+& 481 & amet(1,1,2)*(gxa1(re,6)*gxa2(re,6)+gxa1(im,6)*gxa2(im,6))-& 482 & amet(2,1,2)*(gxa1(re,6)*gxa2(im,6)-gxa1(im,6)*gxa2(re,6))) 483 484 end subroutine cont22so
m_contract/cont24 [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont24
FUNCTION
Contract symmetric rank2 tensor gxa with rank4 symmetric tensor to produce symmetric rank2 tensor.
INPUTS
gxa(2,6)=rank 2 symmetric complex tensor in order 11 22 33 32 31 21 rank4(2,15)=rank 4 complex tensor (symmetric storage)
OUTPUT
rank2(6)=rank 2 real tensor (symmetric storage) 11 22 33 32 31 21.
NOTES
Tensors are in "symmetric" storage mode. for gxa and rank2 this is 11, 22, 33, 32, 31, 21; for the rank 4 tensor rank4 this is 1111 2211 3311 3211 3111 2111 2221 3321 3221 3331 2222 3322 3222 3332 3333. gxa and rank4 are complex; rank2 is real. Want $2 Re[contraction]$. $rank2(a,b)=2 Re[gxa(i,j)^"*" rank4(a,b,i,j)]$. Note that the input gxa is typically the result of $gxa(i,j)=[{3 \over 2} gmet(i,l) gmet(j,m) - {1 \over 2} gmet(i,j) gmet(l,m)] gxa_old(l,m)$ where the subroutine "metcon" already includes weights in the definition of gxa for off-diagonal elements (weight of 2 for symmetry). Components 4, 5, and 6 of gxa have already been multiplied by 2 so the expressions below do not carry the 2.
SOURCE
522 subroutine cont24(gxa,rank4,rank2) 523 524 !Arguments ------------------------------------ 525 !arrays 526 real(dp),intent(in) :: gxa(2,6),rank4(2,15) 527 real(dp),intent(out) :: rank2(6) 528 529 !Local variables------------------------------- 530 !scalars 531 integer,parameter :: im=2,re=1 532 533 ! ************************************************************************* 534 535 !Simply write out index summations 536 537 !a=1, b=1 in rank2(a,b) --> maps to index 1 538 rank2(1)=2.0d0*(& 539 & gxa(re,1)*rank4(re, 1)+gxa(im,1)*rank4(im, 1)+& 540 & gxa(re,2)*rank4(re, 2)+gxa(im,2)*rank4(im, 2)+& 541 & gxa(re,3)*rank4(re, 3)+gxa(im,3)*rank4(im, 3)+& 542 & gxa(re,4)*rank4(re, 4)+gxa(im,4)*rank4(im, 4)+& 543 & gxa(re,5)*rank4(re, 5)+gxa(im,5)*rank4(im, 5)+& 544 & gxa(re,6)*rank4(re, 6)+gxa(im,6)*rank4(im, 6)) 545 546 !a=2, b=2 in rank2(a,b) --> maps to index 2 547 rank2(2)=2.0d0*(& 548 & gxa(re,1)*rank4(re, 2)+gxa(im,1)*rank4(im, 2)+& 549 & gxa(re,2)*rank4(re,11)+gxa(im,2)*rank4(im,11)+& 550 & gxa(re,3)*rank4(re,12)+gxa(im,3)*rank4(im,12)+& 551 & gxa(re,4)*rank4(re,13)+gxa(im,4)*rank4(im,13)+& 552 & gxa(re,5)*rank4(re, 9)+gxa(im,5)*rank4(im, 9)+& 553 & gxa(re,6)*rank4(re, 7)+gxa(im,6)*rank4(im, 7)) 554 555 !a=3, b=3 in rank2(a,b) --> maps to index 3 556 rank2(3)=2.0d0*(& 557 & gxa(re,1)*rank4(re, 3)+gxa(im,1)*rank4(im, 3)+& 558 & gxa(re,2)*rank4(re,12)+gxa(im,2)*rank4(im,12)+& 559 & gxa(re,3)*rank4(re,15)+gxa(im,3)*rank4(im,15)+& 560 & gxa(re,4)*rank4(re,14)+gxa(im,4)*rank4(im,14)+& 561 & gxa(re,5)*rank4(re,10)+gxa(im,5)*rank4(im,10)+& 562 & gxa(re,6)*rank4(re, 8)+gxa(im,6)*rank4(im, 8)) 563 564 !a=3, b=2 in rank2(a,b) --> maps to index 4 565 rank2(4)=2.0d0*(& 566 & gxa(re,1)*rank4(re, 4)+gxa(im,1)*rank4(im, 4)+& 567 & gxa(re,2)*rank4(re,13)+gxa(im,2)*rank4(im,13)+& 568 & gxa(re,3)*rank4(re,14)+gxa(im,3)*rank4(im,14)+& 569 & gxa(re,4)*rank4(re,12)+gxa(im,4)*rank4(im,12)+& 570 & gxa(re,5)*rank4(re, 8)+gxa(im,5)*rank4(im, 8)+& 571 & gxa(re,6)*rank4(re, 9)+gxa(im,6)*rank4(im, 9)) 572 573 !a=3, b=1 in rank2(a,b) --> maps to index 5 574 rank2(5)=2.0d0*(& 575 & gxa(re,1)*rank4(re, 5)+gxa(im,1)*rank4(im, 5)+& 576 & gxa(re,2)*rank4(re, 9)+gxa(im,2)*rank4(im, 9)+& 577 & gxa(re,3)*rank4(re,10)+gxa(im,3)*rank4(im,10)+& 578 & gxa(re,4)*rank4(re, 8)+gxa(im,4)*rank4(im, 8)+& 579 & gxa(re,5)*rank4(re, 3)+gxa(im,5)*rank4(im, 3)+& 580 & gxa(re,6)*rank4(re, 4)+gxa(im,6)*rank4(im, 4)) 581 582 !a=2, b=1 in rank2(a,b) --> maps to index 6 583 rank2(6)=2.0d0*(& 584 & gxa(re,1)*rank4(re, 6)+gxa(im,1)*rank4(im, 6)+& 585 & gxa(re,2)*rank4(re, 7)+gxa(im,2)*rank4(im, 7)+& 586 & gxa(re,3)*rank4(re, 8)+gxa(im,3)*rank4(im, 8)+& 587 & gxa(re,4)*rank4(re, 9)+gxa(im,4)*rank4(im, 9)+& 588 & gxa(re,5)*rank4(re, 4)+gxa(im,5)*rank4(im, 4)+& 589 & gxa(re,6)*rank4(re, 2)+gxa(im,6)*rank4(im, 2)) 590 591 end subroutine cont24
m_contract/cont3 [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont3
FUNCTION
Compute several specialized contractions needed for the l=3 part of the stress tensor.
INPUTS
gxa(2,10)=complex symmetric rank 3 tensor gmet(3,3)=usual metric tensor, a symmetric matrix stored in full storage mode (bohr^-2)
OUTPUT
rank2(6)=2*Re[contraction] given by 2*Re[(15/2)*r3(a,i,j)*r3(b,j,i)-3*r1(i)*r3(a,b,i)-(3/2)*r1(a)*r1(b)] where r3(a,i,j)=gmet(j,k) gxa(a,i,k) and r1(a)=gmet(i,j) gxa(i,j,a). rank2 is stored in the compressed form 11 22 33 32 31 21.
NOTES
Input gxa is a completely symmetric rank 3 tensor (complex) in compressed storage: 111 221 331 321 311 211 222 332 322 333. The output tensor is completely symmetric rank 2, real, and is given by $2 Re[{15 \over 2} r3(a,i,j) r3(b,j,i) - 3 r1(i) r3(a,b,i) - {3 \over 2} r1(a) r1(b)]$ where $r3(a,i,j)=gmet(j,k) gxa(a,i,k)$ and $r1(a)=gmet(i,j) gxa(i,j,a)$. rank2 is stored in the compressed form 11 22 33 32 31 21.
SOURCE
623 subroutine cont3(gxa,gmet,rank2) 624 625 !Arguments ------------------------------------ 626 !arrays 627 real(dp),intent(in) :: gmet(3,3),gxa(2,10) 628 real(dp),intent(out) :: rank2(6) 629 630 !Local variables------------------------------- 631 !scalars 632 integer,parameter :: im=2,re=1 633 integer :: ii 634 !arrays 635 real(dp) :: r1(2,3),r3(2,18),s13(6),s33(6) 636 637 ! ************************************************************************* 638 639 !Compute r1(a) = gmet(i,j) gxa(i,j,a) 640 641 !Write out components for 3 distinct terms, Re and Im 642 do ii=1,2 643 r1(ii,1)=gmet(1,1)*gxa(ii,1)+gmet(2,2)*gxa(ii,2)+& 644 & gmet(3,3)*gxa(ii,3)+2.d0*(& 645 & gmet(3,2)*gxa(ii,4)+gmet(3,1)*gxa(ii,5)+& 646 & gmet(2,1)*gxa(ii,6)) 647 r1(ii,2)=gmet(1,1)*gxa(ii,6)+gmet(2,2)*gxa(ii,7)+& 648 & gmet(3,3)*gxa(ii,8)+2.d0*(& 649 & gmet(3,2)*gxa(ii,9)+gmet(3,1)*gxa(ii,4)+& 650 & gmet(2,1)*gxa(ii,2)) 651 r1(ii,3)=gmet(1,1)*gxa(ii,5)+gmet(2,2)*gxa(ii,9)+& 652 & gmet(3,3)*gxa(ii,10)+2.d0*(& 653 & gmet(3,2)*gxa(ii,8)+gmet(3,1)*gxa(ii,3)+& 654 & gmet(2,1)*gxa(ii,4)) 655 end do 656 657 !Compute r3(a,b,k)=gmet(k,n) gxa(a,b,n) 658 659 !Write out components for 18 distinct terms, Re and Im 660 !(symmetric in first two indices, not in all permutations) 661 !store as 111 221 331 321 311 211 662 !112 222 332 322 312 212 663 !113 223 333 323 313 213 664 do ii=1,2 665 r3(ii, 1)=gmet(1,1)*gxa(ii,1)+gmet(2,1)*gxa(ii,6)+& 666 & gmet(3,1)*gxa(ii,5) 667 r3(ii, 2)=gmet(1,1)*gxa(ii,2)+gmet(2,1)*gxa(ii,7)+& 668 & gmet(3,1)*gxa(ii,9) 669 r3(ii, 3)=gmet(1,1)*gxa(ii,3)+gmet(2,1)*gxa(ii,8)+& 670 & gmet(3,1)*gxa(ii,10) 671 r3(ii, 4)=gmet(1,1)*gxa(ii,4)+gmet(2,1)*gxa(ii,9)+& 672 & gmet(3,1)*gxa(ii,8) 673 r3(ii, 5)=gmet(1,1)*gxa(ii,5)+gmet(2,1)*gxa(ii,4)+& 674 & gmet(3,1)*gxa(ii,3) 675 r3(ii, 6)=gmet(1,1)*gxa(ii,6)+gmet(2,1)*gxa(ii,2)+& 676 & gmet(3,1)*gxa(ii,4) 677 r3(ii, 7)=gmet(2,1)*gxa(ii,1)+gmet(2,2)*gxa(ii,6)+& 678 & gmet(3,2)*gxa(ii,5) 679 r3(ii, 8)=gmet(2,1)*gxa(ii,2)+gmet(2,2)*gxa(ii,7)+& 680 & gmet(3,2)*gxa(ii,9) 681 r3(ii, 9)=gmet(2,1)*gxa(ii,3)+gmet(2,2)*gxa(ii,8)+& 682 & gmet(3,2)*gxa(ii,10) 683 r3(ii,10)=gmet(2,1)*gxa(ii,4)+gmet(2,2)*gxa(ii,9)+& 684 & gmet(3,2)*gxa(ii,8) 685 r3(ii,11)=gmet(2,1)*gxa(ii,5)+gmet(2,2)*gxa(ii,4)+& 686 & gmet(3,2)*gxa(ii,3) 687 r3(ii,12)=gmet(2,1)*gxa(ii,6)+gmet(2,2)*gxa(ii,2)+& 688 & gmet(3,2)*gxa(ii,4) 689 r3(ii,13)=gmet(3,1)*gxa(ii,1)+gmet(3,2)*gxa(ii,6)+& 690 & gmet(3,3)*gxa(ii,5) 691 r3(ii,14)=gmet(3,1)*gxa(ii,2)+gmet(3,2)*gxa(ii,7)+& 692 & gmet(3,3)*gxa(ii,9) 693 r3(ii,15)=gmet(3,1)*gxa(ii,3)+gmet(3,2)*gxa(ii,8)+& 694 & gmet(3,3)*gxa(ii,10) 695 r3(ii,16)=gmet(3,1)*gxa(ii,4)+gmet(3,2)*gxa(ii,9)+& 696 & gmet(3,3)*gxa(ii,8) 697 r3(ii,17)=gmet(3,1)*gxa(ii,5)+gmet(3,2)*gxa(ii,4)+& 698 & gmet(3,3)*gxa(ii,3) 699 r3(ii,18)=gmet(3,1)*gxa(ii,6)+gmet(3,2)*gxa(ii,2)+& 700 & gmet(3,3)*gxa(ii,4) 701 702 end do 703 704 !Now need 705 !2*Re[(15/2)*r3(a,i,j)*r3(b,j,i)-3*r1(i)*r3(a,b,i)-(3/2)*r1(a)*r1(b)]. 706 707 !Write out s33(a,b)=2*Re[r3(a,i,j)*r3(b,j,i)] 708 709 s33(1)=2.d0*(r3(re, 1)*r3(re, 1)+r3(im, 1)*r3(im, 1)+& 710 & r3(re,12)*r3(re,12)+r3(im,12)*r3(im,12)+& 711 & r3(re,17)*r3(re,17)+r3(im,17)*r3(im,17)+& 712 & r3(re,11)*r3(re,18)+r3(im,11)*r3(im,18)+& 713 & r3(re,18)*r3(re,11)+r3(im,18)*r3(im,11)+& 714 & r3(re, 5)*r3(re,13)+r3(im, 5)*r3(im,13)+& 715 & r3(re,13)*r3(re, 5)+r3(im,13)*r3(im, 5)+& 716 & r3(re, 6)*r3(re, 7)+r3(im, 6)*r3(im, 7)+& 717 & r3(re, 7)*r3(re, 6)+r3(im, 7)*r3(im, 6)) 718 719 s33(2)=2.d0*(r3(re, 6)*r3(re, 6)+r3(im, 6)*r3(im, 6)+& 720 & r3(re, 8)*r3(re, 8)+r3(im, 8)*r3(im, 8)+& 721 & r3(re,16)*r3(re,16)+r3(im,16)*r3(im,16)+& 722 & r3(re,10)*r3(re,14)+r3(im,10)*r3(im,14)+& 723 & r3(re,14)*r3(re,10)+r3(im,14)*r3(im,10)+& 724 & r3(re, 4)*r3(re,18)+r3(im, 4)*r3(im,18)+& 725 & r3(re,18)*r3(re, 4)+r3(im,18)*r3(im, 4)+& 726 & r3(re, 2)*r3(re,12)+r3(im, 2)*r3(im,12)+& 727 & r3(re,12)*r3(re, 2)+r3(im,12)*r3(im, 2)) 728 729 s33(3)=2.d0*(r3(re, 5)*r3(re, 5)+r3(im, 5)*r3(im, 5)+& 730 & r3(re,10)*r3(re,10)+r3(im,10)*r3(im,10)+& 731 & r3(re,15)*r3(re,15)+r3(im,15)*r3(im,15)+& 732 & r3(re, 9)*r3(re,16)+r3(im, 9)*r3(im,16)+& 733 & r3(re,16)*r3(re, 9)+r3(im,16)*r3(im, 9)+& 734 & r3(re, 3)*r3(re,17)+r3(im, 3)*r3(im,17)+& 735 & r3(re,17)*r3(re, 3)+r3(im,17)*r3(im, 3)+& 736 & r3(re, 4)*r3(re,11)+r3(im, 4)*r3(im,11)+& 737 & r3(re,11)*r3(re, 4)+r3(im,11)*r3(im, 4)) 738 739 s33(4)=2.d0*(r3(re, 5)*r3(re, 6)+r3(im, 5)*r3(im, 6)+& 740 & r3(re,10)*r3(re, 8)+r3(im,10)*r3(im, 8)+& 741 & r3(re,15)*r3(re,16)+r3(im,15)*r3(im,16)+& 742 & r3(re, 9)*r3(re,14)+r3(im, 9)*r3(im,14)+& 743 & r3(re,16)*r3(re,10)+r3(im,16)*r3(im,10)+& 744 & r3(re, 3)*r3(re,18)+r3(im, 3)*r3(im,18)+& 745 & r3(re,17)*r3(re, 4)+r3(im,17)*r3(im, 4)+& 746 & r3(re, 4)*r3(re,12)+r3(im, 4)*r3(im,12)+& 747 & r3(re,11)*r3(re, 2)+r3(im,11)*r3(im, 2)) 748 749 s33(5)=2.d0*(r3(re, 5)*r3(re, 1)+r3(im, 5)*r3(im, 1)+& 750 & r3(re,10)*r3(re,12)+r3(im,10)*r3(im,12)+& 751 & r3(re,15)*r3(re,17)+r3(im,15)*r3(im,17)+& 752 & r3(re, 9)*r3(re,18)+r3(im, 9)*r3(im,18)+& 753 & r3(re,16)*r3(re,11)+r3(im,16)*r3(im,11)+& 754 & r3(re, 3)*r3(re,13)+r3(im, 3)*r3(im,13)+& 755 & r3(re,17)*r3(re, 5)+r3(im,17)*r3(im, 5)+& 756 & r3(re, 4)*r3(re, 7)+r3(im, 4)*r3(im, 7)+& 757 & r3(re,11)*r3(re, 6)+r3(im,11)*r3(im, 6)) 758 759 s33(6)=2.d0*(r3(re, 6)*r3(re, 1)+r3(im, 6)*r3(im, 1)+& 760 & r3(re, 8)*r3(re,12)+r3(im, 8)*r3(im,12)+& 761 & r3(re,16)*r3(re,17)+r3(im,16)*r3(im,17)+& 762 & r3(re,10)*r3(re,18)+r3(im,10)*r3(im,18)+& 763 & r3(re,14)*r3(re,11)+r3(im,14)*r3(im,11)+& 764 & r3(re, 4)*r3(re,13)+r3(im, 4)*r3(im,13)+& 765 & r3(re,18)*r3(re, 5)+r3(im,18)*r3(im, 5)+& 766 & r3(re, 2)*r3(re, 7)+r3(im, 2)*r3(im, 7)+& 767 & r3(re,12)*r3(re, 6)+r3(im,12)*r3(im, 6)) 768 769 770 !Write out s13(a,b)=2*Re[r1(i)*r3(a,b,i)] 771 772 s13(1)=2.d0*(r1(re,1)*r3(re, 1)+r1(im,1)*r3(im, 1)+& 773 & r1(re,2)*r3(re, 7)+r1(im,2)*r3(im, 7)+& 774 & r1(re,3)*r3(re,13)+r1(im,3)*r3(im,13)) 775 s13(2)=2.d0*(r1(re,1)*r3(re, 2)+r1(im,1)*r3(im, 2)+& 776 & r1(re,2)*r3(re, 8)+r1(im,2)*r3(im, 8)+& 777 & r1(re,3)*r3(re,14)+r1(im,3)*r3(im,14)) 778 s13(3)=2.d0*(r1(re,1)*r3(re, 3)+r1(im,1)*r3(im, 3)+& 779 & r1(re,2)*r3(re, 9)+r1(im,2)*r3(im, 9)+& 780 & r1(re,3)*r3(re,15)+r1(im,3)*r3(im,15)) 781 s13(4)=2.d0*(r1(re,1)*r3(re, 4)+r1(im,1)*r3(im, 4)+& 782 & r1(re,2)*r3(re,10)+r1(im,2)*r3(im,10)+& 783 & r1(re,3)*r3(re,16)+r1(im,3)*r3(im,16)) 784 s13(5)=2.d0*(r1(re,1)*r3(re, 5)+r1(im,1)*r3(im, 5)+& 785 & r1(re,2)*r3(re,11)+r1(im,2)*r3(im,11)+& 786 & r1(re,3)*r3(re,17)+r1(im,3)*r3(im,17)) 787 s13(6)=2.d0*(r1(re,1)*r3(re, 6)+r1(im,1)*r3(im, 6)+& 788 & r1(re,2)*r3(re,12)+r1(im,2)*r3(im,12)+& 789 & r1(re,3)*r3(re,18)+r1(im,3)*r3(im,18)) 790 791 !Finally, write out the six terms as final answer 792 !rank2(a,b)=(15/2)*s33(a,b)-3*s13(a,b)-(3/2)*2*Re[r1(a)*r1(b)] 793 794 rank2(1)=7.5d0*s33(1)-3.d0*s13(1)& 795 & -3.d0*(r1(re,1)*r1(re,1)+r1(im,1)*r1(im,1)) 796 rank2(2)=7.5d0*s33(2)-3.d0*s13(2)& 797 & -3.d0*(r1(re,2)*r1(re,2)+r1(im,2)*r1(im,2)) 798 rank2(3)=7.5d0*s33(3)-3.d0*s13(3)& 799 & -3.d0*(r1(re,3)*r1(re,3)+r1(im,3)*r1(im,3)) 800 rank2(4)=7.5d0*s33(4)-3.d0*s13(4)& 801 & -3.d0*(r1(re,3)*r1(re,2)+r1(im,3)*r1(im,2)) 802 rank2(5)=7.5d0*s33(5)-3.d0*s13(5)& 803 & -3.d0*(r1(re,3)*r1(re,1)+r1(im,3)*r1(im,1)) 804 rank2(6)=7.5d0*s33(6)-3.d0*s13(6)& 805 & -3.d0*(r1(re,2)*r1(re,1)+r1(im,2)*r1(im,1)) 806 807 end subroutine cont3
m_contract/cont33cso [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont33cso
FUNCTION
Contract symmetric rank 3 tensor gxa1 with symmetric rank 3 tensor gxa2 using metric tensor gmet to produce rank 2 complex tensor.
INPUTS
gxa1(2,10)=rank 3 complex symmetric tensor gxa2(2,10)=rank 3 complex symmetric tensor gmet(3,3)=usual metric tensor (symmetric, real)
OUTPUT
rank2c(2,6)=rank 2 complex tensor (pseudo-symmetric storage)
NOTES
This contraction is used for spin-orbit correction in non-local contribution to stresses. Symmetric gxa1, gxa2 are stored as 111 221 331 321 311 211 222 332 322 333; gmet(3,3) is symmetric but stored fully (9 elements); Output rank2c is not symmetric but since $rank2c_{gxa1,gxa2}(a,b)=conjg(rank2c_{gxa2,gxa1}(b,a))$ it is stored as 11 22 33 32 31 21. rank2c(1,1), rank2c(2,2), rank3c(3,3) are not needed; They are not calculated. rank2c(a,b)=7.5 conjg(gxa1(a,i,j))*r_3(i,j,b) - 1.5 r_{11}(a)*r_{12}(b) where: r_3(i,j,b) & = & gxa2(b,l,m) gmet(i,l) gmet(j,m) \nonumber r_{11}(a) & = & conjg(gxa1(a,l,m)) gmet(l,m) \nonumber r_{12}(b) & = & gxa2(b,l,m) gmet(l,m) \end{eqnarray} }}
SOURCE
852 subroutine cont33cso(gxa1,gxa2,gmet,rank2c) 853 854 !Arguments ------------------------------------ 855 !arrays 856 real(dp),intent(in) :: gmet(3,3),gxa1(2,10),gxa2(2,10) 857 real(dp),intent(out) :: rank2c(2,6) 858 859 !Local variables------------------------------- 860 !scalars 861 integer,parameter :: im=2,re=1 862 !arrays 863 real(dp) :: r11(2,3),r12(2,3),r3(2,6,3),r3a(2,6,3) 864 865 ! ************************************************************************* 866 867 !Initialize output tensor 868 rank2c(:,:)=0.d0 869 870 !Compute r3(i,j,b)=gxa2(b,l,m)*gmet(i,l)*gmet(j,m) 871 !stored as 11 22 33 32 31 21 for (i,j) 872 !First compute r3a(b,l,j)=gxa2(b,l,m)*gmet(j,m) 873 !stored as r3a(re/im,(b,l),j) 874 r3a(:,1,1)=gxa2(:,1 )*gmet(1,1)+gxa2(:,6 )*gmet(1,2)+gxa2(:,5 )*gmet(1,3) 875 r3a(:,2,1)=gxa2(:,2 )*gmet(1,1)+gxa2(:,7 )*gmet(1,2)+gxa2(:,9 )*gmet(1,3) 876 r3a(:,3,1)=gxa2(:,3 )*gmet(1,1)+gxa2(:,8 )*gmet(1,2)+gxa2(:,10)*gmet(1,3) 877 r3a(:,4,1)=gxa2(:,4 )*gmet(1,1)+gxa2(:,9 )*gmet(1,2)+gxa2(:,8 )*gmet(1,3) 878 r3a(:,5,1)=gxa2(:,5 )*gmet(1,1)+gxa2(:,4 )*gmet(1,2)+gxa2(:,3 )*gmet(1,3) 879 r3a(:,6,1)=gxa2(:,6 )*gmet(1,1)+gxa2(:,2 )*gmet(1,2)+gxa2(:,4 )*gmet(1,3) 880 r3a(:,1,2)=gxa2(:,1 )*gmet(2,1)+gxa2(:,6 )*gmet(2,2)+gxa2(:,5 )*gmet(2,3) 881 r3a(:,2,2)=gxa2(:,2 )*gmet(2,1)+gxa2(:,7 )*gmet(2,2)+gxa2(:,9 )*gmet(2,3) 882 r3a(:,3,2)=gxa2(:,3 )*gmet(2,1)+gxa2(:,8 )*gmet(2,2)+gxa2(:,10)*gmet(2,3) 883 r3a(:,4,2)=gxa2(:,4 )*gmet(2,1)+gxa2(:,9 )*gmet(2,2)+gxa2(:,8 )*gmet(2,3) 884 r3a(:,5,2)=gxa2(:,5 )*gmet(2,1)+gxa2(:,4 )*gmet(2,2)+gxa2(:,3 )*gmet(2,3) 885 r3a(:,6,2)=gxa2(:,6 )*gmet(2,1)+gxa2(:,2 )*gmet(2,2)+gxa2(:,4 )*gmet(2,3) 886 r3a(:,1,3)=gxa2(:,1 )*gmet(3,1)+gxa2(:,6 )*gmet(3,2)+gxa2(:,5 )*gmet(3,3) 887 r3a(:,2,3)=gxa2(:,2 )*gmet(3,1)+gxa2(:,7 )*gmet(3,2)+gxa2(:,9 )*gmet(3,3) 888 r3a(:,3,3)=gxa2(:,3 )*gmet(3,1)+gxa2(:,8 )*gmet(3,2)+gxa2(:,10)*gmet(3,3) 889 r3a(:,4,3)=gxa2(:,4 )*gmet(3,1)+gxa2(:,9 )*gmet(3,2)+gxa2(:,8 )*gmet(3,3) 890 r3a(:,5,3)=gxa2(:,5 )*gmet(3,1)+gxa2(:,4 )*gmet(3,2)+gxa2(:,3 )*gmet(3,3) 891 r3a(:,6,3)=gxa2(:,6 )*gmet(3,1)+gxa2(:,2 )*gmet(3,2)+gxa2(:,4 )*gmet(3,3) 892 893 !Then compute r3(i,j,b)=r3a(b,l,j)*gmet(i,l) 894 !stored as r3(re/im,(i,j),b) 895 r3(:,1,1)=r3a(:,1,1)*gmet(1,1)+r3a(:,6,1)*gmet(1,2)+r3a(:,5,1)*gmet(1,3) 896 r3(:,1,2)=r3a(:,6,1)*gmet(1,1)+r3a(:,2,1)*gmet(1,2)+r3a(:,4,1)*gmet(1,3) 897 r3(:,1,3)=r3a(:,5,1)*gmet(1,1)+r3a(:,4,1)*gmet(1,2)+r3a(:,3,1)*gmet(1,3) 898 r3(:,2,1)=r3a(:,1,2)*gmet(2,1)+r3a(:,6,2)*gmet(2,2)+r3a(:,5,2)*gmet(2,3) 899 r3(:,2,2)=r3a(:,6,2)*gmet(2,1)+r3a(:,2,2)*gmet(2,2)+r3a(:,4,2)*gmet(2,3) 900 r3(:,2,3)=r3a(:,5,2)*gmet(2,1)+r3a(:,4,2)*gmet(2,2)+r3a(:,3,2)*gmet(2,3) 901 r3(:,3,1)=r3a(:,1,3)*gmet(3,1)+r3a(:,6,3)*gmet(3,2)+r3a(:,5,3)*gmet(3,3) 902 r3(:,3,2)=r3a(:,6,3)*gmet(3,1)+r3a(:,2,3)*gmet(3,2)+r3a(:,4,3)*gmet(3,3) 903 r3(:,3,3)=r3a(:,5,3)*gmet(3,1)+r3a(:,4,3)*gmet(3,2)+r3a(:,3,3)*gmet(3,3) 904 r3(:,4,1)=r3a(:,1,2)*gmet(3,1)+r3a(:,6,2)*gmet(3,2)+r3a(:,5,2)*gmet(3,3) 905 r3(:,4,2)=r3a(:,6,2)*gmet(3,1)+r3a(:,2,2)*gmet(3,2)+r3a(:,4,2)*gmet(3,3) 906 r3(:,4,3)=r3a(:,5,2)*gmet(3,1)+r3a(:,4,2)*gmet(3,2)+r3a(:,3,2)*gmet(3,3) 907 r3(:,5,1)=r3a(:,1,1)*gmet(3,1)+r3a(:,6,1)*gmet(3,2)+r3a(:,5,1)*gmet(3,3) 908 r3(:,5,2)=r3a(:,6,1)*gmet(3,1)+r3a(:,2,1)*gmet(3,2)+r3a(:,4,1)*gmet(3,3) 909 r3(:,5,3)=r3a(:,5,1)*gmet(3,1)+r3a(:,4,1)*gmet(3,2)+r3a(:,3,1)*gmet(3,3) 910 r3(:,6,1)=r3a(:,1,1)*gmet(2,1)+r3a(:,6,1)*gmet(2,2)+r3a(:,5,1)*gmet(2,3) 911 r3(:,6,2)=r3a(:,6,1)*gmet(2,1)+r3a(:,2,1)*gmet(2,2)+r3a(:,4,1)*gmet(2,3) 912 r3(:,6,3)=r3a(:,5,1)*gmet(2,1)+r3a(:,4,1)*gmet(2,2)+r3a(:,3,1)*gmet(2,3) 913 914 !Compute r11(a)=conjg(gxa1(a,l,m))*gmet(l,m) 915 r11(:,1)=gxa1(:,1)*gmet(1,1)+gxa1(:,2)*gmet(2,2)+gxa1(:,3 )*gmet(3,3)& 916 & +2.d0*(gxa1(:,4)*gmet(3,2)+gxa1(:,5)*gmet(3,1)+gxa1(:,6 )*gmet(2,1)) 917 r11(:,2)=gxa1(:,6)*gmet(1,1)+gxa1(:,7)*gmet(2,2)+gxa1(:,8 )*gmet(3,3)& 918 & +2.d0*(gxa1(:,9)*gmet(3,2)+gxa1(:,4)*gmet(3,1)+gxa1(:,2 )*gmet(2,1)) 919 r11(:,3)=gxa1(:,5)*gmet(1,1)+gxa1(:,9)*gmet(2,2)+gxa1(:,10)*gmet(3,3)& 920 & +2.d0*(gxa1(:,8)*gmet(3,2)+gxa1(:,3)*gmet(3,1)+gxa1(:,4 )*gmet(2,1)) 921 r11(im,1)=-r11(im,1);r11(im,2)=-r11(im,2);r11(im,3)=-r11(im,3) 922 923 !Compute r12(b)=gxa2(b,l,m)*gmet(l,m) 924 r12(:,1)=gxa2(:,1)*gmet(1,1)+gxa2(:,2)*gmet(2,2)+gxa2(:,3 )*gmet(3,3)& 925 & +2.d0*(gxa2(:,4)*gmet(3,2)+gxa2(:,5)*gmet(3,1)+gxa2(:,6 )*gmet(2,1)) 926 r12(:,2)=gxa2(:,6)*gmet(1,1)+gxa2(:,7)*gmet(2,2)+gxa2(:,8 )*gmet(3,3)& 927 & +2.d0*(gxa2(:,9)*gmet(3,2)+gxa2(:,4)*gmet(3,1)+gxa2(:,2 )*gmet(2,1)) 928 r12(:,3)=gxa2(:,5)*gmet(1,1)+gxa2(:,9)*gmet(2,2)+gxa2(:,10)*gmet(3,3)& 929 & +2.d0*(gxa2(:,8)*gmet(3,2)+gxa2(:,3)*gmet(3,1)+gxa2(:,4 )*gmet(2,1)) 930 931 !Finally compute rank2c(a,b)=7.5*conjg(gxa1(a,i,j))*r3(i,j,b) - 1.5*r11(a)*r12(b) 932 !rank2c(re,1)=7.5d0*(gxa1(re,1 )*r3(re,1,1)+gxa1(im,1 )*r3(im,1,1)& 933 !& +gxa1(re,2 )*r3(re,2,1)+gxa1(im,2 )*r3(im,2,1)& 934 !& +gxa1(re,3 )*r3(re,3,1)+gxa1(im,3 )*r3(im,3,1))& 935 !& +15.d0*(gxa1(re,4 )*r3(re,4,1)+gxa1(im,4 )*r3(im,4,1)& 936 !& +gxa1(re,5 )*r3(re,5,1)+gxa1(im,5 )*r3(im,5,1)& 937 !& +gxa1(re,6 )*r3(re,6,1)+gxa1(im,6 )*r3(im,6,1))& 938 !& -1.5d0*(r11(re,1)*r12(re,1)-r11(im,1)*r12(im,1)) 939 !rank2c(re,2)=7.5d0*(gxa1(re,6 )*r3(re,1,2)+gxa1(im,6 )*r3(im,1,2)& 940 !& +gxa1(re,7 )*r3(re,2,2)+gxa1(im,7 )*r3(im,2,2)& 941 !& +gxa1(re,8 )*r3(re,3,2)+gxa1(im,8 )*r3(im,3,2))& 942 !& +15.d0*(gxa1(re,9 )*r3(re,4,2)+gxa1(im,9 )*r3(im,4,2)& 943 !& +gxa1(re,4 )*r3(re,5,2)+gxa1(im,4 )*r3(im,5,2)& 944 !& +gxa1(re,2 )*r3(re,6,2)+gxa1(im,2 )*r3(im,6,2))& 945 !& -1.5d0*(r11(re,2)*r12(re,2)-r11(im,2)*r12(im,2)) 946 !rank2c(re,3)=7.5d0*(gxa1(re,5 )*r3(re,1,3)+gxa1(im,5 )*r3(im,1,3)& 947 !& +gxa1(re,9 )*r3(re,2,3)+gxa1(im,9 )*r3(im,2,3)& 948 !& +gxa1(re,10)*r3(re,3,3)+gxa1(im,10)*r3(im,3,3))& 949 !& +15.d0*(gxa1(re,8 )*r3(re,4,3)+gxa1(im,8 )*r3(im,4,3)& 950 !& +gxa1(re,3 )*r3(re,5,3)+gxa1(im,3 )*r3(im,5,3)& 951 !& +gxa1(re,4 )*r3(re,6,3)+gxa1(im,4 )*r3(im,6,3))& 952 !& -1.5d0*(r11(re,3)*r12(re,3)-r11(im,3)*r12(im,3)) 953 rank2c(re,4)=7.5d0*(gxa1(re,5 )*r3(re,1,2)+gxa1(im,5 )*r3(im,1,2)& 954 & +gxa1(re,9 )*r3(re,2,2)+gxa1(im,9 )*r3(im,2,2)& 955 & +gxa1(re,10)*r3(re,3,2)+gxa1(im,10)*r3(im,3,2))& 956 & +15.d0*(gxa1(re,8 )*r3(re,4,2)+gxa1(im,8 )*r3(im,4,2)& 957 & +gxa1(re,3 )*r3(re,5,2)+gxa1(im,3 )*r3(im,5,2)& 958 & +gxa1(re,4 )*r3(re,6,2)+gxa1(im,4 )*r3(im,6,2))& 959 & -1.5d0*(r11(re,3)*r12(re,2)-r11(im,3)*r12(im,2)) 960 rank2c(re,5)=7.5d0*(gxa1(re,5 )*r3(re,1,1)+gxa1(im,5 )*r3(im,1,1)& 961 & +gxa1(re,9 )*r3(re,2,1)+gxa1(im,9 )*r3(im,2,1)& 962 & +gxa1(re,10)*r3(re,3,1)+gxa1(im,10)*r3(im,3,1))& 963 & +15.d0*(gxa1(re,8 )*r3(re,4,1)+gxa1(im,8 )*r3(im,4,1)& 964 & +gxa1(re,3 )*r3(re,5,1)+gxa1(im,3 )*r3(im,5,1)& 965 & +gxa1(re,4 )*r3(re,6,1)+gxa1(im,4 )*r3(im,6,1))& 966 & -1.5d0*(r11(re,3)*r12(re,1)-r11(im,3)*r12(im,1)) 967 rank2c(re,6)=7.5d0*(gxa1(re,6 )*r3(re,1,1)+gxa1(im,6 )*r3(im,1,1)& 968 & +gxa1(re,7 )*r3(re,2,1)+gxa1(im,7 )*r3(im,2,1)& 969 & +gxa1(re,8 )*r3(re,3,1)+gxa1(im,8 )*r3(im,3,1))& 970 & +15.d0*(gxa1(re,9 )*r3(re,4,1)+gxa1(im,9 )*r3(im,4,1)& 971 & +gxa1(re,4 )*r3(re,5,1)+gxa1(im,4 )*r3(im,5,1)& 972 & +gxa1(re,2 )*r3(re,6,1)+gxa1(im,2 )*r3(im,6,1))& 973 & -1.5d0*(r11(re,2)*r12(re,1)-r11(im,2)*r12(im,1)) 974 !rank2c(im,1)=7.5d0*(gxa1(re,1 )*r3(im,1,1)-gxa1(im,1 )*r3(re,1,1)& 975 !& +gxa1(re,2 )*r3(im,2,1)-gxa1(im,2 )*r3(re,2,1)& 976 !& +gxa1(re,3 )*r3(im,3,1)-gxa1(im,3 )*r3(re,3,1))& 977 !& +15.d0*(gxa1(re,4 )*r3(im,4,1)-gxa1(im,4 )*r3(re,4,1)& 978 !& +gxa1(re,5 )*r3(im,5,1)-gxa1(im,5 )*r3(re,5,1)& 979 !& +gxa1(re,6 )*r3(im,6,1)-gxa1(im,6 )*r3(re,6,1))& 980 !& -1.5d0*(r11(re,1)*r12(im,1)+r11(im,1)*r12(re,1)) 981 !rank2c(im,2)=7.5d0*(gxa1(re,6 )*r3(im,1,2)-gxa1(im,6 )*r3(re,1,2)& 982 !& +gxa1(re,7 )*r3(im,2,2)-gxa1(im,7 )*r3(re,2,2)& 983 !& +gxa1(re,8 )*r3(im,3,2)-gxa1(im,8 )*r3(re,3,2))& 984 !& +15.d0*(gxa1(re,9 )*r3(im,4,2)-gxa1(im,9 )*r3(re,4,2)& 985 !& +gxa1(re,4 )*r3(im,5,2)-gxa1(im,4 )*r3(re,5,2)& 986 !& +gxa1(re,2 )*r3(im,6,2)-gxa1(im,2 )*r3(re,6,2))& 987 !& -1.5d0*(r11(re,2)*r12(im,2)+r11(im,2)*r12(re,2)) 988 !rank2c(im,3)=7.5d0*(gxa1(re,5 )*r3(im,1,3)-gxa1(im,5 )*r3(re,1,3)& 989 !& +gxa1(re,9 )*r3(im,2,3)-gxa1(im,9 )*r3(re,2,3)& 990 !& +gxa1(re,10)*r3(im,3,3)-gxa1(im,10)*r3(re,3,3))& 991 !& +15.d0*(gxa1(re,8 )*r3(im,4,3)-gxa1(im,8 )*r3(re,4,3)& 992 !& +gxa1(re,3 )*r3(im,5,3)-gxa1(im,3 )*r3(re,5,3)& 993 !& +gxa1(re,4 )*r3(im,6,3)-gxa1(im,4 )*r3(re,6,3))& 994 !& -1.5d0*(r11(re,3)*r12(im,3)+r11(im,3)*r12(re,3)) 995 rank2c(im,4)=7.5d0*(gxa1(re,5 )*r3(im,1,2)-gxa1(im,5 )*r3(re,1,2)& 996 & +gxa1(re,9 )*r3(im,2,2)-gxa1(im,9 )*r3(re,2,2)& 997 & +gxa1(re,10)*r3(im,3,2)-gxa1(im,10)*r3(re,3,2))& 998 & +15.d0*(gxa1(re,8 )*r3(im,4,2)-gxa1(im,8 )*r3(re,4,2)& 999 & +gxa1(re,3 )*r3(im,5,2)-gxa1(im,3 )*r3(re,5,2)& 1000 & +gxa1(re,4 )*r3(im,6,2)-gxa1(im,4 )*r3(re,6,2))& 1001 & -1.5d0*(r11(re,3)*r12(im,2)+r11(im,3)*r12(re,2)) 1002 rank2c(im,5)=7.5d0*(gxa1(re,5 )*r3(im,1,1)-gxa1(im,5 )*r3(re,1,1)& 1003 & +gxa1(re,9 )*r3(im,2,1)-gxa1(im,9 )*r3(re,2,1)& 1004 & +gxa1(re,10)*r3(im,3,1)-gxa1(im,10)*r3(re,3,1))& 1005 & +15.d0*(gxa1(re,8 )*r3(im,4,1)-gxa1(im,8 )*r3(re,4,1)& 1006 & +gxa1(re,3 )*r3(im,5,1)-gxa1(im,3 )*r3(re,5,1)& 1007 & +gxa1(re,4 )*r3(im,6,1)-gxa1(im,4 )*r3(re,6,1))& 1008 & -1.5d0*(r11(re,3)*r12(im,1)+r11(im,3)*r12(re,1)) 1009 rank2c(im,6)=7.5d0*(gxa1(re,6 )*r3(im,1,1)-gxa1(im,6 )*r3(re,1,1)& 1010 & +gxa1(re,7 )*r3(im,2,1)-gxa1(im,7 )*r3(re,2,1)& 1011 & +gxa1(re,8 )*r3(im,3,1)-gxa1(im,8 )*r3(re,3,1))& 1012 & +15.d0*(gxa1(re,9 )*r3(im,4,1)-gxa1(im,9 )*r3(re,4,1)& 1013 & +gxa1(re,4 )*r3(im,5,1)-gxa1(im,4 )*r3(re,5,1)& 1014 & +gxa1(re,2 )*r3(im,6,1)-gxa1(im,2 )*r3(re,6,1))& 1015 & -1.5d0*(r11(re,2)*r12(im,1)+r11(im,2)*r12(re,1)) 1016 1017 end subroutine cont33cso
m_contract/cont33so [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont33so
FUNCTION
Contract symmetric rank 3 tensor gxa1 with symmetric rank 3 tensor gxa2 using metric tensor gmet and antisymmetric tensor amet to produce rank 2 real tensor.
INPUTS
gxa1(2,10)=rank 3 complex symmetric tensor gxa2(2,10)=rank 3 complex symmetric tensor gmet(3,3)=usual metric tensor (symmetric, real) amet(2,3,3)=antisymmetric complex tensor used for spin-orbit
OUTPUT
rank2(6)=rank 2 real tensor (pseudo-symmetric storage)
NOTES
This contraction is used for spin-orbit correction in non-local contribution to stresses. Symmetric gxa1, gxa2 are stored as 111 221 331 321 311 211 222 332 322 333; gmet(3,3) is symmetric but stored fully (9 elements); amet(3,3) is antisymmetric but stored fully (9 elements); Output rank2 is not symmetric but since $rank2_{gxa1,gxa2}(a,b)=conjg(rank2_{gxa2,gxa1}(b,a))$ it is stored as 11 22 33 32 31 21. Want 2*Re[contraction]. rank2(a,b)=2 Re[15 r_{3A}(a,i,j) r_{3G}(b,j,i)-3 r_{3A}(a,b,i) r_1(i)] where: r_1(i) & = & gxa2(i,j,k) gmet(j,k) \nonumber r_{3A}(i,j,k) & = & conjg(gxa1(p,i,j)) amet(p,k) \nonumber r_{3G}(i,j,k) & = & gxa2(p,i,j) gmet(p,k) \end{eqnarray} }}
SOURCE
1064 subroutine cont33so(gxa1,gxa2,gmet,amet,rank2) 1065 1066 !Arguments ------------------------------------ 1067 !arrays 1068 real(dp),intent(in) :: amet(2,3,3),gmet(3,3),gxa1(2,10),gxa2(2,10) 1069 real(dp),intent(out) :: rank2(6) 1070 1071 !Local variables------------------------------- 1072 !scalars 1073 integer,parameter :: im=2,re=1 1074 integer :: ii 1075 !arrays 1076 real(dp) :: r1(2,3),r3A(2,18),r3G(2,18),s31(6),s33(6) 1077 1078 ! ************************************************************************* 1079 1080 !Compute r1(i)=gxa2(i,j,k)*gmet(j,k) 1081 1082 do ii=1,2 1083 r1(ii,1)=gmet(1,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,2)+& 1084 & gmet(3,3)*gxa2(ii,3)+2.d0*(& 1085 & gmet(3,2)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,5)+& 1086 & gmet(2,1)*gxa2(ii,6)) 1087 r1(ii,2)=gmet(1,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,7)+& 1088 & gmet(3,3)*gxa2(ii,8)+2.d0*(& 1089 & gmet(3,2)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,4)+& 1090 & gmet(2,1)*gxa2(ii,2)) 1091 r1(ii,3)=gmet(1,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,9)+& 1092 & gmet(3,3)*gxa2(ii,10)+2.d0*(& 1093 & gmet(3,2)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,3)+& 1094 & gmet(2,1)*gxa2(ii,4)) 1095 end do 1096 1097 !Compute r3G(i,j,k)=gxa2(p,i,j)*gmet(p,k) 1098 !Write out components for 18 distinct terms, Re and Im 1099 !(r3G is symmetric in first two indices, not in all permutations) 1100 !Store as 111 221 331 321 311 211 1101 !112 222 332 322 312 212 1102 !113 223 333 323 313 213 1103 do ii=1,2 1104 r3G(ii, 1)=gmet(1,1)*gxa2(ii,1)+gmet(2,1)*gxa2(ii,6)+gmet(3,1)*gxa2(ii,5) 1105 r3G(ii, 2)=gmet(1,1)*gxa2(ii,2)+gmet(2,1)*gxa2(ii,7)+gmet(3,1)*gxa2(ii,9) 1106 r3G(ii, 3)=gmet(1,1)*gxa2(ii,3)+gmet(2,1)*gxa2(ii,8)+gmet(3,1)*gxa2(ii,10) 1107 r3G(ii, 4)=gmet(1,1)*gxa2(ii,4)+gmet(2,1)*gxa2(ii,9)+gmet(3,1)*gxa2(ii,8) 1108 r3G(ii, 5)=gmet(1,1)*gxa2(ii,5)+gmet(2,1)*gxa2(ii,4)+gmet(3,1)*gxa2(ii,3) 1109 r3G(ii, 6)=gmet(1,1)*gxa2(ii,6)+gmet(2,1)*gxa2(ii,2)+gmet(3,1)*gxa2(ii,4) 1110 r3G(ii, 7)=gmet(2,1)*gxa2(ii,1)+gmet(2,2)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,5) 1111 r3G(ii, 8)=gmet(2,1)*gxa2(ii,2)+gmet(2,2)*gxa2(ii,7)+gmet(3,2)*gxa2(ii,9) 1112 r3G(ii, 9)=gmet(2,1)*gxa2(ii,3)+gmet(2,2)*gxa2(ii,8)+gmet(3,2)*gxa2(ii,10) 1113 r3G(ii,10)=gmet(2,1)*gxa2(ii,4)+gmet(2,2)*gxa2(ii,9)+gmet(3,2)*gxa2(ii,8) 1114 r3G(ii,11)=gmet(2,1)*gxa2(ii,5)+gmet(2,2)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,3) 1115 r3G(ii,12)=gmet(2,1)*gxa2(ii,6)+gmet(2,2)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,4) 1116 r3G(ii,13)=gmet(3,1)*gxa2(ii,1)+gmet(3,2)*gxa2(ii,6)+gmet(3,3)*gxa2(ii,5) 1117 r3G(ii,14)=gmet(3,1)*gxa2(ii,2)+gmet(3,2)*gxa2(ii,7)+gmet(3,3)*gxa2(ii,9) 1118 r3G(ii,15)=gmet(3,1)*gxa2(ii,3)+gmet(3,2)*gxa2(ii,8)+gmet(3,3)*gxa2(ii,10) 1119 r3G(ii,16)=gmet(3,1)*gxa2(ii,4)+gmet(3,2)*gxa2(ii,9)+gmet(3,3)*gxa2(ii,8) 1120 r3G(ii,17)=gmet(3,1)*gxa2(ii,5)+gmet(3,2)*gxa2(ii,4)+gmet(3,3)*gxa2(ii,3) 1121 r3G(ii,18)=gmet(3,1)*gxa2(ii,6)+gmet(3,2)*gxa2(ii,2)+gmet(3,3)*gxa2(ii,4) 1122 end do 1123 1124 !Compute r3A(i,j,k)=conjg(gxa1(p,i,j))*amet(p,k) 1125 !Write out components for 18 distinct terms, Re and Im 1126 !(r3A is symmetric in first two indices, not in all permutations) 1127 !Store as 111 221 331 321 311 211 1128 !112 222 332 322 312 212 1129 !113 223 333 323 313 213 1130 !Note that, since amet is antisymmetric, amet(i,i)=0 1131 1132 r3A(re, 1)=amet(re,2,1)*gxa1(re,6 )+amet(im,2,1)*gxa1(im,6 )+& 1133 & amet(re,3,1)*gxa1(re,5 )+amet(im,3,1)*gxa1(im,5 ) 1134 r3A(re, 2)=amet(re,2,1)*gxa1(re,7 )+amet(im,2,1)*gxa1(im,7 )+& 1135 & amet(re,3,1)*gxa1(re,9 )+amet(im,3,1)*gxa1(im,9 ) 1136 r3A(re, 3)=amet(re,2,1)*gxa1(re,8 )+amet(im,2,1)*gxa1(im,8 )+& 1137 & amet(re,3,1)*gxa1(re,10)+amet(im,3,1)*gxa1(im,10) 1138 r3A(re, 4)=amet(re,2,1)*gxa1(re,9 )+amet(im,2,1)*gxa1(im,9 )+& 1139 & amet(re,3,1)*gxa1(re,8 )+amet(im,3,1)*gxa1(im,8 ) 1140 r3A(re, 5)=amet(re,2,1)*gxa1(re,4 )+amet(im,2,1)*gxa1(im,4 )+& 1141 & amet(re,3,1)*gxa1(re,3 )+amet(im,3,1)*gxa1(im,3 ) 1142 r3A(re, 6)=amet(re,2,1)*gxa1(re,2 )+amet(im,2,1)*gxa1(im,2 )+& 1143 & amet(re,3,1)*gxa1(re,4 )+amet(im,3,1)*gxa1(im,4 ) 1144 r3A(re, 7)=amet(re,1,2)*gxa1(re,1 )+amet(im,1,2)*gxa1(im,1 )+& 1145 & amet(re,3,2)*gxa1(re,5 )+amet(im,3,2)*gxa1(im,5 ) 1146 r3A(re, 8)=amet(re,1,2)*gxa1(re,2 )+amet(im,1,2)*gxa1(im,2 )+& 1147 & amet(re,3,2)*gxa1(re,9 )+amet(im,3,2)*gxa1(im,9 ) 1148 r3A(re, 9)=amet(re,1,2)*gxa1(re,3 )+amet(im,1,2)*gxa1(im,3 )+& 1149 & amet(re,3,2)*gxa1(re,10)+amet(im,3,2)*gxa1(im,10) 1150 r3A(re,10)=amet(re,1,2)*gxa1(re,4 )+amet(im,1,2)*gxa1(im,4 )+& 1151 & amet(re,3,2)*gxa1(re,8 )+amet(im,3,2)*gxa1(im,8 ) 1152 r3A(re,11)=amet(re,1,2)*gxa1(re,5 )+amet(im,1,2)*gxa1(im,5 )+& 1153 & amet(re,3,2)*gxa1(re,3 )+amet(im,3,2)*gxa1(im,3 ) 1154 r3A(re,12)=amet(re,1,2)*gxa1(re,6 )+amet(im,1,2)*gxa1(im,6 )+& 1155 & amet(re,3,2)*gxa1(re,4 )+amet(im,3,2)*gxa1(im,4 ) 1156 r3A(re,13)=amet(re,1,3)*gxa1(re,1 )+amet(im,1,3)*gxa1(im,1 )+& 1157 & amet(re,2,3)*gxa1(re,6 )+amet(im,2,3)*gxa1(im,6 ) 1158 r3A(re,14)=amet(re,1,3)*gxa1(re,2 )+amet(im,1,3)*gxa1(im,2 )+& 1159 & amet(re,2,3)*gxa1(re,7 )+amet(im,2,3)*gxa1(im,7 ) 1160 r3A(re,15)=amet(re,1,3)*gxa1(re,3 )+amet(im,1,3)*gxa1(im,3 )+& 1161 & amet(re,2,3)*gxa1(re,8 )+amet(im,2,3)*gxa1(im,8 ) 1162 r3A(re,16)=amet(re,1,3)*gxa1(re,4 )+amet(im,1,3)*gxa1(im,4 )+& 1163 & amet(re,2,3)*gxa1(re,9 )+amet(im,2,3)*gxa1(im,9 ) 1164 r3A(re,17)=amet(re,1,3)*gxa1(re,5 )+amet(im,1,3)*gxa1(im,5 )+& 1165 & amet(re,2,3)*gxa1(re,4 )+amet(im,2,3)*gxa1(im,4 ) 1166 r3A(re,18)=amet(re,1,3)*gxa1(re,6 )+amet(im,1,3)*gxa1(im,6 )+& 1167 & amet(re,2,3)*gxa1(re,2 )+amet(im,2,3)*gxa1(im,2 ) 1168 1169 r3A(im, 1)=amet(im,2,1)*gxa1(re,6 )-amet(re,2,1)*gxa1(im,6 )+& 1170 & amet(im,3,1)*gxa1(re,5 )-amet(re,3,1)*gxa1(im,5 ) 1171 r3A(im, 2)=amet(im,2,1)*gxa1(re,7 )-amet(re,2,1)*gxa1(im,7 )+& 1172 & amet(im,3,1)*gxa1(re,9 )-amet(re,3,1)*gxa1(im,9 ) 1173 r3A(im, 3)=amet(im,2,1)*gxa1(re,8 )-amet(re,2,1)*gxa1(im,8 )+& 1174 & amet(im,3,1)*gxa1(re,10)-amet(re,3,1)*gxa1(im,10) 1175 r3A(im, 4)=amet(im,2,1)*gxa1(re,9 )-amet(re,2,1)*gxa1(im,9 )+& 1176 & amet(im,3,1)*gxa1(re,8 )-amet(re,3,1)*gxa1(im,8 ) 1177 r3A(im, 5)=amet(im,2,1)*gxa1(re,4 )-amet(re,2,1)*gxa1(im,4 )+& 1178 & amet(im,3,1)*gxa1(re,3 )-amet(re,3,1)*gxa1(im,3 ) 1179 r3A(im, 6)=amet(im,2,1)*gxa1(re,2 )-amet(re,2,1)*gxa1(im,2 )+& 1180 & amet(im,3,1)*gxa1(re,4 )-amet(re,3,1)*gxa1(im,4 ) 1181 r3A(im, 7)=amet(im,1,2)*gxa1(re,1 )-amet(re,1,2)*gxa1(im,1 )+& 1182 & amet(im,3,2)*gxa1(re,5 )-amet(re,3,2)*gxa1(im,5 ) 1183 r3A(im, 8)=amet(im,1,2)*gxa1(re,2 )-amet(re,1,2)*gxa1(im,2 )+& 1184 & amet(im,3,2)*gxa1(re,9 )-amet(re,3,2)*gxa1(im,9 ) 1185 r3A(im, 9)=amet(im,1,2)*gxa1(re,3 )-amet(re,1,2)*gxa1(im,3 )+& 1186 & amet(im,3,2)*gxa1(re,10)-amet(re,3,2)*gxa1(im,10) 1187 r3A(im,10)=amet(im,1,2)*gxa1(re,4 )-amet(re,1,2)*gxa1(im,4 )+& 1188 & amet(im,3,2)*gxa1(re,8 )-amet(re,3,2)*gxa1(im,8 ) 1189 r3A(im,11)=amet(im,1,2)*gxa1(re,5 )-amet(re,1,2)*gxa1(im,5 )+& 1190 & amet(im,3,2)*gxa1(re,3 )-amet(re,3,2)*gxa1(im,3 ) 1191 r3A(im,12)=amet(im,1,2)*gxa1(re,6 )-amet(re,1,2)*gxa1(im,6 )+& 1192 & amet(im,3,2)*gxa1(re,4 )-amet(re,3,2)*gxa1(im,4 ) 1193 r3A(im,13)=amet(im,1,3)*gxa1(re,1 )-amet(re,1,3)*gxa1(im,1 )+& 1194 & amet(im,2,3)*gxa1(re,6 )-amet(re,2,3)*gxa1(im,6 ) 1195 r3A(im,14)=amet(im,1,3)*gxa1(re,2 )-amet(re,1,3)*gxa1(im,2 )+& 1196 & amet(im,2,3)*gxa1(re,7 )-amet(re,2,3)*gxa1(im,7 ) 1197 r3A(im,15)=amet(im,1,3)*gxa1(re,3 )-amet(re,1,3)*gxa1(im,3 )+& 1198 & amet(im,2,3)*gxa1(re,8 )-amet(re,2,3)*gxa1(im,8 ) 1199 r3A(im,16)=amet(im,1,3)*gxa1(re,4 )-amet(re,1,3)*gxa1(im,4 )+& 1200 & amet(im,2,3)*gxa1(re,9 )-amet(re,2,3)*gxa1(im,9 ) 1201 r3A(im,17)=amet(im,1,3)*gxa1(re,5 )-amet(re,1,3)*gxa1(im,5 )+& 1202 & amet(im,2,3)*gxa1(re,4 )-amet(re,2,3)*gxa1(im,4 ) 1203 r3A(im,18)=amet(im,1,3)*gxa1(re,6 )-amet(re,1,3)*gxa1(im,6 )+& 1204 & amet(im,2,3)*gxa1(re,2 )-amet(re,2,3)*gxa1(im,2 ) 1205 1206 !Compute s33(a,b)=2*Re[r3A(a,i,j)*r3G(b,j,i)] 1207 1208 s33(1)=2.d0*(r3A(re, 1)*r3G(re, 1)-r3A(im, 1)*r3G(im, 1)+& 1209 & r3A(re,12)*r3G(re,12)-r3A(im,12)*r3G(im,12)+& 1210 & r3A(re,17)*r3G(re,17)-r3A(im,17)*r3G(im,17)+& 1211 & r3A(re,11)*r3G(re,18)-r3A(im,11)*r3G(im,18)+& 1212 & r3A(re,18)*r3G(re,11)-r3A(im,18)*r3G(im,11)+& 1213 & r3A(re, 5)*r3G(re,13)-r3A(im, 5)*r3G(im,13)+& 1214 & r3A(re,13)*r3G(re, 5)-r3A(im,13)*r3G(im, 5)+& 1215 & r3A(re, 6)*r3G(re, 7)-r3A(im, 6)*r3G(im, 7)+& 1216 & r3A(re, 7)*r3G(re, 6)-r3A(im, 7)*r3G(im, 6)) 1217 s33(2)=2.d0*(r3A(re, 6)*r3G(re, 6)-r3A(im, 6)*r3G(im, 6)+& 1218 & r3A(re, 8)*r3G(re, 8)-r3A(im, 8)*r3G(im, 8)+& 1219 & r3A(re,16)*r3G(re,16)-r3A(im,16)*r3G(im,16)+& 1220 & r3A(re,10)*r3G(re,14)-r3A(im,10)*r3G(im,14)+& 1221 & r3A(re,14)*r3G(re,10)-r3A(im,14)*r3G(im,10)+& 1222 & r3A(re, 4)*r3G(re,18)-r3A(im, 4)*r3G(im,18)+& 1223 & r3A(re,18)*r3G(re, 4)-r3A(im,18)*r3G(im, 4)+& 1224 & r3A(re, 2)*r3G(re,12)-r3A(im, 2)*r3G(im,12)+& 1225 & r3A(re,12)*r3G(re, 2)-r3A(im,12)*r3G(im, 2)) 1226 s33(3)=2.d0*(r3A(re, 5)*r3G(re, 5)-r3A(im, 5)*r3G(im, 5)+& 1227 & r3A(re,10)*r3G(re,10)-r3A(im,10)*r3G(im,10)+& 1228 & r3A(re,15)*r3G(re,15)-r3A(im,15)*r3G(im,15)+& 1229 & r3A(re, 9)*r3G(re,16)-r3A(im, 9)*r3G(im,16)+& 1230 & r3A(re,16)*r3G(re, 9)-r3A(im,16)*r3G(im, 9)+& 1231 & r3A(re, 3)*r3G(re,17)-r3A(im, 3)*r3G(im,17)+& 1232 & r3A(re,17)*r3G(re, 3)-r3A(im,17)*r3G(im, 3)+& 1233 & r3A(re, 4)*r3G(re,11)-r3A(im, 4)*r3G(im,11)+& 1234 & r3A(re,11)*r3G(re, 4)-r3A(im,11)*r3G(im, 4)) 1235 s33(4)=2.d0*(r3A(re, 5)*r3G(re, 6)-r3A(im, 5)*r3G(im, 6)+& 1236 & r3A(re,10)*r3G(re, 8)-r3A(im,10)*r3G(im, 8)+& 1237 & r3A(re,15)*r3G(re,16)-r3A(im,15)*r3G(im,16)+& 1238 & r3A(re, 9)*r3G(re,14)-r3A(im, 9)*r3G(im,14)+& 1239 & r3A(re,16)*r3G(re,10)-r3A(im,16)*r3G(im,10)+& 1240 & r3A(re, 3)*r3G(re,18)-r3A(im, 3)*r3G(im,18)+& 1241 & r3A(re,17)*r3G(re, 4)-r3A(im,17)*r3G(im, 4)+& 1242 & r3A(re, 4)*r3G(re,12)-r3A(im, 4)*r3G(im,12)+& 1243 & r3A(re,11)*r3G(re, 2)-r3A(im,11)*r3G(im, 2)) 1244 s33(5)=2.d0*(r3A(re, 5)*r3G(re, 1)-r3A(im, 5)*r3G(im, 1)+& 1245 & r3A(re,10)*r3G(re,12)-r3A(im,10)*r3G(im,12)+& 1246 & r3A(re,15)*r3G(re,17)-r3A(im,15)*r3G(im,17)+& 1247 & r3A(re, 9)*r3G(re,18)-r3A(im, 9)*r3G(im,18)+& 1248 & r3A(re,16)*r3G(re,11)-r3A(im,16)*r3G(im,11)+& 1249 & r3A(re, 3)*r3G(re,13)-r3A(im, 3)*r3G(im,13)+& 1250 & r3A(re,17)*r3G(re, 5)-r3A(im,17)*r3G(im, 5)+& 1251 & r3A(re, 4)*r3G(re, 7)-r3A(im, 4)*r3G(im, 7)+& 1252 & r3A(re,11)*r3G(re, 6)-r3A(im,11)*r3G(im, 6)) 1253 s33(6)=2.d0*(r3A(re, 6)*r3G(re, 1)-r3A(im, 6)*r3G(im, 1)+& 1254 & r3A(re, 8)*r3G(re,12)-r3A(im, 8)*r3G(im,12)+& 1255 & r3A(re,16)*r3G(re,17)-r3A(im,16)*r3G(im,17)+& 1256 & r3A(re,10)*r3G(re,18)-r3A(im,10)*r3G(im,18)+& 1257 & r3A(re,14)*r3G(re,11)-r3A(im,14)*r3G(im,11)+& 1258 & r3A(re, 4)*r3G(re,13)-r3A(im, 4)*r3G(im,13)+& 1259 & r3A(re,18)*r3G(re, 5)-r3A(im,18)*r3G(im, 5)+& 1260 & r3A(re, 2)*r3G(re, 7)-r3A(im, 2)*r3G(im, 7)+& 1261 & r3A(re,12)*r3G(re, 6)-r3A(im,12)*r3G(im, 6)) 1262 1263 !Compute s31(a,b)=2*Re[r3A(a,b,i)*r1(i)] 1264 1265 s31(1)=2.d0*(r1(re,1)*r3A(re, 1)-r1(im,1)*r3A(im, 1)+& 1266 & r1(re,2)*r3A(re, 7)-r1(im,2)*r3A(im, 7)+& 1267 & r1(re,3)*r3A(re,13)-r1(im,3)*r3A(im,13)) 1268 s31(2)=2.d0*(r1(re,1)*r3A(re, 2)-r1(im,1)*r3A(im, 2)+& 1269 & r1(re,2)*r3A(re, 8)-r1(im,2)*r3A(im, 8)+& 1270 & r1(re,3)*r3A(re,14)-r1(im,3)*r3A(im,14)) 1271 s31(3)=2.d0*(r1(re,1)*r3A(re, 3)-r1(im,1)*r3A(im, 3)+& 1272 & r1(re,2)*r3A(re, 9)-r1(im,2)*r3A(im, 9)+& 1273 & r1(re,3)*r3A(re,15)-r1(im,3)*r3A(im,15)) 1274 s31(4)=2.d0*(r1(re,1)*r3A(re, 4)-r1(im,1)*r3A(im, 4)+& 1275 & r1(re,2)*r3A(re,10)-r1(im,2)*r3A(im,10)+& 1276 & r1(re,3)*r3A(re,16)-r1(im,3)*r3A(im,16)) 1277 s31(5)=2.d0*(r1(re,1)*r3A(re, 5)-r1(im,1)*r3A(im, 5)+& 1278 & r1(re,2)*r3A(re,11)-r1(im,2)*r3A(im,11)+& 1279 & r1(re,3)*r3A(re,17)-r1(im,3)*r3A(im,17)) 1280 s31(6)=2.d0*(r1(re,1)*r3A(re, 6)-r1(im,1)*r3A(im, 6)+& 1281 & r1(re,2)*r3A(re,12)-r1(im,2)*r3A(im,12)+& 1282 & r1(re,3)*r3A(re,18)-r1(im,3)*r3A(im,18)) 1283 1284 !Finally, compute rank2(a,b)=-15*s33(a,b)+3*s31(a,b) 1285 1286 rank2(:)=15.d0*s33(:)-3.d0*s31(:) 1287 1288 end subroutine cont33so
m_contract/cont35 [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
cont35
FUNCTION
Contract symmetric rank3 tensor gxa with rank5 symmetric tensor to produce symmetric rank2 tensor.
INPUTS
gxa(2,10)=rank 3 symmetric complex tensor in order rank5(2,21)=rank 5 complex tensor (symmetric storage)
OUTPUT
rank2(6)=rank 2 real tensor (symmetric storage) 11 22 33 32 31 21.
NOTES
Tensors are in "symmetric" storage mode. For rank 3 tensor gxa this is 111 221 331 321 311 211 222 332 322 333; For rank 5 tensor rank5 this is 11111 22111 33111 32111 31111 21111 22211 33211 32211 33311 22221 33221 32221 33321 33331 22222 33222 32222 33322 33332 33333; For rank 2 tensor rank2 this is 11, 22, 33, 32, 31, 21; gxa and rank5 are complex; rank2 is real. Want $2 Re[contraction]$. $rank2(a,b)=2 Re[gxa(i,j,k)^"*" rank5(a,b,i,j,k)]$. Note that the input gxa is typically the result of gxa(i,j,k)=[{5 \over 2} gmet(i,l) gmet(j,m) gmet(k,n) - {3 \over 2} gmet(i,j) gmet(l,m) gmet(k,n)] gxa_old(l,m) where the subroutine "metcon" already includes weights in the definition of gxa for off-diagonal elements.
SOURCE
1327 subroutine cont35(gxa,rank5,rank2) 1328 1329 !Arguments ------------------------------------ 1330 !arrays 1331 real(dp),intent(in) :: gxa(2,10),rank5(2,21) 1332 real(dp),intent(out) :: rank2(6) 1333 1334 !Local variables------------------------------- 1335 !scalars 1336 integer,parameter :: im=2,re=1 1337 1338 ! ************************************************************************* 1339 1340 !Simply write out index summations 1341 1342 !a=1, b=1 in rank2(a,b) --> maps to index 1 1343 rank2(1)=2.0d0*(& 1344 & gxa(re, 1)*rank5(re, 1)+gxa(im, 1)*rank5(im, 1)+& 1345 & gxa(re, 7)*rank5(re, 7)+gxa(im, 7)*rank5(im, 7)+& 1346 & gxa(re,10)*rank5(re,10)+gxa(im,10)*rank5(im,10)+& 1347 & gxa(re, 2)*rank5(re, 2)+gxa(im, 2)*rank5(im, 2)+& 1348 & gxa(re, 3)*rank5(re, 3)+gxa(im, 3)*rank5(im, 3)+& 1349 & gxa(re, 5)*rank5(re, 5)+gxa(im, 5)*rank5(im, 5)+& 1350 & gxa(re, 6)*rank5(re, 6)+gxa(im, 6)*rank5(im, 6)+& 1351 & gxa(re, 8)*rank5(re, 8)+gxa(im, 8)*rank5(im, 8)+& 1352 & gxa(re, 9)*rank5(re, 9)+gxa(im, 9)*rank5(im, 9)+& 1353 & gxa(re, 4)*rank5(re, 4)+gxa(im, 4)*rank5(im, 4)) 1354 1355 1356 !a=2, b=2 in rank2(a,b) --> maps to index 2 1357 rank2(2)=2.0d0*(& 1358 & gxa(re, 1)*rank5(re, 2)+gxa(im, 1)*rank5(im, 2)+& 1359 & gxa(re, 7)*rank5(re,16)+gxa(im, 7)*rank5(im,16)+& 1360 & gxa(re,10)*rank5(re,19)+gxa(im,10)*rank5(im,19)+& 1361 & gxa(re, 2)*rank5(re,11)+gxa(im, 2)*rank5(im,11)+& 1362 & gxa(re, 3)*rank5(re,12)+gxa(im, 3)*rank5(im,12)+& 1363 & gxa(re, 5)*rank5(re, 9)+gxa(im, 5)*rank5(im, 9)+& 1364 & gxa(re, 6)*rank5(re, 7)+gxa(im, 6)*rank5(im, 7)+& 1365 & gxa(re, 8)*rank5(re,17)+gxa(im, 8)*rank5(im,17)+& 1366 & gxa(re, 9)*rank5(re,18)+gxa(im, 9)*rank5(im,18)+& 1367 & gxa(re, 4)*rank5(re,13)+gxa(im, 4)*rank5(im,13)) 1368 1369 !a=3, b=3 in rank2(a,b) --> maps to index 3 1370 rank2(3)=2.0d0*(& 1371 & gxa(re, 1)*rank5(re, 3)+gxa(im, 1)*rank5(im, 3)+& 1372 & gxa(re, 7)*rank5(re,17)+gxa(im, 7)*rank5(im,17)+& 1373 & gxa(re,10)*rank5(re,21)+gxa(im,10)*rank5(im,21)+& 1374 & gxa(re, 2)*rank5(re,12)+gxa(im, 2)*rank5(im,12)+& 1375 & gxa(re, 3)*rank5(re,15)+gxa(im, 3)*rank5(im,15)+& 1376 & gxa(re, 5)*rank5(re,10)+gxa(im, 5)*rank5(im,10)+& 1377 & gxa(re, 6)*rank5(re, 8)+gxa(im, 6)*rank5(im, 8)+& 1378 & gxa(re, 8)*rank5(re,20)+gxa(im, 8)*rank5(im,20)+& 1379 & gxa(re, 9)*rank5(re,19)+gxa(im, 9)*rank5(im,19)+& 1380 & gxa(re, 4)*rank5(re,14)+gxa(im, 4)*rank5(im,14)) 1381 1382 !a=3, b=2 in rank2(a,b) --> maps to index 4 1383 rank2(4)=2.0d0*(& 1384 & gxa(re, 1)*rank5(re, 4)+gxa(im, 1)*rank5(im, 4)+& 1385 & gxa(re, 7)*rank5(re,18)+gxa(im, 7)*rank5(im,18)+& 1386 & gxa(re,10)*rank5(re,20)+gxa(im,10)*rank5(im,20)+& 1387 & gxa(re, 2)*rank5(re,13)+gxa(im, 2)*rank5(im,13)+& 1388 & gxa(re, 3)*rank5(re,14)+gxa(im, 3)*rank5(im,14)+& 1389 & gxa(re, 5)*rank5(re, 8)+gxa(im, 5)*rank5(im, 8)+& 1390 & gxa(re, 6)*rank5(re, 9)+gxa(im, 6)*rank5(im, 9)+& 1391 & gxa(re, 8)*rank5(re,19)+gxa(im, 8)*rank5(im,19)+& 1392 & gxa(re, 9)*rank5(re,17)+gxa(im, 9)*rank5(im,17)+& 1393 & gxa(re, 4)*rank5(re,12)+gxa(im, 4)*rank5(im,12)) 1394 1395 !a=3, b=1 in rank2(a,b) --> maps to index 5 1396 rank2(5)=2.0d0*(& 1397 & gxa(re, 1)*rank5(re, 5)+gxa(im, 1)*rank5(im, 5)+& 1398 & gxa(re, 7)*rank5(re,13)+gxa(im, 7)*rank5(im,13)+& 1399 & gxa(re,10)*rank5(re,15)+gxa(im,10)*rank5(im,15)+& 1400 & gxa(re, 2)*rank5(re, 9)+gxa(im, 2)*rank5(im, 9)+& 1401 & gxa(re, 3)*rank5(re,10)+gxa(im, 3)*rank5(im,10)+& 1402 & gxa(re, 5)*rank5(re, 3)+gxa(im, 5)*rank5(im, 3)+& 1403 & gxa(re, 6)*rank5(re, 4)+gxa(im, 6)*rank5(im, 4)+& 1404 & gxa(re, 8)*rank5(re,14)+gxa(im, 8)*rank5(im,14)+& 1405 & gxa(re, 9)*rank5(re,12)+gxa(im, 9)*rank5(im,12)+& 1406 & gxa(re, 4)*rank5(re, 8)+gxa(im, 4)*rank5(im, 8)) 1407 1408 !a=2, b=1 in rank2(a,b) --> maps to index 6 1409 rank2(6)=2.0d0*(& 1410 & gxa(re, 1)*rank5(re, 6)+gxa(im, 1)*rank5(im, 6)+& 1411 & gxa(re, 7)*rank5(re,11)+gxa(im, 7)*rank5(im,11)+& 1412 & gxa(re,10)*rank5(re,14)+gxa(im,10)*rank5(im,14)+& 1413 & gxa(re, 2)*rank5(re, 7)+gxa(im, 2)*rank5(im, 7)+& 1414 & gxa(re, 3)*rank5(re, 8)+gxa(im, 3)*rank5(im, 8)+& 1415 & gxa(re, 5)*rank5(re, 4)+gxa(im, 5)*rank5(im, 4)+& 1416 & gxa(re, 6)*rank5(re, 2)+gxa(im, 6)*rank5(im, 2)+& 1417 & gxa(re, 8)*rank5(re,12)+gxa(im, 8)*rank5(im,12)+& 1418 & gxa(re, 9)*rank5(re,13)+gxa(im, 9)*rank5(im,13)+& 1419 & gxa(re, 4)*rank5(re, 9)+gxa(im, 4)*rank5(im, 9)) 1420 1421 end subroutine cont35
m_contract/metcon [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
metcon
FUNCTION
Carries out specialized metric tensor contractions needed for l=0,1,2,3 nonlocal Kleinman-Bylander pseudopotential operation. Full advantage is taken of the full permutational symmetry of these tensors.
INPUTS
rank=0,1,2, or 3 = rank of input tensor aa gmet(3,3)=metric tensor (array is symmetric but stored as 3x3) aa(2,(rank+1)*(rank+2)/2)=unique elements of complex input tensor
OUTPUT
bb(2,(rank+1)*(rank+2)/2)=unique elements of complex output tensor
NOTES
All tensors are stored in a compressed storage mode defined below; input and output conform to this scheme. When tensor elements occur repeatedly due to symmetry, the WEIGHT IS INCLUDED in the output tensor element to simplify later contractions with other tensors of the same rank and form, i.e. the next contraction is then simply a dot product over the unique elements. Definitions of the contractions: rank=0: bb = aa (simply copy the scalar value, no use of gmet) rank=1: bb(i)= $gmet(i,l) aa(l)$ (3 elements in, 3 elements out) rank=2: bb(i,j)= $[{3 \over 2} gmet(i,l) gmet(j,m) - {1 \over 2} gmet(i,j) gmet(l,m)] aa(l,m)$ (6 elements in, 6 elements out) rank=3: bb(i,j,k)= $[{5 \over 2} g(i,l) g(j,m) g(k,n) - {3 \over 2} g(i,j) g(l,m) g(k,n)] aa(l,m,n)$ (10 elements in, 10 elements out) In this rank 3 case, the second term is NOT symmetric in all permutations of i,j,k, but the final tensor b(ijk) may be symmetrized over all permutations because it will be contracted with a completely symmetric tensor. The compressed storage scheme is based on storing a symmetric 3x3 matrix as (1 . .) (6 2 .) (5 4 3) which leads to the following mappings for all ranks where the compressed storage index is to the right of the arrow: rank=0 1->1 (only a scalar) rank=1 1->1 2->2 3->3 (standard vector, no compression) rank=2 11->1 22->2 33->3 32->4 31->5 21->6 weights 1 1 1 2 2 2 rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10 weights 1 3 3 6 3 3 1 3 3 1
SOURCE
1481 subroutine metcon(rank,gmet,aa,bb) 1482 1483 !Arguments ------------------------------------ 1484 !scalars 1485 integer,intent(in) :: rank 1486 !arrays 1487 real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),gmet(3,3) 1488 real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2) 1489 1490 !Local variables------------------------------- 1491 !scalars 1492 integer :: ii,jj 1493 real(dp) :: scalar,tmpiii,tmpijk 1494 character(len=500) :: message 1495 !arrays 1496 real(dp) :: vector(3) 1497 1498 ! ************************************************************************* 1499 1500 !This statement function defines the l=3 contraction in 1501 !terms of the free indices of the contracted tensor (re and im) 1502 ! coniii(ii,i1,i2,i3)=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+& 1503 !& gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+& 1504 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10) 1505 ! conijk(ii,i1,i2,i3)=aa(ii,4)*& 1506 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+& 1507 !& gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+& 1508 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+& 1509 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+& 1510 !& gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+& 1511 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,3)) 1512 ! con(ii,i1,i2,i3)=coniii(ii,i1,i2,i3)+conijk(ii,i1,i2,i3)+& 1513 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+& 1514 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+& 1515 !& gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+& 1516 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+& 1517 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+& 1518 !& gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+& 1519 !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+& 1520 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+& 1521 !& gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+& 1522 !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+& 1523 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+& 1524 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+& 1525 !& (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+& 1526 !& gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+& 1527 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+& 1528 !& (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+& 1529 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+& 1530 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8) 1531 1532 !DEBUG 1533 !write(std_out,*)' metcon : enter ' 1534 !stop 1535 !ENDDEBUG 1536 if (rank==0) then 1537 ! Simply copy scalar, re and im 1538 bb(1,1)=aa(1,1) 1539 bb(2,1)=aa(2,1) 1540 1541 else if (rank==1) then 1542 ! Apply gmet to input vector, re and im 1543 do ii=1,2 1544 do jj=1,3 1545 bb(ii,jj)=gmet(jj,1)*aa(ii,1)+gmet(jj,2)*aa(ii,2)+gmet(jj,3)*aa(ii,3) 1546 end do 1547 end do 1548 1549 1550 else if (rank==2) then 1551 ! Apply rank2 expression, re and im 1552 do ii=1,2 1553 ! Carry out g(c,d)*aa(c,d) contraction to get scalar 1554 scalar=gmet(1,1)*aa(ii,1)+gmet(2,2)*aa(ii,2)+& 1555 & gmet(3,3)*aa(ii,3)+2.0d0*(gmet(2,1)*aa(ii,6)+& 1556 & gmet(3,1)*aa(ii,5)+gmet(3,2)*aa(ii,4) ) 1557 ! Write out components of contraction 1558 ! (1,1)->1 1559 bb(ii,1)=0.5d0*(3.0d0*(gmet(1,1)*gmet(1,1)*aa(ii,1)+& 1560 & gmet(1,2)*gmet(1,2)*aa(ii,2)+gmet(1,3)*gmet(1,3)*aa(ii,3)+& 1561 & (gmet(1,2)*gmet(1,1)+gmet(1,1)*gmet(1,2))*aa(ii,6)+& 1562 & (gmet(1,3)*gmet(1,1)+gmet(1,1)*gmet(1,3))*aa(ii,5)+& 1563 & (gmet(1,3)*gmet(1,2)+gmet(1,2)*gmet(1,3))*aa(ii,4) ) & 1564 & - gmet(1,1)*scalar) 1565 ! (2,2)->2 1566 bb(ii,2)=0.5d0*(3.0d0*(gmet(2,1)*gmet(2,1)*aa(ii,1)+& 1567 & gmet(2,2)*gmet(2,2)*aa(ii,2)+gmet(2,3)*gmet(2,3)*aa(ii,3)+& 1568 & (gmet(2,2)*gmet(2,1)+gmet(2,1)*gmet(2,2))*aa(ii,6)+& 1569 & (gmet(2,3)*gmet(2,1)+gmet(2,1)*gmet(2,3))*aa(ii,5)+& 1570 & (gmet(2,3)*gmet(2,2)+gmet(2,2)*gmet(2,3))*aa(ii,4) )& 1571 & - gmet(2,2)*scalar) 1572 ! (3,3)->3 1573 bb(ii,3)=0.5d0*(3.0d0*(gmet(3,1)*gmet(3,1)*aa(ii,1)+& 1574 & gmet(3,2)*gmet(3,2)*aa(ii,2)+gmet(3,3)*gmet(3,3)*aa(ii,3)+& 1575 & (gmet(3,2)*gmet(3,1)+gmet(3,1)*gmet(3,2))*aa(ii,6)+& 1576 & (gmet(3,3)*gmet(3,1)+gmet(3,1)*gmet(3,3))*aa(ii,5)+& 1577 & (gmet(3,3)*gmet(3,2)+gmet(3,2)*gmet(3,3))*aa(ii,4) )& 1578 & - gmet(3,3)*scalar) 1579 ! (3,2)->4 1580 bb(ii,4)=0.5d0*(3.0d0*(gmet(3,1)*gmet(2,1)*aa(ii,1)+& 1581 & gmet(3,2)*gmet(2,2)*aa(ii,2)+gmet(3,3)*gmet(2,3)*aa(ii,3)+& 1582 & (gmet(3,2)*gmet(2,1)+gmet(3,1)*gmet(2,2))*aa(ii,6)+& 1583 & (gmet(3,3)*gmet(2,1)+gmet(3,1)*gmet(2,3))*aa(ii,5)+& 1584 & (gmet(3,3)*gmet(2,2)+gmet(3,2)*gmet(2,3))*aa(ii,4) )& 1585 & - gmet(3,2)*scalar) 1586 ! (3,1)->5 1587 bb(ii,5)=0.5d0*(3.0d0*(gmet(3,1)*gmet(1,1)*aa(ii,1)+& 1588 & gmet(3,2)*gmet(1,2)*aa(ii,2)+gmet(3,3)*gmet(1,3)*aa(ii,3)+& 1589 & (gmet(3,2)*gmet(1,1)+gmet(3,1)*gmet(1,2))*aa(ii,6)+& 1590 & (gmet(3,3)*gmet(1,1)+gmet(3,1)*gmet(1,3))*aa(ii,5)+& 1591 & (gmet(3,3)*gmet(1,2)+gmet(3,2)*gmet(1,3))*aa(ii,4) )& 1592 & - gmet(3,1)*scalar) 1593 ! (2,1)->6 1594 bb(ii,6)=0.5d0*(3.0d0*(gmet(2,1)*gmet(1,1)*aa(ii,1)+& 1595 & gmet(2,2)*gmet(1,2)*aa(ii,2)+gmet(2,3)*gmet(1,3)*aa(ii,3)+& 1596 & (gmet(2,2)*gmet(1,1)+gmet(2,1)*gmet(1,2))*aa(ii,6)+& 1597 & (gmet(2,3)*gmet(1,1)+gmet(2,1)*gmet(1,3))*aa(ii,5)+& 1598 & (gmet(2,3)*gmet(1,2)+gmet(2,2)*gmet(1,3))*aa(ii,4) ) & 1599 & - gmet(2,1)*scalar) 1600 ! Include appropriate weights for multiplicity 1601 bb(ii,4)=2.d0*bb(ii,4) 1602 bb(ii,5)=2.d0*bb(ii,5) 1603 bb(ii,6)=2.d0*bb(ii,6) 1604 end do 1605 1606 else if (rank==3) then 1607 ! Apply rank2 expression, re and im 1608 do ii=1,2 1609 ! Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j) 1610 do jj=1,3 1611 tmpiii= gmet(1,1)*gmet(jj,1)*aa(ii,1)+& 1612 & gmet(2,2)*gmet(jj,2)*aa(ii,7)+& 1613 & gmet(3,3)*gmet(jj,3)*aa(ii,10) 1614 tmpijk= (gmet(1,2)*gmet(jj,3)+& 1615 & gmet(3,1)*gmet(jj,2)+& 1616 & gmet(2,3)*gmet(jj,1)+& 1617 & gmet(3,2)*gmet(jj,1)+& 1618 & gmet(1,3)*gmet(jj,2)+& 1619 & gmet(2,1)*gmet(jj,3)) *aa(ii,4) 1620 vector(jj)=tmpiii + tmpijk +& 1621 & (gmet(1,2)*gmet(jj,1)+& 1622 & gmet(2,1)*gmet(jj,1)+& 1623 & gmet(1,1)*gmet(jj,2)) *aa(ii,6)+& 1624 & (gmet(1,2)*gmet(jj,2)+& 1625 & gmet(2,1)*gmet(jj,2)+& 1626 & gmet(2,2)*gmet(jj,1)) *aa(ii,2)+& 1627 & (gmet(1,3)*gmet(jj,1)+& 1628 & gmet(3,1)*gmet(jj,1)+& 1629 & gmet(1,1)*gmet(jj,3)) *aa(ii,5)+& 1630 & (gmet(1,3)*gmet(jj,3)+& 1631 & gmet(3,1)*gmet(jj,3)+& 1632 & gmet(3,3)*gmet(jj,1)) *aa(ii,3)+& 1633 & (gmet(2,3)*gmet(jj,2)+& 1634 & gmet(3,2)*gmet(jj,2)+& 1635 & gmet(2,2)*gmet(jj,3)) *aa(ii,9)+& 1636 & (gmet(2,3)*gmet(jj,3)+& 1637 & gmet(3,2)*gmet(jj,3)+& 1638 & gmet(3,3)*gmet(jj,2)) *aa(ii,8) 1639 end do 1640 ! Write out components of contraction 1641 ! (111)->1 1642 bb(ii,1) =2.5d0*con_met(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1)) 1643 ! (221)->2 1644 bb(ii,2) =2.5d0*con_met(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+& 1645 & gmet(1,2)*vector(2)+gmet(2,2)*vector(1)) 1646 ! (331)->3 1647 bb(ii,3) =2.5d0*con_met(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+& 1648 & gmet(1,3)*vector(3)+gmet(3,3)*vector(1)) 1649 ! (321)->4 1650 bb(ii,4) =2.5d0*con_met(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+& 1651 & gmet(1,2)*vector(3)+gmet(3,2)*vector(1)) 1652 ! (311)->5 1653 bb(ii,5) =2.5d0*con_met(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+& 1654 & gmet(1,1)*vector(3)+gmet(3,1)*vector(1)) 1655 ! (211)->6 1656 bb(ii,6) =2.5d0*con_met(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+& 1657 & gmet(1,1)*vector(2)+gmet(2,1)*vector(1)) 1658 ! (222)->7 1659 bb(ii,7) =2.5d0*con_met(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2)) 1660 1661 ! (332)->8 1662 bb(ii,8) =2.5d0*con_met(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+& 1663 & gmet(2,3)*vector(3)+gmet(3,3)*vector(2)) 1664 ! (322)->9 1665 bb(ii,9) =2.5d0*con_met(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+& 1666 & gmet(2,2)*vector(3)+gmet(3,2)*vector(2)) 1667 ! (333)->10 1668 bb(ii,10)=2.5d0*con_met(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3)) 1669 ! Include appropriate weights for multiplicity 1670 bb(ii,2)=3.d0*bb(ii,2) 1671 bb(ii,3)=3.d0*bb(ii,3) 1672 bb(ii,4)=6.d0*bb(ii,4) 1673 bb(ii,5)=3.d0*bb(ii,5) 1674 bb(ii,6)=3.d0*bb(ii,6) 1675 bb(ii,8)=3.d0*bb(ii,8) 1676 bb(ii,9)=3.d0*bb(ii,9) 1677 end do 1678 1679 else 1680 write(message, '(a,i0,a,a,a)' )& 1681 & 'Input rank=',rank,' not allowed.',ch10,& 1682 & 'Possible values are 0,1,2,3 only.' 1683 ABI_BUG(message) 1684 end if 1685 1686 contains 1687 1688 function con_met(ii,i1,i2,i3) 1689 1690 real(dp) :: con_met 1691 integer :: ii,i1,i2,i3 1692 real(dp)::coniii,conijk 1693 1694 coniii=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+& 1695 & gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+& 1696 & gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10) 1697 conijk=aa(ii,4)*& 1698 & (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+& 1699 & gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+& 1700 & gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+& 1701 & gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+& 1702 & gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+& 1703 & gmet(i1,2)*gmet(i2,1)*gmet(i3,3)) 1704 con_met=coniii+conijk+& 1705 & (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+& 1706 & gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+& 1707 & gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+& 1708 & (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+& 1709 & gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+& 1710 & gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+& 1711 & (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+& 1712 & gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+& 1713 & gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+& 1714 & (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+& 1715 & gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+& 1716 & gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+& 1717 & (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+& 1718 & gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+& 1719 & gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+& 1720 & (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+& 1721 & gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+& 1722 & gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8) 1723 1724 end function con_met 1725 1726 end subroutine metcon
m_contract/metcon_so [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
metcon_so
FUNCTION
Carries out specialized metric tensor contractions needed for l=0,1,2,3 nonlocal Kleinman-Bylander pseudopotential operation in the spin-orbit case. Full advantage is taken of the full permutational symmetry of these tensors.
INPUTS
rank=0,1,2, or 3 = rank of input tensor aa gmet(3,3)=metric tensor (array is symmetric but stored as 3x3) amet(3,3)=real or imaginary part of one spin matrix element of the "spin metric" tensor aa(2,(rank+1)*(rank+2)/2)=unique elements of complex input tensor
OUTPUT
bb(2,(rank+1)*(rank+2)/2)=unique elements of complex output tensor
NOTES
All tensors are stored in a compressed storage mode defined below; input and output conform to this scheme. When tensor elements occur repeatedly due to symmetry, the WEIGHT IS INCLUDED in the output tensor element to simplify later contractions with other tensors of the same rank and form, i.e. the next contraction is then simply a dot product over the unique elements. Definitions of the contractions: rank=0: bb=0 rank=1: bb(i)= $amet(i,l) aa(l)$ (3 elements in, 3 elements out) rank=2: bb(i,j)= $[3 gmet(i,l) amet(j,m)] aa(l,m)$ (6 elements in, 6 elements out) rank=3: bb(i,j,k)= $[{15 \over 2} g(i,l) g(j,m) a(k,n) - {3 \over 2} g(i,j) g(l,m) a(k,n)] aa(l,m,n)$ (10 elements in, 10 elements out) In this rank 3 case, the second term is NOT symmetric in all permutations of i,j,k, but the final tensor b(ijk) may be symmetrized over all permutations because it will be contracted with a completely symmetric tensor. The compressed storage scheme is based on storing a symmetric 3x3 matrix as (1 . .) (6 2 .) (5 4 3) which leads to the following mappings for all ranks where the compressed storage index is to the right of the arrow: rank=0 1->1 (only a scalar) rank=1 1->1 2->2 3->3 (standard vector, no compression) rank=2 11->1 22->2 33->3 32->4 31->5 21->6 weights 1 1 1 2 2 2 rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10 weights 1 3 3 6 3 3 1 3 3 1
SOURCE
1789 subroutine metcon_so(rank,gmet,amet,aa,bb) 1790 1791 !Arguments ------------------------------------ 1792 !scalars 1793 integer,intent(in) :: rank 1794 !arrays 1795 real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),amet(3,3),gmet(3,3) 1796 real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2) 1797 1798 !Local variables------------------------------- 1799 !scalars 1800 integer :: ii,jj 1801 real(dp) :: tmpiii,tmpijk 1802 character(len=500) :: message 1803 !arrays 1804 real(dp) :: vector(3) 1805 1806 ! ************************************************************************* 1807 1808 if (rank==0) then 1809 ! Simply copy scalar, re and im 1810 bb(1,1)=0.d0 1811 bb(2,1)=0.d0 1812 ! 1813 else if (rank==1) then 1814 ! Apply gmet to input vector, re and im 1815 do ii=1,2 1816 do jj=1,3 1817 bb(ii,jj)=amet(jj,1)*aa(ii,1)+amet(jj,2)*aa(ii,2)+amet(jj,3)*aa(ii,3) 1818 end do 1819 end do 1820 1821 else if (rank==2) then 1822 ! Apply rank2 expression, re and im 1823 do ii=1,2 1824 ! Write out components of contraction 1825 ! (1,1)->1 1826 bb(ii,1)=3.0d0*(gmet(1,1)*amet(1,1)*aa(ii,1)+& 1827 & gmet(1,2)*amet(1,2)*aa(ii,2)+gmet(1,3)*amet(1,3)*aa(ii,3)+& 1828 & (gmet(1,2)*amet(1,1)+gmet(1,1)*amet(1,2))*aa(ii,6)+& 1829 & (gmet(1,3)*amet(1,1)+gmet(1,1)*amet(1,3))*aa(ii,5)+& 1830 & (gmet(1,3)*amet(1,2)+gmet(1,2)*amet(1,3))*aa(ii,4)) 1831 ! (2,2)->2 1832 bb(ii,2)=3.0d0*(gmet(2,1)*amet(2,1)*aa(ii,1)+& 1833 & gmet(2,2)*amet(2,2)*aa(ii,2)+gmet(2,3)*amet(2,3)*aa(ii,3)+& 1834 & (gmet(2,2)*amet(2,1)+gmet(2,1)*amet(2,2))*aa(ii,6)+& 1835 & (gmet(2,3)*amet(2,1)+gmet(2,1)*amet(2,3))*aa(ii,5)+& 1836 & (gmet(2,3)*amet(2,2)+gmet(2,2)*amet(2,3))*aa(ii,4) ) 1837 ! (3,3)->3 1838 bb(ii,3)=3.0d0*(gmet(3,1)*amet(3,1)*aa(ii,1)+& 1839 & gmet(3,2)*amet(3,2)*aa(ii,2)+gmet(3,3)*amet(3,3)*aa(ii,3)+& 1840 & (gmet(3,2)*amet(3,1)+gmet(3,1)*amet(3,2))*aa(ii,6)+& 1841 & (gmet(3,3)*amet(3,1)+gmet(3,1)*amet(3,3))*aa(ii,5)+& 1842 & (gmet(3,3)*amet(3,2)+gmet(3,2)*amet(3,3))*aa(ii,4) ) 1843 ! (3,2)->4 1844 bb(ii,4)=3.0d0*(gmet(3,1)*amet(2,1)*aa(ii,1)+& 1845 & gmet(3,2)*amet(2,2)*aa(ii,2)+gmet(3,3)*amet(2,3)*aa(ii,3)+& 1846 & (gmet(3,2)*amet(2,1)+gmet(3,1)*amet(2,2))*aa(ii,6)+& 1847 & (gmet(3,3)*amet(2,1)+gmet(3,1)*amet(2,3))*aa(ii,5)+& 1848 & (gmet(3,3)*amet(2,2)+gmet(3,2)*amet(2,3))*aa(ii,4) ) 1849 bb(ii,4)=bb(ii,4)+3.0d0*(amet(3,1)*gmet(2,1)*aa(ii,1)+& 1850 & amet(3,2)*gmet(2,2)*aa(ii,2)+amet(3,3)*gmet(2,3)*aa(ii,3)+& 1851 & (amet(3,2)*gmet(2,1)+amet(3,1)*gmet(2,2))*aa(ii,6)+& 1852 & (amet(3,3)*gmet(2,1)+amet(3,1)*gmet(2,3))*aa(ii,5)+& 1853 & (amet(3,3)*gmet(2,2)+amet(3,2)*gmet(2,3))*aa(ii,4) ) 1854 ! (3,1)->5 1855 bb(ii,5)=3.0d0*(gmet(3,1)*amet(1,1)*aa(ii,1)+& 1856 & gmet(3,2)*amet(1,2)*aa(ii,2)+gmet(3,3)*amet(1,3)*aa(ii,3)+& 1857 & (gmet(3,2)*amet(1,1)+gmet(3,1)*amet(1,2))*aa(ii,6)+& 1858 & (gmet(3,3)*amet(1,1)+gmet(3,1)*amet(1,3))*aa(ii,5)+& 1859 & (gmet(3,3)*amet(1,2)+gmet(3,2)*amet(1,3))*aa(ii,4) ) 1860 bb(ii,5)=bb(ii,5)+3.0d0*(amet(3,1)*gmet(1,1)*aa(ii,1)+& 1861 & amet(3,2)*gmet(1,2)*aa(ii,2)+amet(3,3)*gmet(1,3)*aa(ii,3)+& 1862 & (amet(3,2)*gmet(1,1)+amet(3,1)*gmet(1,2))*aa(ii,6)+& 1863 & (amet(3,3)*gmet(1,1)+amet(3,1)*gmet(1,3))*aa(ii,5)+& 1864 & (amet(3,3)*gmet(1,2)+amet(3,2)*gmet(1,3))*aa(ii,4) ) 1865 ! (2,1)->6 1866 bb(ii,6)=3.0d0*(gmet(2,1)*amet(1,1)*aa(ii,1)+& 1867 & gmet(2,2)*amet(1,2)*aa(ii,2)+gmet(2,3)*amet(1,3)*aa(ii,3)+& 1868 & (gmet(2,2)*amet(1,1)+gmet(2,1)*amet(1,2))*aa(ii,6)+& 1869 & (gmet(2,3)*amet(1,1)+gmet(2,1)*amet(1,3))*aa(ii,5)+& 1870 & (gmet(2,3)*amet(1,2)+gmet(2,2)*amet(1,3))*aa(ii,4) ) 1871 bb(ii,6)=bb(ii,6)+3.0d0*(amet(2,1)*gmet(1,1)*aa(ii,1)+& 1872 & amet(2,2)*gmet(1,2)*aa(ii,2)+amet(2,3)*gmet(1,3)*aa(ii,3)+& 1873 & (amet(2,2)*gmet(1,1)+amet(2,1)*gmet(1,2))*aa(ii,6)+& 1874 & (amet(2,3)*gmet(1,1)+amet(2,1)*gmet(1,3))*aa(ii,5)+& 1875 & (amet(2,3)*gmet(1,2)+amet(2,2)*gmet(1,3))*aa(ii,4) ) 1876 ! Include appropriate weights for multiplicity 1877 bb(ii,4)=bb(ii,4) 1878 bb(ii,5)=bb(ii,5) 1879 bb(ii,6)=bb(ii,6) 1880 end do 1881 1882 else if (rank==3) then 1883 ! Apply rank2 expression, re and im 1884 do ii=1,2 1885 ! Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j) 1886 do jj=1,3 1887 tmpiii= gmet(1,1)*amet(jj,1)*aa(ii,1)+& 1888 & gmet(2,2)*amet(jj,2)*aa(ii,7)+& 1889 & gmet(3,3)*amet(jj,3)*aa(ii,10) 1890 tmpijk= (gmet(1,2)*amet(jj,3)+& 1891 & gmet(3,1)*amet(jj,2)+& 1892 & gmet(2,3)*amet(jj,1)+& 1893 & gmet(3,2)*amet(jj,1)+& 1894 & gmet(1,3)*amet(jj,2)+& 1895 & gmet(2,1)*amet(jj,3)) *aa(ii,4) 1896 vector(jj)=tmpiii + tmpijk +& 1897 & (gmet(1,2)*amet(jj,1)+& 1898 & gmet(2,1)*amet(jj,1)+& 1899 & gmet(1,1)*amet(jj,2)) *aa(ii,6)+& 1900 & (gmet(1,2)*amet(jj,2)+& 1901 & gmet(2,1)*amet(jj,2)+& 1902 & gmet(2,2)*amet(jj,1)) *aa(ii,2)+& 1903 & (gmet(1,3)*amet(jj,1)+& 1904 & gmet(3,1)*amet(jj,1)+& 1905 & gmet(1,1)*amet(jj,3)) *aa(ii,5)+& 1906 & (gmet(1,3)*amet(jj,3)+& 1907 & gmet(3,1)*amet(jj,3)+& 1908 & gmet(3,3)*amet(jj,1)) *aa(ii,3)+& 1909 & (gmet(2,3)*amet(jj,2)+& 1910 & gmet(3,2)*amet(jj,2)+& 1911 & gmet(2,2)*amet(jj,3)) *aa(ii,9)+& 1912 & (gmet(2,3)*amet(jj,3)+& 1913 & gmet(3,2)*amet(jj,3)+& 1914 & gmet(3,3)*amet(jj,2)) *aa(ii,8) 1915 end do 1916 ! Write out components of contraction 1917 ! (111)->1 1918 bb(ii,1) =7.5d0*con_metso(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1)) 1919 ! (221)->2 1920 bb(ii,2) =7.5d0*con_metso(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+& 1921 & gmet(1,2)*vector(2)+gmet(2,2)*vector(1)) 1922 ! (331)->3 1923 bb(ii,3) =7.5d0*con_metso(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+& 1924 & gmet(1,3)*vector(3)+gmet(3,3)*vector(1)) 1925 ! (321)->4 1926 bb(ii,4) =7.5d0*con_metso(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+& 1927 & gmet(1,2)*vector(3)+gmet(3,2)*vector(1)) 1928 ! (311)->5 1929 bb(ii,5) =7.5d0*con_metso(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+& 1930 & gmet(1,1)*vector(3)+gmet(3,1)*vector(1)) 1931 ! (211)->6 1932 bb(ii,6) =7.5d0*con_metso(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+& 1933 & gmet(1,1)*vector(2)+gmet(2,1)*vector(1)) 1934 ! (222)->7 1935 bb(ii,7) =7.5d0*con_metso(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2)) 1936 1937 ! (332)->8 1938 bb(ii,8) =7.5d0*con_metso(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+& 1939 & gmet(2,3)*vector(3)+gmet(3,3)*vector(2)) 1940 ! (322)->9 1941 bb(ii,9) =7.5d0*con_metso(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+& 1942 & gmet(2,2)*vector(3)+gmet(3,2)*vector(2)) 1943 ! (333)->10 1944 bb(ii,10)=7.5d0*con_metso(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3)) 1945 ! Include appropriate weights for multiplicity 1946 bb(ii,2)=3.d0*bb(ii,2) 1947 bb(ii,3)=3.d0*bb(ii,3) 1948 bb(ii,4)=6.d0*bb(ii,4) 1949 bb(ii,5)=3.d0*bb(ii,5) 1950 bb(ii,6)=3.d0*bb(ii,6) 1951 bb(ii,8)=3.d0*bb(ii,8) 1952 bb(ii,9)=3.d0*bb(ii,9) 1953 end do 1954 1955 else 1956 write(message, '(a,i0,a,a,a)' )& 1957 & 'Input rank=',rank,' not allowed.',ch10,& 1958 & 'Possible values are 0,1,2,3 only.' 1959 ABI_BUG(message) 1960 end if 1961 1962 contains 1963 1964 !This function defines the l=3 contraction in 1965 !terms of the free indices of the contracted tensor (re and im) 1966 1967 function cona_metso(ii,i1,i2,i3) 1968 1969 real(dp) :: cona_metso 1970 integer,intent(in) :: ii,i1,i2,i3 1971 real(dp) :: coniii, conijk 1972 1973 coniii=gmet(i1,1)*gmet(i2,1)*amet(i3,1)*aa(ii,1)+& 1974 & gmet(i1,2)*gmet(i2,2)*amet(i3,2)*aa(ii,7)+& 1975 & gmet(i1,3)*gmet(i2,3)*amet(i3,3)*aa(ii,10) 1976 conijk=aa(ii,4)*& 1977 & (gmet(i1,1)*gmet(i2,2)*amet(i3,3)+& 1978 & gmet(i1,2)*gmet(i2,3)*amet(i3,1)+& 1979 & gmet(i1,3)*gmet(i2,1)*amet(i3,2)+& 1980 & gmet(i1,3)*gmet(i2,2)*amet(i3,1)+& 1981 & gmet(i1,1)*gmet(i2,3)*amet(i3,2)+& 1982 & gmet(i1,2)*gmet(i2,1)*amet(i3,3)) 1983 cona_metso=coniii+conijk+& 1984 & (gmet(i1,1)*gmet(i2,2)*amet(i3,1)+& 1985 & gmet(i1,2)*gmet(i2,1)*amet(i3,1)+& 1986 & gmet(i1,1)*gmet(i2,1)*amet(i3,2))*aa(ii,6)+& 1987 & (gmet(i1,1)*gmet(i2,2)*amet(i3,2)+& 1988 & gmet(i1,2)*gmet(i2,1)*amet(i3,2)+& 1989 & gmet(i1,2)*gmet(i2,2)*amet(i3,1))*aa(ii,2)+& 1990 & (gmet(i1,1)*gmet(i2,3)*amet(i3,1)+& 1991 & gmet(i1,3)*gmet(i2,1)*amet(i3,1)+& 1992 & gmet(i1,1)*gmet(i2,1)*amet(i3,3))*aa(ii,5)+& 1993 & (gmet(i1,1)*gmet(i2,3)*amet(i3,3)+& 1994 & gmet(i1,3)*gmet(i2,1)*amet(i3,3)+& 1995 & gmet(i1,3)*gmet(i2,3)*amet(i3,1))*aa(ii,3)+& 1996 & (gmet(i1,2)*gmet(i2,2)*amet(i3,3)+& 1997 & gmet(i1,2)*gmet(i2,3)*amet(i3,2)+& 1998 & gmet(i1,3)*gmet(i2,2)*amet(i3,2))*aa(ii,9)+& 1999 & (gmet(i1,2)*gmet(i2,3)*amet(i3,3)+& 2000 & gmet(i1,3)*gmet(i2,2)*amet(i3,3)+& 2001 & gmet(i1,3)*gmet(i2,3)*amet(i3,2))*aa(ii,8) 2002 end function cona_metso 2003 2004 2005 function con_metso(ii,i1,i2,i3) 2006 2007 real(dp) :: con_metso 2008 integer,intent(in) :: ii,i1,i2,i3 2009 2010 con_metso=(cona_metso(ii,i3,i1,i2)+cona_metso(ii,i2,i3,i1)+cona_metso(ii,i1,i2,i3))/3.d0 2011 2012 end function con_metso 2013 2014 end subroutine metcon_so
m_contract/metric_so [ Functions ]
[ Top ] [ m_contract ] [ Functions ]
NAME
metric_so
FUNCTION
Computes Pauli matrices and antisymmetric tensor needed for spin-orbit.
INPUTS
gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
OUTPUT
amet(2,3,3,2,2)=the antisymmetric tensor A(Re/Im,y,y'',s,s'') pauli(2,2,2,3)=Pauli matrixes
NOTES
(the preprocessing based on CPP does not allow isolated quotes, so that they have been doubled in the following latex equations) $2 Pauli(s,s'',1) & = & 0 & 1 \nonumber & & 1 & 0 \nonumber 2 Pauli(s,s'',2) & = & 0 & -i \nonumber & & i & 0 \nonumber 2 Pauli(s,s'',3) & = & 1 & 0 \nonumber & & 0 & -1 \end{eqnarray} }} $Amet(y,y'',s,s'') & = & -i Pauli(s,s'',n) E(n,m,m'') Gprimd(m,y) Gprimd(m'',y'') \end{eqnarray} }} E(n,m,m''): full antisymetric tensor s,s'': spin indices (1..2) y,y'',m,m'',n: metric indices (1..3) a,b: strain indices (1..3) Amet and Pauli are complex
SOURCE
2056 subroutine metric_so(amet,gprimd,pauli) 2057 2058 !Arguments ------------------------------------ 2059 !arrays 2060 real(dp),intent(in) :: gprimd(3,3) 2061 real(dp),intent(out) :: amet(2,3,3,2,2),pauli(2,2,2,3) 2062 2063 !Local variables------------------------------- 2064 !scalars 2065 integer :: iy1,iy2,m1,m2,n 2066 !arrays 2067 real(dp) :: buffer1(3,3,2,2) !,buffer2(3,3,3,3,2,2) 2068 2069 ! ********************************************************************** 2070 2071 !Fill in Pauli matrices and make them spin matrices: 2072 !Pauli(Re/Im,up,down,coord): 2073 pauli(:,:,:,:)=0.d0 2074 pauli(1,1,2,1)= 1.d0;pauli(1,2,1,1)= 1.d0 2075 pauli(2,1,2,2)=-1.d0;pauli(2,2,1,2)= 1.d0 2076 pauli(1,1,1,3)= 1.d0;pauli(1,2,2,3)=-1.d0 2077 pauli(:,:,:,:)= 0.5d0*pauli(:,:,:,:) 2078 2079 !Construct the antisymmetric tensor: 2080 amet(:,:,:,:,:)=0.d0 2081 do iy2=1,3 2082 do iy1=1,3 2083 do n=1,3 2084 m1=mod(n ,3)+1 ! n,m1,m2 is an even permutation 2085 m2=mod(m1,3)+1 2086 amet(1:2,iy1,iy2,1:2,1:2) = amet(:,iy1,iy2,:,:) & 2087 & + pauli(:,:,:,n) & 2088 & *(gprimd(m1,iy1)*gprimd(m2,iy2) & 2089 & -gprimd(m2,iy1)*gprimd(m1,iy2)) 2090 end do 2091 end do 2092 end do 2093 !amet= -i amet 2094 buffer1(:,:,:,:)=-amet(1,:,:,:,:) 2095 amet(1,:,:,:,:) = amet(2,:,:,:,:) 2096 amet(2,:,:,:,:) =buffer1(:,:,:,:) 2097 2098 !DEBUG 2099 !!Eventually construct the gradients of Amet wrt strains: 2100 !! DAmet(y,y'',a,b,s,s'')= d[Amet(y,y'',s,s'')]/d[strain(a,b)] 2101 !! -i Pauli(s,s'',n)* 2102 !! ( E(n,a,m'')*Gprimd(b,y )*Gprimd(m'',y'') 2103 !! +E(n,m,a )*Gprimd(b,y'')*Gprimd(m ,y ) ) 2104 !damet(:,:,:,:,:,:,:)=0.d0 2105 !do ib=1,3 2106 !do ia=1,3 2107 !m1=mod(ia,3)+1 ! ia,m1,m2 is an even permutation 2108 !m2=mod(m1,3)+1 2109 !do iy2=1,3 2110 !do iy1=1,3 2111 !damet(:,iy1,iy2,ia,ib,:,:) = damet(:,iy1,iy2,ia,ib,:,:) & 2112 !& + (pauli(:,:,:,m2)*gprimd(m1,iy2) & 2113 !-pauli(:,:,:,m1)*gprimd(m2,iy2))*gprimd(ib,iy1) & 2114 !& + (pauli(:,:,:,m1)*gprimd(m2,iy1) & 2115 !-pauli(:,:,:,m2)*gprimd(m1,iy1))*gprimd(ib,iy2) 2116 !end do 2117 !end do 2118 !end do 2119 !end do 2120 !! damet= i damet 2121 !buffer2(:,:,:,:,:,:)= damet(1,:,:,:,:,:,:) 2122 !damet(1,:,:,:,:,:,:)= -damet(2,:,:,:,:,:,:) 2123 !damet(2,:,:,:,:,:,:)=buffer2(:,:,:,:,:,:) 2124 !! Symetrize damet(:,:,:,a,b,:,:) 2125 !damet(:,:,:,1,2,:,:)=0.5d0*(damet(:,:,:,1,2,:,:)+damet(:,:,:,2,1,:,:)) 2126 !damet(:,:,:,1,3,:,:)=0.5d0*(damet(:,:,:,1,3,:,:)+damet(:,:,:,3,1,:,:)) 2127 !damet(:,:,:,2,3,:,:)=0.5d0*(damet(:,:,:,2,3,:,:)+damet(:,:,:,3,2,:,:)) 2128 !damet(:,:,:,2,1,:,:)=damet(:,:,:,1,2,:,:) 2129 !damet(:,:,:,3,1,:,:)=damet(:,:,:,1,3,:,:) 2130 !damet(:,:,:,3,2,:,:)=damet(:,:,:,2,3,:,:) 2131 !ENDDEBUG 2132 2133 end subroutine metric_so