TABLE OF CONTENTS
ABINIT/dfpt_mkffkg [ Functions ]
NAME
dfpt_mkffkg
FUNCTION
Prepare the application of the projectors to the shifted wavefunctions, by precomputing the k+G factors and their product with the form factors Do this on a block of plane waves.
INPUTS
choice=governs the combination of k+G vectors to be computed ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere. gmet(3,3)=metric tensor for G vecs (in bohr**-2) nffnl=3rd dimension of ffnl(2, conventional, or 3 for 2nd derivatives) idir=direction of the perturbation (needed if choice==2 and ndgxdt==1, or if choice==5) indlmn(6,i,ntypat)=array giving l,m,n,lm,ln,spin for i=ln ipw1 = index of the first plane wave treated in this block ispinor=1 or 2, gives the spinorial component of ffnl to be used itypat = type of atom, needed for ffnl kg_k(3,npw)=integer coords of planewaves in basis sphere kpg_k(npw,npkg)= (k+G) components and related data kpt(3)=real components of k point in terms of recip. translations lmnmax=max. number of (l,n) components over all type of psps mblkpw=first dimension of kpgx ndgxdt=number of components of first order derivative nffkg=number of products of ffnls with combinations of k+G nincpw=number of plane waves in the block nkpg=second size of array kpg_k nlang = number of angular momenta to be treated = 1 + highest ang. mom. nloalg(3)=governs the choice of the algorithm for non-local operator. npw = total number of plane waves in reciprocal space ntens=second dimension of kpgx, number of distinct tensorial products ntypat = number of type of atoms, dimension needed for ffnl
OUTPUT
kpgx(mblkpw,ntens)=different tensorial products of k+G ffkg(nffkg,mblkpw)=different products of ffnls with k+G parity(nffkg)=parity of the tensorial product of k+G (2 if even, 1 of odd)
NOTES
This routine must be thread-safe as it is called inside loops that are OpenMP parallelized. Please, do not add variables with the save attribute or SIDE EFFECTS.
SOURCE
85 subroutine dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 86 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 87 & npw,ntens,ntypat,parity) 88 89 !Arguments ------------------------------------ 90 !scalars 91 integer,intent(in) :: choice,idir,ipw1,ispinor,itypat,lmnmax,mblkpw,ndgxdt 92 integer,intent(in) :: nffkg,nffnl,nincpw,nkpg,nlang,npw,ntens,ntypat 93 !arrays 94 integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw) 95 integer,intent(out) :: parity(nffkg) 96 real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg) 97 real(dp),intent(in) :: kpt(3) 98 real(dp),intent(out) :: ffkg(nffkg,mblkpw),kpgx(mblkpw,ntens) 99 100 !Local variables------------------------------- 101 !scalars 102 integer :: iffkg,ig,ii,ilang,ilang2,ilangx,ilmn,iln,iln0,iproj,ipw,jj 103 integer :: nffkge 104 real(dp) :: ffkg_now,kpg_x,kpg_y,kpg_z 105 106 ! ************************************************************************* 107 108 jj=0;ilangx=0 109 110 !This will be useless after all the modifications have been done 111 do ipw=1,nincpw 112 kpgx(ipw,1)=1.0d0 113 end do 114 115 !Initialize kpgx array related to tensors defined below 116 if ( nlang>=2 .or. choice==2 .or. choice==3 .or. choice==4 .or. choice==5& 117 & .or. choice==6 .or. choice==23) then 118 if (nkpg>=3) then 119 kpgx(1:nincpw,2)=kpg_k(ipw1+1:ipw1+nincpw,1) 120 kpgx(1:nincpw,3)=kpg_k(ipw1+1:ipw1+nincpw,2) 121 kpgx(1:nincpw,4)=kpg_k(ipw1+1:ipw1+nincpw,3) 122 else 123 ig=ipw1 124 do ipw=1,nincpw 125 kpgx(ipw,2)=kpt(1)+dble(kg_k(1,ig)) 126 kpgx(ipw,3)=kpt(2)+dble(kg_k(2,ig)) 127 kpgx(ipw,4)=kpt(3)+dble(kg_k(3,ig)) 128 ig=ig+1 129 end do 130 end if 131 end if 132 if (nlang>=3 .or. choice==3 .or. choice==6 .or. choice==23) then 133 ! Define (k+G) part of rank 2 symmetric tensor (6 components), l=2 134 ! Compressed storage is 11 22 33 32 31 21 135 if (nkpg>=9) then 136 kpgx(1:nincpw,5) =kpg_k(ipw1+1:ipw1+nincpw,4) 137 kpgx(1:nincpw,6) =kpg_k(ipw1+1:ipw1+nincpw,5) 138 kpgx(1:nincpw,7) =kpg_k(ipw1+1:ipw1+nincpw,6) 139 kpgx(1:nincpw,8) =kpg_k(ipw1+1:ipw1+nincpw,7) 140 kpgx(1:nincpw,9) =kpg_k(ipw1+1:ipw1+nincpw,8) 141 kpgx(1:nincpw,10)=kpg_k(ipw1+1:ipw1+nincpw,9) 142 else 143 do ipw=1,nincpw 144 kpgx(ipw, 5) = kpgx(ipw, 2)*kpgx(ipw, 2) 145 kpgx(ipw, 6) = kpgx(ipw, 3)*kpgx(ipw, 3) 146 kpgx(ipw, 7) = kpgx(ipw, 4)*kpgx(ipw, 4) 147 kpgx(ipw, 8) = kpgx(ipw, 4)*kpgx(ipw, 3) 148 kpgx(ipw, 9) = kpgx(ipw, 4)*kpgx(ipw, 2) 149 kpgx(ipw,10) = kpgx(ipw, 3)*kpgx(ipw, 2) 150 end do 151 end if 152 end if 153 if (nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) .or. choice==6) then 154 ! Define (k+G) part of rank 3 symmetric tensor (10 components), l=3 155 ! Compressed storage is 111 221 331 321 311 211 222 332 322 333 156 do ipw=1,nincpw 157 kpgx(ipw,11) = kpgx(ipw, 5)*kpgx(ipw, 2) 158 kpgx(ipw,12) = kpgx(ipw, 6)*kpgx(ipw, 2) 159 kpgx(ipw,13) = kpgx(ipw, 7)*kpgx(ipw, 2) 160 kpgx(ipw,14) = kpgx(ipw, 8)*kpgx(ipw, 2) 161 kpgx(ipw,15) = kpgx(ipw, 9)*kpgx(ipw, 2) 162 kpgx(ipw,16) = kpgx(ipw,10)*kpgx(ipw, 2) 163 kpgx(ipw,17) = kpgx(ipw, 6)*kpgx(ipw, 3) 164 kpgx(ipw,18) = kpgx(ipw, 7)*kpgx(ipw, 3) 165 kpgx(ipw,19) = kpgx(ipw, 8)*kpgx(ipw, 3) 166 kpgx(ipw,20) = kpgx(ipw, 7)*kpgx(ipw, 4) 167 end do 168 end if 169 if (((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6) then 170 ! Add additional tensors for strain gradients 171 ! Define (k+G) part of rank 4 symmetric tensor (15 components), l=2 172 ! Compressed storage is 1111 2211 3311 3211 3111 2111 2221 3321 3221 173 ! 3331 2222 3322 3222 3332 3333 174 do ipw=1,nincpw 175 kpgx(ipw,21) = kpgx(ipw, 5)*kpgx(ipw, 5) 176 kpgx(ipw,22) = kpgx(ipw, 6)*kpgx(ipw, 5) 177 kpgx(ipw,23) = kpgx(ipw, 7)*kpgx(ipw, 5) 178 kpgx(ipw,24) = kpgx(ipw, 8)*kpgx(ipw, 5) 179 kpgx(ipw,25) = kpgx(ipw, 9)*kpgx(ipw, 5) 180 kpgx(ipw,26) = kpgx(ipw,10)*kpgx(ipw, 5) 181 kpgx(ipw,27) = kpgx(ipw, 6)*kpgx(ipw,10) 182 kpgx(ipw,28) = kpgx(ipw, 7)*kpgx(ipw,10) 183 kpgx(ipw,29) = kpgx(ipw, 8)*kpgx(ipw,10) 184 kpgx(ipw,30) = kpgx(ipw, 7)*kpgx(ipw, 9) 185 kpgx(ipw,31) = kpgx(ipw, 6)*kpgx(ipw, 6) 186 kpgx(ipw,32) = kpgx(ipw, 7)*kpgx(ipw, 6) 187 kpgx(ipw,33) = kpgx(ipw, 8)*kpgx(ipw, 6) 188 kpgx(ipw,34) = kpgx(ipw, 7)*kpgx(ipw, 8) 189 kpgx(ipw,35) = kpgx(ipw, 7)*kpgx(ipw, 7) 190 end do 191 end if 192 if (((choice==3.or.choice==23) .and. nlang>=4) .or. (choice==6 .and. nlang>=2)) then 193 ! Define (k+G) part of rank 5 symmetric tensor (21 components), l=3 194 ! Compressed storage is 11111 22111 33111 32111 31111 21111 195 ! 22211 33211 32211 33311 22221 33221 32221 33321 33331 196 ! 22222 33222 32222 33322 33332 33333 197 do ipw=1,nincpw 198 kpgx(ipw,36) = kpgx(ipw,21)*kpgx(ipw, 2) 199 kpgx(ipw,37) = kpgx(ipw,22)*kpgx(ipw, 2) 200 kpgx(ipw,38) = kpgx(ipw,23)*kpgx(ipw, 2) 201 kpgx(ipw,39) = kpgx(ipw,24)*kpgx(ipw, 2) 202 kpgx(ipw,40) = kpgx(ipw,25)*kpgx(ipw, 2) 203 kpgx(ipw,41) = kpgx(ipw,26)*kpgx(ipw, 2) 204 kpgx(ipw,42) = kpgx(ipw,27)*kpgx(ipw, 2) 205 kpgx(ipw,43) = kpgx(ipw,28)*kpgx(ipw, 2) 206 kpgx(ipw,44) = kpgx(ipw,29)*kpgx(ipw, 2) 207 kpgx(ipw,45) = kpgx(ipw,30)*kpgx(ipw, 2) 208 kpgx(ipw,46) = kpgx(ipw,31)*kpgx(ipw, 2) 209 kpgx(ipw,47) = kpgx(ipw,32)*kpgx(ipw, 2) 210 kpgx(ipw,48) = kpgx(ipw,33)*kpgx(ipw, 2) 211 kpgx(ipw,49) = kpgx(ipw,34)*kpgx(ipw, 2) 212 kpgx(ipw,50) = kpgx(ipw,35)*kpgx(ipw, 2) 213 kpgx(ipw,51) = kpgx(ipw,31)*kpgx(ipw, 3) 214 kpgx(ipw,52) = kpgx(ipw,32)*kpgx(ipw, 3) 215 kpgx(ipw,53) = kpgx(ipw,33)*kpgx(ipw, 3) 216 kpgx(ipw,54) = kpgx(ipw,34)*kpgx(ipw, 3) 217 kpgx(ipw,55) = kpgx(ipw,35)*kpgx(ipw, 3) 218 kpgx(ipw,56) = kpgx(ipw,35)*kpgx(ipw, 4) 219 end do 220 end if 221 if (choice==6 .and. nlang>=3) then 222 ! Define (k+G) part of rank 6 symmetric tensor (28 components) 223 ! Compressed storage is 224 ! 111111 221111 331111 321111 311111 211111 222111 332111 322111 225 ! 333111 222211 332211 322211 333211 333311 222221 332221 322221 226 ! 333221 333321 333331 222222 332222 322222 333222 333322 333332 227 ! 333333 228 do ipw=1,nincpw 229 kpgx(ipw,57) = kpgx(ipw,36)*kpgx(ipw, 2) 230 kpgx(ipw,58) = kpgx(ipw,37)*kpgx(ipw, 2) 231 kpgx(ipw,59) = kpgx(ipw,38)*kpgx(ipw, 2) 232 kpgx(ipw,60) = kpgx(ipw,39)*kpgx(ipw, 2) 233 kpgx(ipw,61) = kpgx(ipw,40)*kpgx(ipw, 2) 234 kpgx(ipw,62) = kpgx(ipw,41)*kpgx(ipw, 2) 235 kpgx(ipw,63) = kpgx(ipw,42)*kpgx(ipw, 2) 236 kpgx(ipw,64) = kpgx(ipw,43)*kpgx(ipw, 2) 237 kpgx(ipw,65) = kpgx(ipw,44)*kpgx(ipw, 2) 238 kpgx(ipw,66) = kpgx(ipw,45)*kpgx(ipw, 2) 239 kpgx(ipw,67) = kpgx(ipw,46)*kpgx(ipw, 2) 240 kpgx(ipw,68) = kpgx(ipw,47)*kpgx(ipw, 2) 241 kpgx(ipw,69) = kpgx(ipw,48)*kpgx(ipw, 2) 242 kpgx(ipw,70) = kpgx(ipw,49)*kpgx(ipw, 2) 243 kpgx(ipw,71) = kpgx(ipw,50)*kpgx(ipw, 2) 244 kpgx(ipw,72) = kpgx(ipw,51)*kpgx(ipw, 2) 245 kpgx(ipw,73) = kpgx(ipw,52)*kpgx(ipw, 2) 246 kpgx(ipw,74) = kpgx(ipw,53)*kpgx(ipw, 2) 247 kpgx(ipw,75) = kpgx(ipw,54)*kpgx(ipw, 2) 248 kpgx(ipw,76) = kpgx(ipw,55)*kpgx(ipw, 2) 249 kpgx(ipw,77) = kpgx(ipw,56)*kpgx(ipw, 2) 250 kpgx(ipw,78) = kpgx(ipw,51)*kpgx(ipw, 3) 251 kpgx(ipw,79) = kpgx(ipw,52)*kpgx(ipw, 3) 252 kpgx(ipw,80) = kpgx(ipw,53)*kpgx(ipw, 3) 253 kpgx(ipw,81) = kpgx(ipw,54)*kpgx(ipw, 3) 254 kpgx(ipw,82) = kpgx(ipw,55)*kpgx(ipw, 3) 255 kpgx(ipw,83) = kpgx(ipw,56)*kpgx(ipw, 3) 256 kpgx(ipw,84) = kpgx(ipw,56)*kpgx(ipw, 4) 257 end do 258 end if 259 if (choice==6 .and. nlang==4) then 260 ! Define (k+G) part of rank 7 symmetric tensor (36 components) 261 ! Compressed storage is 262 ! 1111111 2211111 3311111 3211111 3111111 2111111 2221111 3321111 3221111 263 ! 3331111 2222111 3322111 3222111 3332111 3333111 2222211 3322211 3222211 264 ! 3332211 3333211 3333311 2222221 3322221 3222221 3332221 3333221 3333321 265 ! 3333331 2222222 3322222 3222222 3332222 3333222 3333322 3333332 3333333 266 do ipw=1,nincpw 267 kpgx(ipw,85) = kpgx(ipw,57)*kpgx(ipw, 2) 268 kpgx(ipw,86) = kpgx(ipw,58)*kpgx(ipw, 2) 269 kpgx(ipw,87) = kpgx(ipw,59)*kpgx(ipw, 2) 270 kpgx(ipw,88) = kpgx(ipw,60)*kpgx(ipw, 2) 271 kpgx(ipw,89) = kpgx(ipw,61)*kpgx(ipw, 2) 272 kpgx(ipw,90) = kpgx(ipw,62)*kpgx(ipw, 2) 273 kpgx(ipw,91) = kpgx(ipw,63)*kpgx(ipw, 2) 274 kpgx(ipw,92) = kpgx(ipw,64)*kpgx(ipw, 2) 275 kpgx(ipw,93) = kpgx(ipw,65)*kpgx(ipw, 2) 276 kpgx(ipw,94) = kpgx(ipw,66)*kpgx(ipw, 2) 277 kpgx(ipw,95) = kpgx(ipw,67)*kpgx(ipw, 2) 278 kpgx(ipw,96) = kpgx(ipw,68)*kpgx(ipw, 2) 279 kpgx(ipw,97) = kpgx(ipw,69)*kpgx(ipw, 2) 280 kpgx(ipw,98) = kpgx(ipw,70)*kpgx(ipw, 2) 281 kpgx(ipw,99) = kpgx(ipw,71)*kpgx(ipw, 2) 282 kpgx(ipw,100) = kpgx(ipw,72)*kpgx(ipw, 2) 283 kpgx(ipw,101) = kpgx(ipw,73)*kpgx(ipw, 2) 284 kpgx(ipw,102) = kpgx(ipw,74)*kpgx(ipw, 2) 285 kpgx(ipw,103) = kpgx(ipw,75)*kpgx(ipw, 2) 286 kpgx(ipw,104) = kpgx(ipw,76)*kpgx(ipw, 2) 287 kpgx(ipw,105) = kpgx(ipw,77)*kpgx(ipw, 2) 288 kpgx(ipw,106) = kpgx(ipw,78)*kpgx(ipw, 2) 289 kpgx(ipw,107) = kpgx(ipw,79)*kpgx(ipw, 2) 290 kpgx(ipw,108) = kpgx(ipw,80)*kpgx(ipw, 2) 291 kpgx(ipw,109) = kpgx(ipw,81)*kpgx(ipw, 2) 292 kpgx(ipw,110) = kpgx(ipw,82)*kpgx(ipw, 2) 293 kpgx(ipw,111) = kpgx(ipw,83)*kpgx(ipw, 2) 294 kpgx(ipw,112) = kpgx(ipw,84)*kpgx(ipw, 2) 295 kpgx(ipw,113) = kpgx(ipw,78)*kpgx(ipw, 3) 296 kpgx(ipw,114) = kpgx(ipw,79)*kpgx(ipw, 3) 297 kpgx(ipw,115) = kpgx(ipw,80)*kpgx(ipw, 3) 298 kpgx(ipw,116) = kpgx(ipw,81)*kpgx(ipw, 3) 299 kpgx(ipw,117) = kpgx(ipw,82)*kpgx(ipw, 3) 300 kpgx(ipw,118) = kpgx(ipw,83)*kpgx(ipw, 3) 301 kpgx(ipw,119) = kpgx(ipw,84)*kpgx(ipw, 3) 302 kpgx(ipw,120) = kpgx(ipw,84)*kpgx(ipw, 4) 303 end do 304 end if 305 306 !***************************************************************************** 307 ! 308 !Packing of composite projectors in ffkg 309 310 iffkg=0 311 312 !Treat composite projectors for the energy 313 iln0=0 314 do ilmn=1,lmnmax 315 iln=indlmn(5,ilmn,itypat) 316 if (iln>iln0) then 317 iln0=iln 318 ilang=1+indlmn(1,ilmn,itypat) 319 iproj=indlmn(3,ilmn,itypat) 320 if(iproj>0)then 321 ilang2=(ilang*(ilang+1))/2 322 323 if(ilang==1)then 324 ! Treat s-component separately 325 ig=ipw1 326 iffkg=iffkg+1 327 do ipw=1,nincpw 328 ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat) 329 ig=ig+1 330 end do 331 parity(iffkg)=2 332 else 333 ! Treat other components (could be made faster by treating explicitely 334 ! each angular momentum) 335 do ii=1,ilang2 336 ! Get the starting address for the relevant tensor 337 jj=ii+((ilang-1)*ilang*(ilang+1))/6 338 ig=ipw1 339 iffkg=iffkg+1 340 do ipw=1,nincpw 341 ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 342 ig=ig+1 343 end do 344 if(ilang==2 .or. ilang==4)parity(iffkg)=1 345 if(ilang==3)parity(iffkg)=2 346 end do 347 end if 348 349 ! End condition if(iproj>0) 350 end if 351 352 ! End loop on ilang (ilmn) 353 end if 354 end do 355 356 !This is the number of composite projectors for the energy 357 nffkge=iffkg 358 359 !Second, treat forces : actually, this part could be rationalized, 360 !since the outcome is a multiplication by three of the number 361 !of composite projectors for the energy, while less should be needed 362 if((choice==2.or.choice==23) .and. ndgxdt/=1)then 363 do ii=1,nffkge 364 do ipw=1,nincpw 365 ffkg(iffkg+1,ipw)=ffkg(ii,ipw)*kpgx(ipw,2) 366 ffkg(iffkg+2,ipw)=ffkg(ii,ipw)*kpgx(ipw,3) 367 ffkg(iffkg+3,ipw)=ffkg(ii,ipw)*kpgx(ipw,4) 368 end do 369 parity(iffkg+1)=3-parity(ii) 370 parity(iffkg+2)=parity(iffkg+1) 371 parity(iffkg+3)=parity(iffkg+1) 372 iffkg=iffkg+3 373 end do 374 end if 375 !Note that the additional number of projectors for forces is 3*nffkge 376 377 !Third, treat first-derivative of the non-local operator 378 !with respect to an atomic displacement in one direction : 379 if(choice==2 .and. ndgxdt==1)then 380 do ii=1,nffkge 381 do ipw=1,nincpw 382 ffkg(iffkg+1,ipw)=ffkg(ii,ipw)*kpgx(ipw,idir+1) 383 end do 384 parity(iffkg+1)=3-parity(ii) 385 iffkg=iffkg+1 386 end do 387 end if 388 !Note that the additional number of projectors for this case is nffkge 389 390 391 !Fourth, treat dynamical matrices : like forces, this part could be rationalized. 392 if(choice==4)then 393 do ii=1,nffkge 394 do ipw=1,nincpw 395 kpg_x=kpgx(ipw,2) ; kpg_y=kpgx(ipw,3) ; kpg_z=kpgx(ipw,4) 396 ffkg_now=ffkg(ii,ipw) 397 ffkg(iffkg+1,ipw)=ffkg_now*kpg_x 398 ffkg(iffkg+2,ipw)=ffkg_now*kpg_y 399 ffkg(iffkg+3,ipw)=ffkg_now*kpg_z 400 ffkg(iffkg+4,ipw)=ffkg_now*kpg_x*kpg_x 401 ffkg(iffkg+5,ipw)=ffkg_now*kpg_y*kpg_y 402 ffkg(iffkg+6,ipw)=ffkg_now*kpg_z*kpg_z 403 ffkg(iffkg+7,ipw)=ffkg_now*kpg_z*kpg_y 404 ffkg(iffkg+8,ipw)=ffkg_now*kpg_z*kpg_x 405 ffkg(iffkg+9,ipw)=ffkg_now*kpg_y*kpg_x 406 end do 407 parity(iffkg+1:iffkg+3)=3-parity(ii) 408 parity(iffkg+4:iffkg+9)=parity(ii) 409 iffkg=iffkg+9 410 end do 411 end if 412 !Note that the additional number of projectors for dynamical matrices is 9*nffkge 413 414 !Treat composite projectors for the stress or 1st derivative contribution 415 !to frozen-wavefunction part of elastic tensor 416 !as well as, for ddk perturbation, the part that depend on ffnl(:,2,..) 417 if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then 418 419 iln0=0 420 do ilmn=1,lmnmax 421 if (ispinor==indlmn(6,ilmn,itypat)) then 422 iln=indlmn(5,ilmn,itypat) 423 if (iln>iln0) then 424 iln0=iln 425 ilang=1+indlmn(1,ilmn,itypat) 426 iproj=indlmn(3,ilmn,itypat) 427 if(iproj>0)then 428 ! number of unique tensor components 429 if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2 430 if(choice==5)ilangx=(ilang*(ilang+1))/2 431 432 do ii=1,ilangx 433 ! Get the starting address for the relevant tensor 434 if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 435 if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6 436 ig=ipw1 437 iffkg=iffkg+1 438 if(choice==3 .or. choice==6 .or. choice==23)then 439 do ipw=1,nincpw 440 ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 441 ig=ig+1 442 end do 443 else 444 do ipw=1,nincpw 445 ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)*& 446 & (kpgx(ipw,2)*gmet(1,idir)+ & 447 & kpgx(ipw,3)*gmet(2,idir)+ & 448 & kpgx(ipw,4)*gmet(3,idir) ) 449 ig=ig+1 450 end do 451 end if 452 if(ilang==1 .or. ilang==3)parity(iffkg)=2 453 if(ilang==2 .or. ilang==4)parity(iffkg)=1 454 if(choice==5)parity(iffkg)=3-parity(iffkg) 455 end do 456 457 ! End condition if(iproj>0) 458 end if 459 ! End condition on iln 460 end if 461 ! End condition if(ispinor=indlmn(6,...)) 462 end if 463 ! End loop on ilmn 464 end do 465 466 ! End condition of stress 467 end if 468 469 !Treat composite projectors for the 2nd derivative wrt 2 strains 470 !and wrt one strain and one atomic displacement (internal strain) 471 !contributions to frozen-wavefunction part of (generalized) elastic tensor. 472 !There are 3 sets on terms (in historical order): 473 !first, terms with ffnl(:,3,...) and rank+4 tensors. 474 !second, terms with ffnl(:,1,...) and rank+1 tensors. 475 !third, terms with ffnl(:,2,...) and rank+3 tensors. 476 477 if(choice==6)then 478 479 iln0=0 480 do ilmn=1,lmnmax 481 if (ispinor==indlmn(6,ilmn,itypat)) then 482 iln=indlmn(5,ilmn,itypat) 483 if (iln>iln0) then 484 iln0=iln 485 ilang=1+indlmn(1,ilmn,itypat) 486 iproj=indlmn(3,ilmn,itypat) 487 488 if(iproj>0)then 489 ! First set of terms 490 ! number of unique tensor components 491 ilangx=((ilang+4)*(ilang+5))/2 492 493 do ii=1,ilangx 494 ! Get the starting address for the relevant tensor 495 jj=ii+((ilang+3)*(ilang+4)*(ilang+5))/6 496 ig=ipw1 497 iffkg=iffkg+1 498 do ipw=1,nincpw 499 ffkg(iffkg,ipw)=ffnl(ig,3,ilmn,itypat)*kpgx(ipw,jj) 500 ig=ig+1 501 end do 502 if(ilang==1 .or. ilang==3)parity(iffkg)=2 503 if(ilang==2 .or. ilang==4)parity(iffkg)=1 504 end do 505 506 ! Second set of terms 507 ! number of unique tensor components 508 ilangx=((ilang+1)*(ilang+2))/2 509 510 do ii=1,ilangx 511 ! Get the starting address for the relevant tensor 512 jj=ii+((ilang)*(ilang+1)*(ilang+2))/6 513 ig=ipw1 514 iffkg=iffkg+1 515 do ipw=1,nincpw 516 ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 517 ig=ig+1 518 end do 519 if(ilang==1 .or. ilang==3)parity(iffkg)=1 520 if(ilang==2 .or. ilang==4)parity(iffkg)=2 521 end do 522 523 ! Third set of terms 524 ! number of unique tensor components 525 ilangx=((ilang+3)*(ilang+4))/2 526 527 do ii=1,ilangx 528 ! Get the starting address for the relevant tensor 529 jj=ii+((ilang+2)*(ilang+3)*(ilang+4))/6 530 ig=ipw1 531 iffkg=iffkg+1 532 do ipw=1,nincpw 533 ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 534 ig=ig+1 535 end do 536 if(ilang==1 .or. ilang==3)parity(iffkg)=1 537 if(ilang==2 .or. ilang==4)parity(iffkg)=2 538 end do 539 540 ! End condition if(iproj>0) 541 end if 542 ! End condition on iln 543 end if 544 ! End condition if(ispinor=indlmn(6,...)) 545 end if 546 ! End loop on ilmn 547 end do 548 549 ! End condition of 2nd strain derivatives 550 end if 551 552 !For ddk perturbation, treat the part that depend on ffnl(:,1,..) 553 !no contribution from s state 554 if(nlang>=2 .and. choice==5)then 555 iln0=0 556 do ilmn=1,lmnmax 557 if (ispinor==indlmn(6,ilmn,itypat)) then 558 iln=indlmn(5,ilmn,itypat) 559 if (iln>iln0) then 560 iln0=iln 561 ilang=1+indlmn(1,ilmn,itypat) 562 if (ilang>=2) then 563 iproj=indlmn(3,ilmn,itypat) 564 if(iproj>0)then 565 ilang2=(ilang*(ilang-1))/2 566 567 do ii=1,ilang2 568 ! Get the starting address for the relevant tensor 569 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 570 ig=ipw1 571 iffkg=iffkg+1 572 do ipw=1,nincpw 573 ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 574 ig=ig+1 575 end do 576 if(ilang==2 .or. ilang==4)parity(iffkg)=2 577 if(ilang==3)parity(iffkg)=1 578 end do 579 580 ! End condition if(iproj>0) 581 end if 582 ! End condition if(ilang>=2) 583 end if 584 ! End condition if(iln>iln0) 585 end if 586 ! End condition if(ispinor=indlmn(6,...)) 587 end if 588 ! End loop on ilmn 589 end do 590 ! End condition of p,d or f state 591 end if 592 593 !DEBUG 594 !write(std_out,*)' dfpt_mkffkg : exit ' 595 !ENDDEBUG 596 597 end subroutine dfpt_mkffkg
ABINIT/m_mkffkg [ Modules ]
NAME
m_mkffkg
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, MT, DRH) 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_mkffkg 23 24 use defs_basis 25 use m_abicore 26 27 implicit none 28 29 private
ABINIT/mkffkg [ Functions ]
NAME
mkffkg
FUNCTION
Prepare the application of the projectors to the shifted wavefunctions, by precomputing the k+G factors and their product with the form factors Do this on a block of plane waves.
INPUTS
choice=governs the combination of k+G vectors to be computed ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere. gmet(3,3)=metric tensor for G vecs (in bohr**-2) nffnl=3rd dimension of ffnl(2, conventional, or 3 for 2nd derivatives) idir=direction of the perturbation (needed if choice==2 and ndgxdt==1, or if choice==5) indlmn(6,i,ntypat)=array giving l,m,n,lm,ln,spin for i=ln ipw1 = index of the first plane wave treated in this block ispinor=1 or 2, gives the spinorial component of ffnl to be used itypat = type of atom, needed for ffnl kg_k(3,npw)=integer coords of planewaves in basis sphere kpg_k(npw,npkg)= (k+G) components and related data kpt(3)=real components of k point in terms of recip. translations lmnmax=max. number of (l,n) components over all type of psps mblkpw=first dimension of kpgx ndgxdt=number of components of first order derivative nffkg=number of products of ffnls with combinations of k+G nincpw=number of plane waves in the block nkpg=second size of array kpg_k nlang = number of angular momenta to be treated = 1 + highest ang. mom. nloalg(3)=governs the choice of the algorithm for non-local operator. npw = total number of plane waves in reciprocal space ntens=second dimension of kpgx, number of distinct tensorial products ntypat = number of type of atoms, dimension needed for ffnl
OUTPUT
kpgx(mblkpw,ntens)=different tensorial products of k+G ffkg(mblkpw,nffkg)=different products of ffnls with k+G parity(nffkg)=parity of the tensorial product of k+G (2 if even, 1 of odd)
NOTES
This routine must be thread-safe as it is called inside loops that are OpenMP parallelized. Please, do not add variables with the save attribute or SIDE EFFECTS.
SOURCE
646 subroutine mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 647 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 648 & npw,ntens,ntypat,parity) 649 650 !Arguments ------------------------------------ 651 !scalars 652 integer,intent(in) :: choice,idir,ipw1,ispinor,itypat,lmnmax,mblkpw,ndgxdt 653 integer,intent(in) :: nffkg,nffnl,nincpw,nkpg,nlang,npw,ntens,ntypat 654 !arrays 655 integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw) 656 integer,intent(out) :: parity(nffkg) 657 real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg) 658 real(dp),intent(in) :: kpt(3) 659 real(dp),intent(out) :: ffkg(mblkpw,nffkg),kpgx(mblkpw,ntens) 660 661 !Local variables------------------------------- 662 !scalars 663 integer :: iffkg,ig,ii,ilang,ilang2,ilangx,ilmn,iln,iln0,iproj,ipw,jj 664 integer :: nffkge 665 real(dp) :: ffkg_now,kpg_x,kpg_y,kpg_z 666 !arrays 667 668 ! ************************************************************************* 669 670 jj=0;ilangx=0 671 672 !This will be useless after all the modifications have been done 673 do ipw=1,nincpw 674 kpgx(ipw,1)=1.0d0 675 end do 676 677 !Initialize kpgx array related to tensors defined below 678 if ( nlang>=2 .or. choice==2 .or. choice==3 .or. choice==4 .or. choice==5& 679 & .or. choice==6 .or. choice==23) then 680 if (nkpg>=3) then 681 kpgx(1:nincpw,2)=kpg_k(ipw1+1:ipw1+nincpw,1) 682 kpgx(1:nincpw,3)=kpg_k(ipw1+1:ipw1+nincpw,2) 683 kpgx(1:nincpw,4)=kpg_k(ipw1+1:ipw1+nincpw,3) 684 else 685 ig=ipw1 686 do ipw=1,nincpw 687 kpgx(ipw,2)=kpt(1)+dble(kg_k(1,ig)) 688 kpgx(ipw,3)=kpt(2)+dble(kg_k(2,ig)) 689 kpgx(ipw,4)=kpt(3)+dble(kg_k(3,ig)) 690 ig=ig+1 691 end do 692 end if 693 end if 694 if (nlang>=3 .or. choice==3 .or. choice==6 .or. choice==23) then 695 ! Define (k+G) part of rank 2 symmetric tensor (6 components), l=2 696 ! Compressed storage is 11 22 33 32 31 21 697 if (nkpg>=9) then 698 kpgx(1:nincpw,5) =kpg_k(ipw1+1:ipw1+nincpw,4) 699 kpgx(1:nincpw,6) =kpg_k(ipw1+1:ipw1+nincpw,5) 700 kpgx(1:nincpw,7) =kpg_k(ipw1+1:ipw1+nincpw,6) 701 kpgx(1:nincpw,8) =kpg_k(ipw1+1:ipw1+nincpw,7) 702 kpgx(1:nincpw,9) =kpg_k(ipw1+1:ipw1+nincpw,8) 703 kpgx(1:nincpw,10)=kpg_k(ipw1+1:ipw1+nincpw,9) 704 else 705 do ipw=1,nincpw 706 kpgx(ipw, 5) = kpgx(ipw, 2)*kpgx(ipw, 2) 707 kpgx(ipw, 6) = kpgx(ipw, 3)*kpgx(ipw, 3) 708 kpgx(ipw, 7) = kpgx(ipw, 4)*kpgx(ipw, 4) 709 kpgx(ipw, 8) = kpgx(ipw, 4)*kpgx(ipw, 3) 710 kpgx(ipw, 9) = kpgx(ipw, 4)*kpgx(ipw, 2) 711 kpgx(ipw,10) = kpgx(ipw, 3)*kpgx(ipw, 2) 712 end do 713 end if 714 end if 715 if (nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) .or. choice==6) then 716 ! Define (k+G) part of rank 3 symmetric tensor (10 components), l=3 717 ! Compressed storage is 111 221 331 321 311 211 222 332 322 333 718 do ipw=1,nincpw 719 kpgx(ipw,11) = kpgx(ipw, 5)*kpgx(ipw, 2) 720 kpgx(ipw,12) = kpgx(ipw, 6)*kpgx(ipw, 2) 721 kpgx(ipw,13) = kpgx(ipw, 7)*kpgx(ipw, 2) 722 kpgx(ipw,14) = kpgx(ipw, 8)*kpgx(ipw, 2) 723 kpgx(ipw,15) = kpgx(ipw, 9)*kpgx(ipw, 2) 724 kpgx(ipw,16) = kpgx(ipw,10)*kpgx(ipw, 2) 725 kpgx(ipw,17) = kpgx(ipw, 6)*kpgx(ipw, 3) 726 kpgx(ipw,18) = kpgx(ipw, 7)*kpgx(ipw, 3) 727 kpgx(ipw,19) = kpgx(ipw, 8)*kpgx(ipw, 3) 728 kpgx(ipw,20) = kpgx(ipw, 7)*kpgx(ipw, 4) 729 end do 730 end if 731 if (((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6) then 732 ! Add additional tensors for strain gradients 733 ! Define (k+G) part of rank 4 symmetric tensor (15 components), l=2 734 ! Compressed storage is 1111 2211 3311 3211 3111 2111 2221 3321 3221 735 ! 3331 2222 3322 3222 3332 3333 736 do ipw=1,nincpw 737 kpgx(ipw,21) = kpgx(ipw, 5)*kpgx(ipw, 5) 738 kpgx(ipw,22) = kpgx(ipw, 6)*kpgx(ipw, 5) 739 kpgx(ipw,23) = kpgx(ipw, 7)*kpgx(ipw, 5) 740 kpgx(ipw,24) = kpgx(ipw, 8)*kpgx(ipw, 5) 741 kpgx(ipw,25) = kpgx(ipw, 9)*kpgx(ipw, 5) 742 kpgx(ipw,26) = kpgx(ipw,10)*kpgx(ipw, 5) 743 kpgx(ipw,27) = kpgx(ipw, 6)*kpgx(ipw,10) 744 kpgx(ipw,28) = kpgx(ipw, 7)*kpgx(ipw,10) 745 kpgx(ipw,29) = kpgx(ipw, 8)*kpgx(ipw,10) 746 kpgx(ipw,30) = kpgx(ipw, 7)*kpgx(ipw, 9) 747 kpgx(ipw,31) = kpgx(ipw, 6)*kpgx(ipw, 6) 748 kpgx(ipw,32) = kpgx(ipw, 7)*kpgx(ipw, 6) 749 kpgx(ipw,33) = kpgx(ipw, 8)*kpgx(ipw, 6) 750 kpgx(ipw,34) = kpgx(ipw, 7)*kpgx(ipw, 8) 751 kpgx(ipw,35) = kpgx(ipw, 7)*kpgx(ipw, 7) 752 end do 753 end if 754 if (((choice==3.or.choice==23) .and. nlang>=4) .or. (choice==6 .and. nlang>=2)) then 755 ! Define (k+G) part of rank 5 symmetric tensor (21 components), l=3 756 ! Compressed storage is 11111 22111 33111 32111 31111 21111 757 ! 22211 33211 32211 33311 22221 33221 32221 33321 33331 758 ! 22222 33222 32222 33322 33332 33333 759 do ipw=1,nincpw 760 kpgx(ipw,36) = kpgx(ipw,21)*kpgx(ipw, 2) 761 kpgx(ipw,37) = kpgx(ipw,22)*kpgx(ipw, 2) 762 kpgx(ipw,38) = kpgx(ipw,23)*kpgx(ipw, 2) 763 kpgx(ipw,39) = kpgx(ipw,24)*kpgx(ipw, 2) 764 kpgx(ipw,40) = kpgx(ipw,25)*kpgx(ipw, 2) 765 kpgx(ipw,41) = kpgx(ipw,26)*kpgx(ipw, 2) 766 kpgx(ipw,42) = kpgx(ipw,27)*kpgx(ipw, 2) 767 kpgx(ipw,43) = kpgx(ipw,28)*kpgx(ipw, 2) 768 kpgx(ipw,44) = kpgx(ipw,29)*kpgx(ipw, 2) 769 kpgx(ipw,45) = kpgx(ipw,30)*kpgx(ipw, 2) 770 kpgx(ipw,46) = kpgx(ipw,31)*kpgx(ipw, 2) 771 kpgx(ipw,47) = kpgx(ipw,32)*kpgx(ipw, 2) 772 kpgx(ipw,48) = kpgx(ipw,33)*kpgx(ipw, 2) 773 kpgx(ipw,49) = kpgx(ipw,34)*kpgx(ipw, 2) 774 kpgx(ipw,50) = kpgx(ipw,35)*kpgx(ipw, 2) 775 kpgx(ipw,51) = kpgx(ipw,31)*kpgx(ipw, 3) 776 kpgx(ipw,52) = kpgx(ipw,32)*kpgx(ipw, 3) 777 kpgx(ipw,53) = kpgx(ipw,33)*kpgx(ipw, 3) 778 kpgx(ipw,54) = kpgx(ipw,34)*kpgx(ipw, 3) 779 kpgx(ipw,55) = kpgx(ipw,35)*kpgx(ipw, 3) 780 kpgx(ipw,56) = kpgx(ipw,35)*kpgx(ipw, 4) 781 end do 782 end if 783 if (choice==6 .and. nlang>=3) then 784 ! Define (k+G) part of rank 6 symmetric tensor (28 components) 785 ! Compressed storage is 786 ! 111111 221111 331111 321111 311111 211111 222111 332111 322111 787 ! 333111 222211 332211 322211 333211 333311 222221 332221 322221 788 ! 333221 333321 333331 222222 332222 322222 333222 333322 333332 789 ! 333333 790 do ipw=1,nincpw 791 kpgx(ipw,57) = kpgx(ipw,36)*kpgx(ipw, 2) 792 kpgx(ipw,58) = kpgx(ipw,37)*kpgx(ipw, 2) 793 kpgx(ipw,59) = kpgx(ipw,38)*kpgx(ipw, 2) 794 kpgx(ipw,60) = kpgx(ipw,39)*kpgx(ipw, 2) 795 kpgx(ipw,61) = kpgx(ipw,40)*kpgx(ipw, 2) 796 kpgx(ipw,62) = kpgx(ipw,41)*kpgx(ipw, 2) 797 kpgx(ipw,63) = kpgx(ipw,42)*kpgx(ipw, 2) 798 kpgx(ipw,64) = kpgx(ipw,43)*kpgx(ipw, 2) 799 kpgx(ipw,65) = kpgx(ipw,44)*kpgx(ipw, 2) 800 kpgx(ipw,66) = kpgx(ipw,45)*kpgx(ipw, 2) 801 kpgx(ipw,67) = kpgx(ipw,46)*kpgx(ipw, 2) 802 kpgx(ipw,68) = kpgx(ipw,47)*kpgx(ipw, 2) 803 kpgx(ipw,69) = kpgx(ipw,48)*kpgx(ipw, 2) 804 kpgx(ipw,70) = kpgx(ipw,49)*kpgx(ipw, 2) 805 kpgx(ipw,71) = kpgx(ipw,50)*kpgx(ipw, 2) 806 kpgx(ipw,72) = kpgx(ipw,51)*kpgx(ipw, 2) 807 kpgx(ipw,73) = kpgx(ipw,52)*kpgx(ipw, 2) 808 kpgx(ipw,74) = kpgx(ipw,53)*kpgx(ipw, 2) 809 kpgx(ipw,75) = kpgx(ipw,54)*kpgx(ipw, 2) 810 kpgx(ipw,76) = kpgx(ipw,55)*kpgx(ipw, 2) 811 kpgx(ipw,77) = kpgx(ipw,56)*kpgx(ipw, 2) 812 kpgx(ipw,78) = kpgx(ipw,51)*kpgx(ipw, 3) 813 kpgx(ipw,79) = kpgx(ipw,52)*kpgx(ipw, 3) 814 kpgx(ipw,80) = kpgx(ipw,53)*kpgx(ipw, 3) 815 kpgx(ipw,81) = kpgx(ipw,54)*kpgx(ipw, 3) 816 kpgx(ipw,82) = kpgx(ipw,55)*kpgx(ipw, 3) 817 kpgx(ipw,83) = kpgx(ipw,56)*kpgx(ipw, 3) 818 kpgx(ipw,84) = kpgx(ipw,56)*kpgx(ipw, 4) 819 end do 820 end if 821 if (choice==6 .and. nlang==4) then 822 ! Define (k+G) part of rank 7 symmetric tensor (36 components) 823 ! Compressed storage is 824 ! 1111111 2211111 3311111 3211111 3111111 2111111 2221111 3321111 3221111 825 ! 3331111 2222111 3322111 3222111 3332111 3333111 2222211 3322211 3222211 826 ! 3332211 3333211 3333311 2222221 3322221 3222221 3332221 3333221 3333321 827 ! 3333331 2222222 3322222 3222222 3332222 3333222 3333322 3333332 3333333 828 do ipw=1,nincpw 829 kpgx(ipw,85) = kpgx(ipw,57)*kpgx(ipw, 2) 830 kpgx(ipw,86) = kpgx(ipw,58)*kpgx(ipw, 2) 831 kpgx(ipw,87) = kpgx(ipw,59)*kpgx(ipw, 2) 832 kpgx(ipw,88) = kpgx(ipw,60)*kpgx(ipw, 2) 833 kpgx(ipw,89) = kpgx(ipw,61)*kpgx(ipw, 2) 834 kpgx(ipw,90) = kpgx(ipw,62)*kpgx(ipw, 2) 835 kpgx(ipw,91) = kpgx(ipw,63)*kpgx(ipw, 2) 836 kpgx(ipw,92) = kpgx(ipw,64)*kpgx(ipw, 2) 837 kpgx(ipw,93) = kpgx(ipw,65)*kpgx(ipw, 2) 838 kpgx(ipw,94) = kpgx(ipw,66)*kpgx(ipw, 2) 839 kpgx(ipw,95) = kpgx(ipw,67)*kpgx(ipw, 2) 840 kpgx(ipw,96) = kpgx(ipw,68)*kpgx(ipw, 2) 841 kpgx(ipw,97) = kpgx(ipw,69)*kpgx(ipw, 2) 842 kpgx(ipw,98) = kpgx(ipw,70)*kpgx(ipw, 2) 843 kpgx(ipw,99) = kpgx(ipw,71)*kpgx(ipw, 2) 844 kpgx(ipw,100) = kpgx(ipw,72)*kpgx(ipw, 2) 845 kpgx(ipw,101) = kpgx(ipw,73)*kpgx(ipw, 2) 846 kpgx(ipw,102) = kpgx(ipw,74)*kpgx(ipw, 2) 847 kpgx(ipw,103) = kpgx(ipw,75)*kpgx(ipw, 2) 848 kpgx(ipw,104) = kpgx(ipw,76)*kpgx(ipw, 2) 849 kpgx(ipw,105) = kpgx(ipw,77)*kpgx(ipw, 2) 850 kpgx(ipw,106) = kpgx(ipw,78)*kpgx(ipw, 2) 851 kpgx(ipw,107) = kpgx(ipw,79)*kpgx(ipw, 2) 852 kpgx(ipw,108) = kpgx(ipw,80)*kpgx(ipw, 2) 853 kpgx(ipw,109) = kpgx(ipw,81)*kpgx(ipw, 2) 854 kpgx(ipw,110) = kpgx(ipw,82)*kpgx(ipw, 2) 855 kpgx(ipw,111) = kpgx(ipw,83)*kpgx(ipw, 2) 856 kpgx(ipw,112) = kpgx(ipw,84)*kpgx(ipw, 2) 857 kpgx(ipw,113) = kpgx(ipw,78)*kpgx(ipw, 3) 858 kpgx(ipw,114) = kpgx(ipw,79)*kpgx(ipw, 3) 859 kpgx(ipw,115) = kpgx(ipw,80)*kpgx(ipw, 3) 860 kpgx(ipw,116) = kpgx(ipw,81)*kpgx(ipw, 3) 861 kpgx(ipw,117) = kpgx(ipw,82)*kpgx(ipw, 3) 862 kpgx(ipw,118) = kpgx(ipw,83)*kpgx(ipw, 3) 863 kpgx(ipw,119) = kpgx(ipw,84)*kpgx(ipw, 3) 864 kpgx(ipw,120) = kpgx(ipw,84)*kpgx(ipw, 4) 865 end do 866 end if 867 868 !***************************************************************************** 869 ! 870 !Packing of composite projectors in ffkg 871 872 iffkg=0 873 874 !Treat composite projectors for the energy 875 iln0=0 876 do ilmn=1,lmnmax 877 iln=indlmn(5,ilmn,itypat) 878 if (iln>iln0) then 879 iln0=iln 880 ilang=1+indlmn(1,ilmn,itypat) 881 iproj=indlmn(3,ilmn,itypat) 882 if(iproj>0)then 883 ilang2=(ilang*(ilang+1))/2 884 885 if(ilang==1)then 886 ! Treat s-component separately 887 ig=ipw1 888 iffkg=iffkg+1 889 do ipw=1,nincpw 890 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat) 891 ig=ig+1 892 end do 893 parity(iffkg)=2 894 else 895 ! Treat other components (could be made faster by treating explicitely 896 ! each angular momentum) 897 do ii=1,ilang2 898 ! Get the starting address for the relevant tensor 899 jj=ii+((ilang-1)*ilang*(ilang+1))/6 900 ig=ipw1 901 iffkg=iffkg+1 902 do ipw=1,nincpw 903 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 904 ig=ig+1 905 end do 906 if(ilang==2 .or. ilang==4)parity(iffkg)=1 907 if(ilang==3)parity(iffkg)=2 908 end do 909 end if 910 911 ! End condition if(iproj>0) 912 end if 913 ! End loop on ilang (ilmn) 914 end if 915 end do 916 917 !This is the number of composite projectors for the energy 918 nffkge=iffkg 919 920 !Second, treat forces : actually, this part could be rationalized, 921 !since the outcome is a multiplication by three of the number 922 !of composite projectors for the energy, while less should be needed 923 if((choice==2.or.choice==23) .and. ndgxdt/=1)then 924 do ii=1,nffkge 925 do ipw=1,nincpw 926 ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,2) 927 ffkg(ipw,iffkg+2)=ffkg(ipw,ii)*kpgx(ipw,3) 928 ffkg(ipw,iffkg+3)=ffkg(ipw,ii)*kpgx(ipw,4) 929 end do 930 parity(iffkg+1)=3-parity(ii) 931 parity(iffkg+2)=parity(iffkg+1) 932 parity(iffkg+3)=parity(iffkg+1) 933 iffkg=iffkg+3 934 end do 935 end if 936 !Note that the additional number of projectors for forces is 3*nffkge 937 938 !Third, treat first-derivative of the non-local operator 939 !with respect to an atomic displacement in one direction : 940 if(choice==2 .and. ndgxdt==1)then 941 do ii=1,nffkge 942 do ipw=1,nincpw 943 ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,idir+1) 944 end do 945 parity(iffkg+1)=3-parity(ii) 946 iffkg=iffkg+1 947 end do 948 end if 949 !Note that the additional number of projectors for this case is nffkge 950 951 952 !Fourth, treat dynamical matrices : like forces, this part could be rationalized. 953 if(choice==4)then 954 do ii=1,nffkge 955 do ipw=1,nincpw 956 kpg_x=kpgx(ipw,2) ; kpg_y=kpgx(ipw,3) ; kpg_z=kpgx(ipw,4) 957 ffkg_now=ffkg(ipw,ii) 958 ffkg(ipw,iffkg+1)=ffkg_now*kpg_x 959 ffkg(ipw,iffkg+2)=ffkg_now*kpg_y 960 ffkg(ipw,iffkg+3)=ffkg_now*kpg_z 961 ffkg(ipw,iffkg+4)=ffkg_now*kpg_x*kpg_x 962 ffkg(ipw,iffkg+5)=ffkg_now*kpg_y*kpg_y 963 ffkg(ipw,iffkg+6)=ffkg_now*kpg_z*kpg_z 964 ffkg(ipw,iffkg+7)=ffkg_now*kpg_z*kpg_y 965 ffkg(ipw,iffkg+8)=ffkg_now*kpg_z*kpg_x 966 ffkg(ipw,iffkg+9)=ffkg_now*kpg_y*kpg_x 967 end do 968 parity(iffkg+1:iffkg+3)=3-parity(ii) 969 parity(iffkg+4:iffkg+9)=parity(ii) 970 iffkg=iffkg+9 971 end do 972 end if 973 !Note that the additional number of projectors for dynamical matrices is 9*nffkge 974 975 !Treat composite projectors for the stress or 1st derivative contribution 976 !to frozen-wavefunction part of elastic tensor 977 !as well as, for ddk perturbation, the part that depend on ffnl(:,2,..) 978 if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then 979 980 iln0=0 981 do ilmn=1,lmnmax 982 if (ispinor==indlmn(6,ilmn,itypat)) then 983 iln=indlmn(5,ilmn,itypat) 984 if (iln>iln0) then 985 iln0=iln 986 ilang=1+indlmn(1,ilmn,itypat) 987 iproj=indlmn(3,ilmn,itypat) 988 if(iproj>0)then 989 ! number of unique tensor components 990 if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2 991 if(choice==5)ilangx=(ilang*(ilang+1))/2 992 993 do ii=1,ilangx 994 ! Get the starting address for the relevant tensor 995 if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 996 if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6 997 ig=ipw1 998 iffkg=iffkg+1 999 if(choice==3 .or. choice==6 .or. choice==23)then 1000 do ipw=1,nincpw 1001 ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 1002 ig=ig+1 1003 end do 1004 else 1005 do ipw=1,nincpw 1006 ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)*& 1007 & (kpgx(ipw,2)*gmet(1,idir)+ & 1008 & kpgx(ipw,3)*gmet(2,idir)+ & 1009 & kpgx(ipw,4)*gmet(3,idir) ) 1010 ig=ig+1 1011 end do 1012 end if 1013 if(ilang==1 .or. ilang==3)parity(iffkg)=2 1014 if(ilang==2 .or. ilang==4)parity(iffkg)=1 1015 if(choice==5)parity(iffkg)=3-parity(iffkg) 1016 end do 1017 1018 ! End condition if(iproj>0) 1019 end if 1020 1021 ! End loop on ilang (ilmn) 1022 end if 1023 end if 1024 end do 1025 1026 ! End condition of stress 1027 end if 1028 1029 !Treat composite projectors for the 2nd derivative wrt 2 strains 1030 !and wrt one strain and one atomic displacement (internal strain) 1031 !contributions to frozen-wavefunction part of (generalized) elastic tensor. 1032 !There are 3 sets on terms (in historical order): 1033 !first, terms with ffnl(:,3,...) and rank+4 tensors. 1034 !second, terms with ffnl(:,1,...) and rank+1 tensors. 1035 !third, terms with ffnl(:,2,...) and rank+3 tensors. 1036 1037 if(choice==6)then 1038 1039 iln0=0 1040 do ilmn=1,lmnmax 1041 if (ispinor==indlmn(6,ilmn,itypat)) then 1042 iln=indlmn(5,ilmn,itypat) 1043 if (iln>iln0) then 1044 iln0=iln 1045 ilang=1+indlmn(1,ilmn,itypat) 1046 iproj=indlmn(3,ilmn,itypat) 1047 if(iproj>0)then 1048 ! First set of terms 1049 ! number of unique tensor components 1050 ilangx=((ilang+4)*(ilang+5))/2 1051 1052 do ii=1,ilangx 1053 ! Get the starting address for the relevant tensor 1054 jj=ii+((ilang+3)*(ilang+4)*(ilang+5))/6 1055 ig=ipw1 1056 iffkg=iffkg+1 1057 do ipw=1,nincpw 1058 ffkg(ipw,iffkg)=ffnl(ig,3,ilmn,itypat)*kpgx(ipw,jj) 1059 ig=ig+1 1060 end do 1061 if(ilang==1 .or. ilang==3)parity(iffkg)=2 1062 if(ilang==2 .or. ilang==4)parity(iffkg)=1 1063 end do 1064 1065 ! Second set of terms 1066 ! number of unique tensor components 1067 ilangx=((ilang+1)*(ilang+2))/2 1068 1069 do ii=1,ilangx 1070 ! Get the starting address for the relevant tensor 1071 jj=ii+((ilang)*(ilang+1)*(ilang+2))/6 1072 ig=ipw1 1073 iffkg=iffkg+1 1074 do ipw=1,nincpw 1075 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 1076 ig=ig+1 1077 end do 1078 if(ilang==1 .or. ilang==3)parity(iffkg)=1 1079 if(ilang==2 .or. ilang==4)parity(iffkg)=2 1080 end do 1081 1082 ! Third set of terms 1083 ! number of unique tensor components 1084 ilangx=((ilang+3)*(ilang+4))/2 1085 1086 do ii=1,ilangx 1087 ! Get the starting address for the relevant tensor 1088 jj=ii+((ilang+2)*(ilang+3)*(ilang+4))/6 1089 ig=ipw1 1090 iffkg=iffkg+1 1091 do ipw=1,nincpw 1092 ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 1093 ig=ig+1 1094 end do 1095 if(ilang==1 .or. ilang==3)parity(iffkg)=1 1096 if(ilang==2 .or. ilang==4)parity(iffkg)=2 1097 end do 1098 1099 ! End condition if(iproj>0) 1100 end if 1101 ! End loop on ilang (ilmn) 1102 end if 1103 end if 1104 end do 1105 1106 ! End condition of 2nd strain derivatives 1107 end if 1108 1109 !For ddk perturbation, treat the part that depend on ffnl(:,1,..) 1110 !no contribution from s state 1111 if(nlang>=2 .and. choice==5)then 1112 iln0=0 1113 do ilmn=1,lmnmax 1114 if (ispinor==indlmn(6,ilmn,itypat)) then 1115 iln=indlmn(5,ilmn,itypat) 1116 if (iln>iln0) then 1117 iln0=iln 1118 ilang=1+indlmn(1,ilmn,itypat) 1119 if (ilang>=2) then 1120 iproj=indlmn(3,ilmn,itypat) 1121 if(iproj>0)then 1122 ilang2=(ilang*(ilang-1))/2 1123 1124 do ii=1,ilang2 1125 ! Get the starting address for the relevant tensor 1126 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 1127 ig=ipw1 1128 iffkg=iffkg+1 1129 do ipw=1,nincpw 1130 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 1131 ig=ig+1 1132 end do 1133 if(ilang==2 .or. ilang==4)parity(iffkg)=2 1134 if(ilang==3)parity(iffkg)=1 1135 end do 1136 1137 ! End condition if(iproj>0) 1138 end if 1139 ! End loop on ilang>=2 1140 end if 1141 end if 1142 end if 1143 end do 1144 ! End condition of p,d or f state 1145 end if 1146 1147 end subroutine mkffkg