TABLE OF CONTENTS
ABINIT/metcon [ 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.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR) 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) aa(2,(rank+1)*(rank+2)/2)=unique elements of complex input tensor
OUTPUT
bb(2,(rank+1)*(rank+2)/2)=unique elements of complex output tensor
NOTES
All tensors are stored in a compressed storage mode defined below;
input and output conform to this scheme.
When tensor elements occur repeatedly due to symmetry, the
WEIGHT IS INCLUDED in the output tensor element to simplify later
contractions with other tensors of the same rank and form, i.e. the
next contraction is then simply a dot product over the unique elements.
Definitions of the contractions:
rank=0: bb = aa (simply copy the scalar value, no use of gmet)
rank=1: bb(i)= $gmet(i,l) aa(l)$ (3 elements in, 3 elements out)
rank=2: bb(i,j)= $[{3 \over 2} gmet(i,l) gmet(j,m) - {1 \over 2} gmet(i,j) gmet(l,m)] aa(l,m)$
(6 elements in, 6 elements out)
rank=3: bb(i,j,k)= $[{5 \over 2} g(i,l) g(j,m) g(k,n) - {3 \over 2} g(i,j) g(l,m) g(k,n)] aa(l,m,n)$
(10 elements in, 10 elements out)
In this rank 3 case, the second term is NOT symmetric in all
permutations of i,j,k, but the final tensor b(ijk) may be
symmetrized over all permutations because it will be
contracted with a completely symmetric tensor.
The compressed storage scheme is based on storing a symmetric 3x3 matrix as
(1 . .)
(6 2 .)
(5 4 3)
which leads to the following mappings for all ranks
where the compressed storage index is to the right of the arrow:
rank=0 1->1 (only a scalar)
rank=1 1->1 2->2 3->3 (standard vector, no compression)
rank=2 11->1 22->2 33->3 32->4 31->5 21->6
weights 1 1 1 2 2 2
rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
weights 1 3 3 6 3 3 1 3 3 1
PARENTS
nonlop_pl
CHILDREN
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 78 subroutine metcon(rank,gmet,aa,bb) 79 80 use defs_basis 81 use m_errors 82 83 !This section has been created automatically by the script Abilint (TD). 84 !Do not modify the following lines by hand. 85 #undef ABI_FUNC 86 #define ABI_FUNC 'metcon' 87 !End of the abilint section 88 89 implicit none 90 91 !Arguments ------------------------------------ 92 !scalars 93 integer,intent(in) :: rank 94 !arrays 95 real(dp),intent(in) :: aa(2,((rank+1)*(rank+2))/2),gmet(3,3) 96 real(dp),intent(out) :: bb(2,((rank+1)*(rank+2))/2) 97 98 !Local variables------------------------------- 99 !scalars 100 integer :: ii,jj 101 real(dp) :: scalar,tmpiii,tmpijk 102 character(len=500) :: message 103 !arrays 104 real(dp) :: vector(3) 105 106 ! ************************************************************************* 107 108 !This statement function defines the l=3 contraction in 109 !terms of the free indices of the contracted tensor (re and im) 110 ! coniii(ii,i1,i2,i3)=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+& 111 !& gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+& 112 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10) 113 ! conijk(ii,i1,i2,i3)=aa(ii,4)*& 114 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+& 115 !& gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+& 116 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+& 117 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+& 118 !& gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+& 119 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,3)) 120 ! con(ii,i1,i2,i3)=coniii(ii,i1,i2,i3)+conijk(ii,i1,i2,i3)+& 121 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+& 122 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+& 123 !& gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+& 124 !& (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+& 125 !& gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+& 126 !& gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+& 127 !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+& 128 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+& 129 !& gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+& 130 !& (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+& 131 !& gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+& 132 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+& 133 !& (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+& 134 !& gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+& 135 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+& 136 !& (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+& 137 !& gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+& 138 !& gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8) 139 140 !DEBUG 141 !write(std_out,*)' metcon : enter ' 142 !stop 143 !ENDDEBUG 144 if (rank==0) then 145 ! Simply copy scalar, re and im 146 bb(1,1)=aa(1,1) 147 bb(2,1)=aa(2,1) 148 149 else if (rank==1) then 150 ! Apply gmet to input vector, re and im 151 do ii=1,2 152 do jj=1,3 153 bb(ii,jj)=gmet(jj,1)*aa(ii,1)+gmet(jj,2)*aa(ii,2)+gmet(jj,3)*aa(ii,3) 154 end do 155 end do 156 157 158 else if (rank==2) then 159 ! Apply rank2 expression, re and im 160 do ii=1,2 161 ! Carry out g(c,d)*aa(c,d) contraction to get scalar 162 scalar=gmet(1,1)*aa(ii,1)+gmet(2,2)*aa(ii,2)+& 163 & gmet(3,3)*aa(ii,3)+2.0d0*(gmet(2,1)*aa(ii,6)+& 164 & gmet(3,1)*aa(ii,5)+gmet(3,2)*aa(ii,4) ) 165 ! Write out components of contraction 166 ! (1,1)->1 167 bb(ii,1)=0.5d0*(3.0d0*(gmet(1,1)*gmet(1,1)*aa(ii,1)+& 168 & gmet(1,2)*gmet(1,2)*aa(ii,2)+gmet(1,3)*gmet(1,3)*aa(ii,3)+& 169 & (gmet(1,2)*gmet(1,1)+gmet(1,1)*gmet(1,2))*aa(ii,6)+& 170 & (gmet(1,3)*gmet(1,1)+gmet(1,1)*gmet(1,3))*aa(ii,5)+& 171 & (gmet(1,3)*gmet(1,2)+gmet(1,2)*gmet(1,3))*aa(ii,4) ) & 172 & - gmet(1,1)*scalar) 173 ! (2,2)->2 174 bb(ii,2)=0.5d0*(3.0d0*(gmet(2,1)*gmet(2,1)*aa(ii,1)+& 175 & gmet(2,2)*gmet(2,2)*aa(ii,2)+gmet(2,3)*gmet(2,3)*aa(ii,3)+& 176 & (gmet(2,2)*gmet(2,1)+gmet(2,1)*gmet(2,2))*aa(ii,6)+& 177 & (gmet(2,3)*gmet(2,1)+gmet(2,1)*gmet(2,3))*aa(ii,5)+& 178 & (gmet(2,3)*gmet(2,2)+gmet(2,2)*gmet(2,3))*aa(ii,4) )& 179 & - gmet(2,2)*scalar) 180 ! (3,3)->3 181 bb(ii,3)=0.5d0*(3.0d0*(gmet(3,1)*gmet(3,1)*aa(ii,1)+& 182 & gmet(3,2)*gmet(3,2)*aa(ii,2)+gmet(3,3)*gmet(3,3)*aa(ii,3)+& 183 & (gmet(3,2)*gmet(3,1)+gmet(3,1)*gmet(3,2))*aa(ii,6)+& 184 & (gmet(3,3)*gmet(3,1)+gmet(3,1)*gmet(3,3))*aa(ii,5)+& 185 & (gmet(3,3)*gmet(3,2)+gmet(3,2)*gmet(3,3))*aa(ii,4) )& 186 & - gmet(3,3)*scalar) 187 ! (3,2)->4 188 bb(ii,4)=0.5d0*(3.0d0*(gmet(3,1)*gmet(2,1)*aa(ii,1)+& 189 & gmet(3,2)*gmet(2,2)*aa(ii,2)+gmet(3,3)*gmet(2,3)*aa(ii,3)+& 190 & (gmet(3,2)*gmet(2,1)+gmet(3,1)*gmet(2,2))*aa(ii,6)+& 191 & (gmet(3,3)*gmet(2,1)+gmet(3,1)*gmet(2,3))*aa(ii,5)+& 192 & (gmet(3,3)*gmet(2,2)+gmet(3,2)*gmet(2,3))*aa(ii,4) )& 193 & - gmet(3,2)*scalar) 194 ! (3,1)->5 195 bb(ii,5)=0.5d0*(3.0d0*(gmet(3,1)*gmet(1,1)*aa(ii,1)+& 196 & gmet(3,2)*gmet(1,2)*aa(ii,2)+gmet(3,3)*gmet(1,3)*aa(ii,3)+& 197 & (gmet(3,2)*gmet(1,1)+gmet(3,1)*gmet(1,2))*aa(ii,6)+& 198 & (gmet(3,3)*gmet(1,1)+gmet(3,1)*gmet(1,3))*aa(ii,5)+& 199 & (gmet(3,3)*gmet(1,2)+gmet(3,2)*gmet(1,3))*aa(ii,4) )& 200 & - gmet(3,1)*scalar) 201 ! (2,1)->6 202 bb(ii,6)=0.5d0*(3.0d0*(gmet(2,1)*gmet(1,1)*aa(ii,1)+& 203 & gmet(2,2)*gmet(1,2)*aa(ii,2)+gmet(2,3)*gmet(1,3)*aa(ii,3)+& 204 & (gmet(2,2)*gmet(1,1)+gmet(2,1)*gmet(1,2))*aa(ii,6)+& 205 & (gmet(2,3)*gmet(1,1)+gmet(2,1)*gmet(1,3))*aa(ii,5)+& 206 & (gmet(2,3)*gmet(1,2)+gmet(2,2)*gmet(1,3))*aa(ii,4) ) & 207 & - gmet(2,1)*scalar) 208 ! Include appropriate weights for multiplicity 209 bb(ii,4)=2.d0*bb(ii,4) 210 bb(ii,5)=2.d0*bb(ii,5) 211 bb(ii,6)=2.d0*bb(ii,6) 212 end do 213 214 else if (rank==3) then 215 ! Apply rank2 expression, re and im 216 do ii=1,2 217 ! Carry out g(l,m)g(j,n)*aa(l,m,n) contraction to get vector(j) 218 do jj=1,3 219 tmpiii= gmet(1,1)*gmet(jj,1)*aa(ii,1)+& 220 & gmet(2,2)*gmet(jj,2)*aa(ii,7)+& 221 & gmet(3,3)*gmet(jj,3)*aa(ii,10) 222 tmpijk= (gmet(1,2)*gmet(jj,3)+& 223 & gmet(3,1)*gmet(jj,2)+& 224 & gmet(2,3)*gmet(jj,1)+& 225 & gmet(3,2)*gmet(jj,1)+& 226 & gmet(1,3)*gmet(jj,2)+& 227 & gmet(2,1)*gmet(jj,3)) *aa(ii,4) 228 vector(jj)=tmpiii + tmpijk +& 229 & (gmet(1,2)*gmet(jj,1)+& 230 & gmet(2,1)*gmet(jj,1)+& 231 & gmet(1,1)*gmet(jj,2)) *aa(ii,6)+& 232 & (gmet(1,2)*gmet(jj,2)+& 233 & gmet(2,1)*gmet(jj,2)+& 234 & gmet(2,2)*gmet(jj,1)) *aa(ii,2)+& 235 & (gmet(1,3)*gmet(jj,1)+& 236 & gmet(3,1)*gmet(jj,1)+& 237 & gmet(1,1)*gmet(jj,3)) *aa(ii,5)+& 238 & (gmet(1,3)*gmet(jj,3)+& 239 & gmet(3,1)*gmet(jj,3)+& 240 & gmet(3,3)*gmet(jj,1)) *aa(ii,3)+& 241 & (gmet(2,3)*gmet(jj,2)+& 242 & gmet(3,2)*gmet(jj,2)+& 243 & gmet(2,2)*gmet(jj,3)) *aa(ii,9)+& 244 & (gmet(2,3)*gmet(jj,3)+& 245 & gmet(3,2)*gmet(jj,3)+& 246 & gmet(3,3)*gmet(jj,2)) *aa(ii,8) 247 end do 248 ! Write out components of contraction 249 ! (111)->1 250 bb(ii,1) =2.5d0*con_met(ii,1,1,1)-1.5d0*(gmet(1,1)*vector(1)) 251 ! (221)->2 252 bb(ii,2) =2.5d0*con_met(ii,2,2,1)-0.5d0*(gmet(1,2)*vector(2)+& 253 & gmet(1,2)*vector(2)+gmet(2,2)*vector(1)) 254 ! (331)->3 255 bb(ii,3) =2.5d0*con_met(ii,3,3,1)-0.5d0*(gmet(1,3)*vector(3)+& 256 & gmet(1,3)*vector(3)+gmet(3,3)*vector(1)) 257 ! (321)->4 258 bb(ii,4) =2.5d0*con_met(ii,3,2,1)-0.5d0*(gmet(1,3)*vector(2)+& 259 & gmet(1,2)*vector(3)+gmet(3,2)*vector(1)) 260 ! (311)->5 261 bb(ii,5) =2.5d0*con_met(ii,3,1,1)-0.5d0*(gmet(1,3)*vector(1)+& 262 & gmet(1,1)*vector(3)+gmet(3,1)*vector(1)) 263 ! (211)->6 264 bb(ii,6) =2.5d0*con_met(ii,2,1,1)-0.5d0*(gmet(1,2)*vector(1)+& 265 & gmet(1,1)*vector(2)+gmet(2,1)*vector(1)) 266 ! (222)->7 267 bb(ii,7) =2.5d0*con_met(ii,2,2,2)-1.5d0*(gmet(2,2)*vector(2)) 268 269 ! (332)->8 270 bb(ii,8) =2.5d0*con_met(ii,3,3,2)-0.5d0*(gmet(2,3)*vector(3)+& 271 & gmet(2,3)*vector(3)+gmet(3,3)*vector(2)) 272 ! (322)->9 273 bb(ii,9) =2.5d0*con_met(ii,3,2,2)-0.5d0*(gmet(2,3)*vector(2)+& 274 & gmet(2,2)*vector(3)+gmet(3,2)*vector(2)) 275 ! (333)->10 276 bb(ii,10)=2.5d0*con_met(ii,3,3,3)-1.5d0*(gmet(3,3)*vector(3)) 277 ! Include appropriate weights for multiplicity 278 bb(ii,2)=3.d0*bb(ii,2) 279 bb(ii,3)=3.d0*bb(ii,3) 280 bb(ii,4)=6.d0*bb(ii,4) 281 bb(ii,5)=3.d0*bb(ii,5) 282 bb(ii,6)=3.d0*bb(ii,6) 283 bb(ii,8)=3.d0*bb(ii,8) 284 bb(ii,9)=3.d0*bb(ii,9) 285 end do 286 287 else 288 write(message, '(a,i0,a,a,a)' )& 289 & 'Input rank=',rank,' not allowed.',ch10,& 290 & 'Possible values are 0,1,2,3 only.' 291 MSG_BUG(message) 292 end if 293 294 contains 295 296 function con_met(ii,i1,i2,i3) 297 298 299 !This section has been created automatically by the script Abilint (TD). 300 !Do not modify the following lines by hand. 301 #undef ABI_FUNC 302 #define ABI_FUNC 'con_met' 303 !End of the abilint section 304 305 real(dp) :: con_met 306 integer :: ii,i1,i2,i3 307 real(dp)::coniii,conijk 308 309 coniii=gmet(i1,1)*gmet(i2,1)*gmet(i3,1)*aa(ii,1)+& 310 & gmet(i1,2)*gmet(i2,2)*gmet(i3,2)*aa(ii,7)+& 311 & gmet(i1,3)*gmet(i2,3)*gmet(i3,3)*aa(ii,10) 312 conijk=aa(ii,4)*& 313 & (gmet(i1,1)*gmet(i2,2)*gmet(i3,3)+& 314 & gmet(i1,2)*gmet(i2,3)*gmet(i3,1)+& 315 & gmet(i1,3)*gmet(i2,1)*gmet(i3,2)+& 316 & gmet(i1,3)*gmet(i2,2)*gmet(i3,1)+& 317 & gmet(i1,1)*gmet(i2,3)*gmet(i3,2)+& 318 & gmet(i1,2)*gmet(i2,1)*gmet(i3,3)) 319 con_met=coniii+conijk+& 320 & (gmet(i1,1)*gmet(i2,2)*gmet(i3,1)+& 321 & gmet(i1,2)*gmet(i2,1)*gmet(i3,1)+& 322 & gmet(i1,1)*gmet(i2,1)*gmet(i3,2))*aa(ii,6)+& 323 & (gmet(i1,1)*gmet(i2,2)*gmet(i3,2)+& 324 & gmet(i1,2)*gmet(i2,1)*gmet(i3,2)+& 325 & gmet(i1,2)*gmet(i2,2)*gmet(i3,1))*aa(ii,2)+& 326 & (gmet(i1,1)*gmet(i2,3)*gmet(i3,1)+& 327 & gmet(i1,3)*gmet(i2,1)*gmet(i3,1)+& 328 & gmet(i1,1)*gmet(i2,1)*gmet(i3,3))*aa(ii,5)+& 329 & (gmet(i1,1)*gmet(i2,3)*gmet(i3,3)+& 330 & gmet(i1,3)*gmet(i2,1)*gmet(i3,3)+& 331 & gmet(i1,3)*gmet(i2,3)*gmet(i3,1))*aa(ii,3)+& 332 & (gmet(i1,2)*gmet(i2,2)*gmet(i3,3)+& 333 & gmet(i1,2)*gmet(i2,3)*gmet(i3,2)+& 334 & gmet(i1,3)*gmet(i2,2)*gmet(i3,2))*aa(ii,9)+& 335 & (gmet(i1,2)*gmet(i2,3)*gmet(i3,3)+& 336 & gmet(i1,3)*gmet(i2,2)*gmet(i3,3)+& 337 & gmet(i1,3)*gmet(i2,3)*gmet(i3,2))*aa(ii,8) 338 339 end function con_met 340 341 end subroutine metcon