TABLE OF CONTENTS
ABINIT/metcon_so [ 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.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
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
75 #if defined HAVE_CONFIG_H 76 #include "config.h" 77 #endif 78 79 #include "abi_common.h" 80 81 82 subroutine metcon_so(rank,gmet,amet,aa,bb) 83 84 use defs_basis 85 use m_errors 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'metcon_so' 91 !End of the abilint section 92 93 implicit none 94 95 !Arguments ------------------------------------ 96 !scalars 97 integer,intent(in) :: rank 98 !arrays 99 real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),amet(3,3),gmet(3,3) 100 real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2) 101 102 !Local variables------------------------------- 103 !scalars 104 integer :: ii,jj 105 real(dp) :: tmpiii,tmpijk 106 character(len=500) :: message 107 !arrays 108 real(dp) :: vector(3) 109 110 ! ************************************************************************* 111 112 if (rank==0) then 113 ! Simply copy scalar, re and im 114 bb(1,1)=0.d0 115 bb(2,1)=0.d0 116 ! 117 else if (rank==1) then 118 ! Apply gmet to input vector, re and im 119 do ii=1,2 120 do jj=1,3 121 bb(ii,jj)=amet(jj,1)*aa(ii,1)+amet(jj,2)*aa(ii,2)+amet(jj,3)*aa(ii,3) 122 end do 123 end do 124 125 else if (rank==2) then 126 ! Apply rank2 expression, re and im 127 do ii=1,2 128 ! Write out components of contraction 129 ! (1,1)->1 130 bb(ii,1)=3.0d0*(gmet(1,1)*amet(1,1)*aa(ii,1)+& 131 & gmet(1,2)*amet(1,2)*aa(ii,2)+gmet(1,3)*amet(1,3)*aa(ii,3)+& 132 & (gmet(1,2)*amet(1,1)+gmet(1,1)*amet(1,2))*aa(ii,6)+& 133 & (gmet(1,3)*amet(1,1)+gmet(1,1)*amet(1,3))*aa(ii,5)+& 134 & (gmet(1,3)*amet(1,2)+gmet(1,2)*amet(1,3))*aa(ii,4)) 135 ! (2,2)->2 136 bb(ii,2)=3.0d0*(gmet(2,1)*amet(2,1)*aa(ii,1)+& 137 & gmet(2,2)*amet(2,2)*aa(ii,2)+gmet(2,3)*amet(2,3)*aa(ii,3)+& 138 & (gmet(2,2)*amet(2,1)+gmet(2,1)*amet(2,2))*aa(ii,6)+& 139 & (gmet(2,3)*amet(2,1)+gmet(2,1)*amet(2,3))*aa(ii,5)+& 140 & (gmet(2,3)*amet(2,2)+gmet(2,2)*amet(2,3))*aa(ii,4) ) 141 ! (3,3)->3 142 bb(ii,3)=3.0d0*(gmet(3,1)*amet(3,1)*aa(ii,1)+& 143 & gmet(3,2)*amet(3,2)*aa(ii,2)+gmet(3,3)*amet(3,3)*aa(ii,3)+& 144 & (gmet(3,2)*amet(3,1)+gmet(3,1)*amet(3,2))*aa(ii,6)+& 145 & (gmet(3,3)*amet(3,1)+gmet(3,1)*amet(3,3))*aa(ii,5)+& 146 & (gmet(3,3)*amet(3,2)+gmet(3,2)*amet(3,3))*aa(ii,4) ) 147 ! (3,2)->4 148 bb(ii,4)=3.0d0*(gmet(3,1)*amet(2,1)*aa(ii,1)+& 149 & gmet(3,2)*amet(2,2)*aa(ii,2)+gmet(3,3)*amet(2,3)*aa(ii,3)+& 150 & (gmet(3,2)*amet(2,1)+gmet(3,1)*amet(2,2))*aa(ii,6)+& 151 & (gmet(3,3)*amet(2,1)+gmet(3,1)*amet(2,3))*aa(ii,5)+& 152 & (gmet(3,3)*amet(2,2)+gmet(3,2)*amet(2,3))*aa(ii,4) ) 153 bb(ii,4)=bb(ii,4)+3.0d0*(amet(3,1)*gmet(2,1)*aa(ii,1)+& 154 & amet(3,2)*gmet(2,2)*aa(ii,2)+amet(3,3)*gmet(2,3)*aa(ii,3)+& 155 & (amet(3,2)*gmet(2,1)+amet(3,1)*gmet(2,2))*aa(ii,6)+& 156 & (amet(3,3)*gmet(2,1)+amet(3,1)*gmet(2,3))*aa(ii,5)+& 157 & (amet(3,3)*gmet(2,2)+amet(3,2)*gmet(2,3))*aa(ii,4) ) 158 ! (3,1)->5 159 bb(ii,5)=3.0d0*(gmet(3,1)*amet(1,1)*aa(ii,1)+& 160 & gmet(3,2)*amet(1,2)*aa(ii,2)+gmet(3,3)*amet(1,3)*aa(ii,3)+& 161 & (gmet(3,2)*amet(1,1)+gmet(3,1)*amet(1,2))*aa(ii,6)+& 162 & (gmet(3,3)*amet(1,1)+gmet(3,1)*amet(1,3))*aa(ii,5)+& 163 & (gmet(3,3)*amet(1,2)+gmet(3,2)*amet(1,3))*aa(ii,4) ) 164 bb(ii,5)=bb(ii,5)+3.0d0*(amet(3,1)*gmet(1,1)*aa(ii,1)+& 165 & amet(3,2)*gmet(1,2)*aa(ii,2)+amet(3,3)*gmet(1,3)*aa(ii,3)+& 166 & (amet(3,2)*gmet(1,1)+amet(3,1)*gmet(1,2))*aa(ii,6)+& 167 & (amet(3,3)*gmet(1,1)+amet(3,1)*gmet(1,3))*aa(ii,5)+& 168 & (amet(3,3)*gmet(1,2)+amet(3,2)*gmet(1,3))*aa(ii,4) ) 169 ! (2,1)->6 170 bb(ii,6)=3.0d0*(gmet(2,1)*amet(1,1)*aa(ii,1)+& 171 & gmet(2,2)*amet(1,2)*aa(ii,2)+gmet(2,3)*amet(1,3)*aa(ii,3)+& 172 & (gmet(2,2)*amet(1,1)+gmet(2,1)*amet(1,2))*aa(ii,6)+& 173 & (gmet(2,3)*amet(1,1)+gmet(2,1)*amet(1,3))*aa(ii,5)+& 174 & (gmet(2,3)*amet(1,2)+gmet(2,2)*amet(1,3))*aa(ii,4) ) 175 bb(ii,6)=bb(ii,6)+3.0d0*(amet(2,1)*gmet(1,1)*aa(ii,1)+& 176 & amet(2,2)*gmet(1,2)*aa(ii,2)+amet(2,3)*gmet(1,3)*aa(ii,3)+& 177 & (amet(2,2)*gmet(1,1)+amet(2,1)*gmet(1,2))*aa(ii,6)+& 178 & (amet(2,3)*gmet(1,1)+amet(2,1)*gmet(1,3))*aa(ii,5)+& 179 & (amet(2,3)*gmet(1,2)+amet(2,2)*gmet(1,3))*aa(ii,4) ) 180 ! Include appropriate weights for multiplicity 181 bb(ii,4)=bb(ii,4) 182 bb(ii,5)=bb(ii,5) 183 bb(ii,6)=bb(ii,6) 184 end do 185 186 else if (rank==3) then 187 ! Apply rank2 expression, re and im 188 do ii=1,2 189 ! Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j) 190 do jj=1,3 191 tmpiii= gmet(1,1)*amet(jj,1)*aa(ii,1)+& 192 & gmet(2,2)*amet(jj,2)*aa(ii,7)+& 193 & gmet(3,3)*amet(jj,3)*aa(ii,10) 194 tmpijk= (gmet(1,2)*amet(jj,3)+& 195 & gmet(3,1)*amet(jj,2)+& 196 & gmet(2,3)*amet(jj,1)+& 197 & gmet(3,2)*amet(jj,1)+& 198 & gmet(1,3)*amet(jj,2)+& 199 & gmet(2,1)*amet(jj,3)) *aa(ii,4) 200 vector(jj)=tmpiii + tmpijk +& 201 & (gmet(1,2)*amet(jj,1)+& 202 & gmet(2,1)*amet(jj,1)+& 203 & gmet(1,1)*amet(jj,2)) *aa(ii,6)+& 204 & (gmet(1,2)*amet(jj,2)+& 205 & gmet(2,1)*amet(jj,2)+& 206 & gmet(2,2)*amet(jj,1)) *aa(ii,2)+& 207 & (gmet(1,3)*amet(jj,1)+& 208 & gmet(3,1)*amet(jj,1)+& 209 & gmet(1,1)*amet(jj,3)) *aa(ii,5)+& 210 & (gmet(1,3)*amet(jj,3)+& 211 & gmet(3,1)*amet(jj,3)+& 212 & gmet(3,3)*amet(jj,1)) *aa(ii,3)+& 213 & (gmet(2,3)*amet(jj,2)+& 214 & gmet(3,2)*amet(jj,2)+& 215 & gmet(2,2)*amet(jj,3)) *aa(ii,9)+& 216 & (gmet(2,3)*amet(jj,3)+& 217 & gmet(3,2)*amet(jj,3)+& 218 & gmet(3,3)*amet(jj,2)) *aa(ii,8) 219 end do 220 ! Write out components of contraction 221 ! (111)->1 222 bb(ii,1) =7.5d0*con_metso(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1)) 223 ! (221)->2 224 bb(ii,2) =7.5d0*con_metso(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+& 225 & gmet(1,2)*vector(2)+gmet(2,2)*vector(1)) 226 ! (331)->3 227 bb(ii,3) =7.5d0*con_metso(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+& 228 & gmet(1,3)*vector(3)+gmet(3,3)*vector(1)) 229 ! (321)->4 230 bb(ii,4) =7.5d0*con_metso(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+& 231 & gmet(1,2)*vector(3)+gmet(3,2)*vector(1)) 232 ! (311)->5 233 bb(ii,5) =7.5d0*con_metso(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+& 234 & gmet(1,1)*vector(3)+gmet(3,1)*vector(1)) 235 ! (211)->6 236 bb(ii,6) =7.5d0*con_metso(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+& 237 & gmet(1,1)*vector(2)+gmet(2,1)*vector(1)) 238 ! (222)->7 239 bb(ii,7) =7.5d0*con_metso(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2)) 240 241 ! (332)->8 242 bb(ii,8) =7.5d0*con_metso(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+& 243 & gmet(2,3)*vector(3)+gmet(3,3)*vector(2)) 244 ! (322)->9 245 bb(ii,9) =7.5d0*con_metso(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+& 246 & gmet(2,2)*vector(3)+gmet(3,2)*vector(2)) 247 ! (333)->10 248 bb(ii,10)=7.5d0*con_metso(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3)) 249 ! Include appropriate weights for multiplicity 250 bb(ii,2)=3.d0*bb(ii,2) 251 bb(ii,3)=3.d0*bb(ii,3) 252 bb(ii,4)=6.d0*bb(ii,4) 253 bb(ii,5)=3.d0*bb(ii,5) 254 bb(ii,6)=3.d0*bb(ii,6) 255 bb(ii,8)=3.d0*bb(ii,8) 256 bb(ii,9)=3.d0*bb(ii,9) 257 end do 258 259 else 260 write(message, '(a,i0,a,a,a)' )& 261 & 'Input rank=',rank,' not allowed.',ch10,& 262 & 'Possible values are 0,1,2,3 only.' 263 MSG_BUG(message) 264 end if 265 266 contains 267 268 !This function defines the l=3 contraction in 269 !terms of the free indices of the contracted tensor (re and im) 270 271 function cona_metso(ii,i1,i2,i3) 272 273 274 !This section has been created automatically by the script Abilint (TD). 275 !Do not modify the following lines by hand. 276 #undef ABI_FUNC 277 #define ABI_FUNC 'cona_metso' 278 !End of the abilint section 279 280 real(dp) :: cona_metso 281 integer,intent(in) :: ii,i1,i2,i3 282 real(dp) :: coniii, conijk 283 284 coniii=gmet(i1,1)*gmet(i2,1)*amet(i3,1)*aa(ii,1)+& 285 & gmet(i1,2)*gmet(i2,2)*amet(i3,2)*aa(ii,7)+& 286 & gmet(i1,3)*gmet(i2,3)*amet(i3,3)*aa(ii,10) 287 conijk=aa(ii,4)*& 288 & (gmet(i1,1)*gmet(i2,2)*amet(i3,3)+& 289 & gmet(i1,2)*gmet(i2,3)*amet(i3,1)+& 290 & gmet(i1,3)*gmet(i2,1)*amet(i3,2)+& 291 & gmet(i1,3)*gmet(i2,2)*amet(i3,1)+& 292 & gmet(i1,1)*gmet(i2,3)*amet(i3,2)+& 293 & gmet(i1,2)*gmet(i2,1)*amet(i3,3)) 294 cona_metso=coniii+conijk+& 295 & (gmet(i1,1)*gmet(i2,2)*amet(i3,1)+& 296 & gmet(i1,2)*gmet(i2,1)*amet(i3,1)+& 297 & gmet(i1,1)*gmet(i2,1)*amet(i3,2))*aa(ii,6)+& 298 & (gmet(i1,1)*gmet(i2,2)*amet(i3,2)+& 299 & gmet(i1,2)*gmet(i2,1)*amet(i3,2)+& 300 & gmet(i1,2)*gmet(i2,2)*amet(i3,1))*aa(ii,2)+& 301 & (gmet(i1,1)*gmet(i2,3)*amet(i3,1)+& 302 & gmet(i1,3)*gmet(i2,1)*amet(i3,1)+& 303 & gmet(i1,1)*gmet(i2,1)*amet(i3,3))*aa(ii,5)+& 304 & (gmet(i1,1)*gmet(i2,3)*amet(i3,3)+& 305 & gmet(i1,3)*gmet(i2,1)*amet(i3,3)+& 306 & gmet(i1,3)*gmet(i2,3)*amet(i3,1))*aa(ii,3)+& 307 & (gmet(i1,2)*gmet(i2,2)*amet(i3,3)+& 308 & gmet(i1,2)*gmet(i2,3)*amet(i3,2)+& 309 & gmet(i1,3)*gmet(i2,2)*amet(i3,2))*aa(ii,9)+& 310 & (gmet(i1,2)*gmet(i2,3)*amet(i3,3)+& 311 & gmet(i1,3)*gmet(i2,2)*amet(i3,3)+& 312 & gmet(i1,3)*gmet(i2,3)*amet(i3,2))*aa(ii,8) 313 end function cona_metso 314 315 316 function con_metso(ii,i1,i2,i3) 317 318 319 !This section has been created automatically by the script Abilint (TD). 320 !Do not modify the following lines by hand. 321 #undef ABI_FUNC 322 #define ABI_FUNC 'con_metso' 323 !End of the abilint section 324 325 real(dp) :: con_metso 326 integer,intent(in) :: ii,i1,i2,i3 327 328 con_metso=(cona_metso(ii,i3,i1,i2)+cona_metso(ii,i2,i3,i1)+cona_metso(ii,i1,i2,i3))/3.d0 329 330 end function con_metso 331 332 end subroutine metcon_so