TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
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.
PARENTS
opernl2
CHILDREN
SOURCE
60 #if defined HAVE_CONFIG_H 61 #include "config.h" 62 #endif 63 64 #include "abi_common.h" 65 66 subroutine mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 67 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 68 & npw,ntens,ntypat,parity) 69 70 use defs_basis 71 use m_profiling_abi 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'mkffkg' 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments ------------------------------------ 82 !scalars 83 integer,intent(in) :: choice,idir,ipw1,ispinor,itypat,lmnmax,mblkpw,ndgxdt 84 integer,intent(in) :: nffkg,nffnl,nincpw,nkpg,nlang,npw,ntens,ntypat 85 !arrays 86 integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw) 87 integer,intent(out) :: parity(nffkg) 88 real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg) 89 real(dp),intent(in) :: kpt(3) 90 real(dp),intent(out) :: ffkg(mblkpw,nffkg),kpgx(mblkpw,ntens) 91 92 !Local variables------------------------------- 93 !scalars 94 integer :: iffkg,ig,ii,ilang,ilang2,ilangx,ilmn,iln,iln0,iproj,ipw,jj 95 integer :: nffkge 96 real(dp) :: ffkg_now,kpg_x,kpg_y,kpg_z 97 !arrays 98 99 ! ************************************************************************* 100 101 jj=0;ilangx=0 102 103 !This will be useless after all the modifications have been done 104 do ipw=1,nincpw 105 kpgx(ipw,1)=1.0d0 106 end do 107 108 !Initialize kpgx array related to tensors defined below 109 if ( nlang>=2 .or. choice==2 .or. choice==3 .or. choice==4 .or. choice==5& 110 & .or. choice==6 .or. choice==23) then 111 if (nkpg>=3) then 112 kpgx(1:nincpw,2)=kpg_k(ipw1+1:ipw1+nincpw,1) 113 kpgx(1:nincpw,3)=kpg_k(ipw1+1:ipw1+nincpw,2) 114 kpgx(1:nincpw,4)=kpg_k(ipw1+1:ipw1+nincpw,3) 115 else 116 ig=ipw1 117 do ipw=1,nincpw 118 kpgx(ipw,2)=kpt(1)+dble(kg_k(1,ig)) 119 kpgx(ipw,3)=kpt(2)+dble(kg_k(2,ig)) 120 kpgx(ipw,4)=kpt(3)+dble(kg_k(3,ig)) 121 ig=ig+1 122 end do 123 end if 124 end if 125 if (nlang>=3 .or. choice==3 .or. choice==6 .or. choice==23) then 126 ! Define (k+G) part of rank 2 symmetric tensor (6 components), l=2 127 ! Compressed storage is 11 22 33 32 31 21 128 if (nkpg>=9) then 129 kpgx(1:nincpw,5) =kpg_k(ipw1+1:ipw1+nincpw,4) 130 kpgx(1:nincpw,6) =kpg_k(ipw1+1:ipw1+nincpw,5) 131 kpgx(1:nincpw,7) =kpg_k(ipw1+1:ipw1+nincpw,6) 132 kpgx(1:nincpw,8) =kpg_k(ipw1+1:ipw1+nincpw,7) 133 kpgx(1:nincpw,9) =kpg_k(ipw1+1:ipw1+nincpw,8) 134 kpgx(1:nincpw,10)=kpg_k(ipw1+1:ipw1+nincpw,9) 135 else 136 do ipw=1,nincpw 137 kpgx(ipw, 5) = kpgx(ipw, 2)*kpgx(ipw, 2) 138 kpgx(ipw, 6) = kpgx(ipw, 3)*kpgx(ipw, 3) 139 kpgx(ipw, 7) = kpgx(ipw, 4)*kpgx(ipw, 4) 140 kpgx(ipw, 8) = kpgx(ipw, 4)*kpgx(ipw, 3) 141 kpgx(ipw, 9) = kpgx(ipw, 4)*kpgx(ipw, 2) 142 kpgx(ipw,10) = kpgx(ipw, 3)*kpgx(ipw, 2) 143 end do 144 end if 145 end if 146 if (nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) .or. choice==6) then 147 ! Define (k+G) part of rank 3 symmetric tensor (10 components), l=3 148 ! Compressed storage is 111 221 331 321 311 211 222 332 322 333 149 do ipw=1,nincpw 150 kpgx(ipw,11) = kpgx(ipw, 5)*kpgx(ipw, 2) 151 kpgx(ipw,12) = kpgx(ipw, 6)*kpgx(ipw, 2) 152 kpgx(ipw,13) = kpgx(ipw, 7)*kpgx(ipw, 2) 153 kpgx(ipw,14) = kpgx(ipw, 8)*kpgx(ipw, 2) 154 kpgx(ipw,15) = kpgx(ipw, 9)*kpgx(ipw, 2) 155 kpgx(ipw,16) = kpgx(ipw,10)*kpgx(ipw, 2) 156 kpgx(ipw,17) = kpgx(ipw, 6)*kpgx(ipw, 3) 157 kpgx(ipw,18) = kpgx(ipw, 7)*kpgx(ipw, 3) 158 kpgx(ipw,19) = kpgx(ipw, 8)*kpgx(ipw, 3) 159 kpgx(ipw,20) = kpgx(ipw, 7)*kpgx(ipw, 4) 160 end do 161 end if 162 if (((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6) then 163 ! Add additional tensors for strain gradients 164 ! Define (k+G) part of rank 4 symmetric tensor (15 components), l=2 165 ! Compressed storage is 1111 2211 3311 3211 3111 2111 2221 3321 3221 166 ! 3331 2222 3322 3222 3332 3333 167 do ipw=1,nincpw 168 kpgx(ipw,21) = kpgx(ipw, 5)*kpgx(ipw, 5) 169 kpgx(ipw,22) = kpgx(ipw, 6)*kpgx(ipw, 5) 170 kpgx(ipw,23) = kpgx(ipw, 7)*kpgx(ipw, 5) 171 kpgx(ipw,24) = kpgx(ipw, 8)*kpgx(ipw, 5) 172 kpgx(ipw,25) = kpgx(ipw, 9)*kpgx(ipw, 5) 173 kpgx(ipw,26) = kpgx(ipw,10)*kpgx(ipw, 5) 174 kpgx(ipw,27) = kpgx(ipw, 6)*kpgx(ipw,10) 175 kpgx(ipw,28) = kpgx(ipw, 7)*kpgx(ipw,10) 176 kpgx(ipw,29) = kpgx(ipw, 8)*kpgx(ipw,10) 177 kpgx(ipw,30) = kpgx(ipw, 7)*kpgx(ipw, 9) 178 kpgx(ipw,31) = kpgx(ipw, 6)*kpgx(ipw, 6) 179 kpgx(ipw,32) = kpgx(ipw, 7)*kpgx(ipw, 6) 180 kpgx(ipw,33) = kpgx(ipw, 8)*kpgx(ipw, 6) 181 kpgx(ipw,34) = kpgx(ipw, 7)*kpgx(ipw, 8) 182 kpgx(ipw,35) = kpgx(ipw, 7)*kpgx(ipw, 7) 183 end do 184 end if 185 if (((choice==3.or.choice==23) .and. nlang>=4) .or. (choice==6 .and. nlang>=2)) then 186 ! Define (k+G) part of rank 5 symmetric tensor (21 components), l=3 187 ! Compressed storage is 11111 22111 33111 32111 31111 21111 188 ! 22211 33211 32211 33311 22221 33221 32221 33321 33331 189 ! 22222 33222 32222 33322 33332 33333 190 do ipw=1,nincpw 191 kpgx(ipw,36) = kpgx(ipw,21)*kpgx(ipw, 2) 192 kpgx(ipw,37) = kpgx(ipw,22)*kpgx(ipw, 2) 193 kpgx(ipw,38) = kpgx(ipw,23)*kpgx(ipw, 2) 194 kpgx(ipw,39) = kpgx(ipw,24)*kpgx(ipw, 2) 195 kpgx(ipw,40) = kpgx(ipw,25)*kpgx(ipw, 2) 196 kpgx(ipw,41) = kpgx(ipw,26)*kpgx(ipw, 2) 197 kpgx(ipw,42) = kpgx(ipw,27)*kpgx(ipw, 2) 198 kpgx(ipw,43) = kpgx(ipw,28)*kpgx(ipw, 2) 199 kpgx(ipw,44) = kpgx(ipw,29)*kpgx(ipw, 2) 200 kpgx(ipw,45) = kpgx(ipw,30)*kpgx(ipw, 2) 201 kpgx(ipw,46) = kpgx(ipw,31)*kpgx(ipw, 2) 202 kpgx(ipw,47) = kpgx(ipw,32)*kpgx(ipw, 2) 203 kpgx(ipw,48) = kpgx(ipw,33)*kpgx(ipw, 2) 204 kpgx(ipw,49) = kpgx(ipw,34)*kpgx(ipw, 2) 205 kpgx(ipw,50) = kpgx(ipw,35)*kpgx(ipw, 2) 206 kpgx(ipw,51) = kpgx(ipw,31)*kpgx(ipw, 3) 207 kpgx(ipw,52) = kpgx(ipw,32)*kpgx(ipw, 3) 208 kpgx(ipw,53) = kpgx(ipw,33)*kpgx(ipw, 3) 209 kpgx(ipw,54) = kpgx(ipw,34)*kpgx(ipw, 3) 210 kpgx(ipw,55) = kpgx(ipw,35)*kpgx(ipw, 3) 211 kpgx(ipw,56) = kpgx(ipw,35)*kpgx(ipw, 4) 212 end do 213 end if 214 if (choice==6 .and. nlang>=3) then 215 ! Define (k+G) part of rank 6 symmetric tensor (28 components) 216 ! Compressed storage is 217 ! 111111 221111 331111 321111 311111 211111 222111 332111 322111 218 ! 333111 222211 332211 322211 333211 333311 222221 332221 322221 219 ! 333221 333321 333331 222222 332222 322222 333222 333322 333332 220 ! 333333 221 do ipw=1,nincpw 222 kpgx(ipw,57) = kpgx(ipw,36)*kpgx(ipw, 2) 223 kpgx(ipw,58) = kpgx(ipw,37)*kpgx(ipw, 2) 224 kpgx(ipw,59) = kpgx(ipw,38)*kpgx(ipw, 2) 225 kpgx(ipw,60) = kpgx(ipw,39)*kpgx(ipw, 2) 226 kpgx(ipw,61) = kpgx(ipw,40)*kpgx(ipw, 2) 227 kpgx(ipw,62) = kpgx(ipw,41)*kpgx(ipw, 2) 228 kpgx(ipw,63) = kpgx(ipw,42)*kpgx(ipw, 2) 229 kpgx(ipw,64) = kpgx(ipw,43)*kpgx(ipw, 2) 230 kpgx(ipw,65) = kpgx(ipw,44)*kpgx(ipw, 2) 231 kpgx(ipw,66) = kpgx(ipw,45)*kpgx(ipw, 2) 232 kpgx(ipw,67) = kpgx(ipw,46)*kpgx(ipw, 2) 233 kpgx(ipw,68) = kpgx(ipw,47)*kpgx(ipw, 2) 234 kpgx(ipw,69) = kpgx(ipw,48)*kpgx(ipw, 2) 235 kpgx(ipw,70) = kpgx(ipw,49)*kpgx(ipw, 2) 236 kpgx(ipw,71) = kpgx(ipw,50)*kpgx(ipw, 2) 237 kpgx(ipw,72) = kpgx(ipw,51)*kpgx(ipw, 2) 238 kpgx(ipw,73) = kpgx(ipw,52)*kpgx(ipw, 2) 239 kpgx(ipw,74) = kpgx(ipw,53)*kpgx(ipw, 2) 240 kpgx(ipw,75) = kpgx(ipw,54)*kpgx(ipw, 2) 241 kpgx(ipw,76) = kpgx(ipw,55)*kpgx(ipw, 2) 242 kpgx(ipw,77) = kpgx(ipw,56)*kpgx(ipw, 2) 243 kpgx(ipw,78) = kpgx(ipw,51)*kpgx(ipw, 3) 244 kpgx(ipw,79) = kpgx(ipw,52)*kpgx(ipw, 3) 245 kpgx(ipw,80) = kpgx(ipw,53)*kpgx(ipw, 3) 246 kpgx(ipw,81) = kpgx(ipw,54)*kpgx(ipw, 3) 247 kpgx(ipw,82) = kpgx(ipw,55)*kpgx(ipw, 3) 248 kpgx(ipw,83) = kpgx(ipw,56)*kpgx(ipw, 3) 249 kpgx(ipw,84) = kpgx(ipw,56)*kpgx(ipw, 4) 250 end do 251 end if 252 if (choice==6 .and. nlang==4) then 253 ! Define (k+G) part of rank 7 symmetric tensor (36 components) 254 ! Compressed storage is 255 ! 1111111 2211111 3311111 3211111 3111111 2111111 2221111 3321111 3221111 256 ! 3331111 2222111 3322111 3222111 3332111 3333111 2222211 3322211 3222211 257 ! 3332211 3333211 3333311 2222221 3322221 3222221 3332221 3333221 3333321 258 ! 3333331 2222222 3322222 3222222 3332222 3333222 3333322 3333332 3333333 259 do ipw=1,nincpw 260 kpgx(ipw,85) = kpgx(ipw,57)*kpgx(ipw, 2) 261 kpgx(ipw,86) = kpgx(ipw,58)*kpgx(ipw, 2) 262 kpgx(ipw,87) = kpgx(ipw,59)*kpgx(ipw, 2) 263 kpgx(ipw,88) = kpgx(ipw,60)*kpgx(ipw, 2) 264 kpgx(ipw,89) = kpgx(ipw,61)*kpgx(ipw, 2) 265 kpgx(ipw,90) = kpgx(ipw,62)*kpgx(ipw, 2) 266 kpgx(ipw,91) = kpgx(ipw,63)*kpgx(ipw, 2) 267 kpgx(ipw,92) = kpgx(ipw,64)*kpgx(ipw, 2) 268 kpgx(ipw,93) = kpgx(ipw,65)*kpgx(ipw, 2) 269 kpgx(ipw,94) = kpgx(ipw,66)*kpgx(ipw, 2) 270 kpgx(ipw,95) = kpgx(ipw,67)*kpgx(ipw, 2) 271 kpgx(ipw,96) = kpgx(ipw,68)*kpgx(ipw, 2) 272 kpgx(ipw,97) = kpgx(ipw,69)*kpgx(ipw, 2) 273 kpgx(ipw,98) = kpgx(ipw,70)*kpgx(ipw, 2) 274 kpgx(ipw,99) = kpgx(ipw,71)*kpgx(ipw, 2) 275 kpgx(ipw,100) = kpgx(ipw,72)*kpgx(ipw, 2) 276 kpgx(ipw,101) = kpgx(ipw,73)*kpgx(ipw, 2) 277 kpgx(ipw,102) = kpgx(ipw,74)*kpgx(ipw, 2) 278 kpgx(ipw,103) = kpgx(ipw,75)*kpgx(ipw, 2) 279 kpgx(ipw,104) = kpgx(ipw,76)*kpgx(ipw, 2) 280 kpgx(ipw,105) = kpgx(ipw,77)*kpgx(ipw, 2) 281 kpgx(ipw,106) = kpgx(ipw,78)*kpgx(ipw, 2) 282 kpgx(ipw,107) = kpgx(ipw,79)*kpgx(ipw, 2) 283 kpgx(ipw,108) = kpgx(ipw,80)*kpgx(ipw, 2) 284 kpgx(ipw,109) = kpgx(ipw,81)*kpgx(ipw, 2) 285 kpgx(ipw,110) = kpgx(ipw,82)*kpgx(ipw, 2) 286 kpgx(ipw,111) = kpgx(ipw,83)*kpgx(ipw, 2) 287 kpgx(ipw,112) = kpgx(ipw,84)*kpgx(ipw, 2) 288 kpgx(ipw,113) = kpgx(ipw,78)*kpgx(ipw, 3) 289 kpgx(ipw,114) = kpgx(ipw,79)*kpgx(ipw, 3) 290 kpgx(ipw,115) = kpgx(ipw,80)*kpgx(ipw, 3) 291 kpgx(ipw,116) = kpgx(ipw,81)*kpgx(ipw, 3) 292 kpgx(ipw,117) = kpgx(ipw,82)*kpgx(ipw, 3) 293 kpgx(ipw,118) = kpgx(ipw,83)*kpgx(ipw, 3) 294 kpgx(ipw,119) = kpgx(ipw,84)*kpgx(ipw, 3) 295 kpgx(ipw,120) = kpgx(ipw,84)*kpgx(ipw, 4) 296 end do 297 end if 298 299 !***************************************************************************** 300 ! 301 !Packing of composite projectors in ffkg 302 303 iffkg=0 304 305 !Treat composite projectors for the energy 306 iln0=0 307 do ilmn=1,lmnmax 308 iln=indlmn(5,ilmn,itypat) 309 if (iln>iln0) then 310 iln0=iln 311 ilang=1+indlmn(1,ilmn,itypat) 312 iproj=indlmn(3,ilmn,itypat) 313 if(iproj>0)then 314 ilang2=(ilang*(ilang+1))/2 315 316 if(ilang==1)then 317 ! Treat s-component separately 318 ig=ipw1 319 iffkg=iffkg+1 320 do ipw=1,nincpw 321 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat) 322 ig=ig+1 323 end do 324 parity(iffkg)=2 325 else 326 ! Treat other components (could be made faster by treating explicitely 327 ! each angular momentum) 328 do ii=1,ilang2 329 ! Get the starting address for the relevant tensor 330 jj=ii+((ilang-1)*ilang*(ilang+1))/6 331 ig=ipw1 332 iffkg=iffkg+1 333 do ipw=1,nincpw 334 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 335 ig=ig+1 336 end do 337 if(ilang==2 .or. ilang==4)parity(iffkg)=1 338 if(ilang==3)parity(iffkg)=2 339 end do 340 end if 341 342 ! End condition if(iproj>0) 343 end if 344 ! End loop on ilang (ilmn) 345 end if 346 end do 347 348 !This is the number of composite projectors for the energy 349 nffkge=iffkg 350 351 !Second, treat forces : actually, this part could be rationalized, 352 !since the outcome is a multiplication by three of the number 353 !of composite projectors for the energy, while less should be needed 354 if((choice==2.or.choice==23) .and. ndgxdt/=1)then 355 do ii=1,nffkge 356 do ipw=1,nincpw 357 ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,2) 358 ffkg(ipw,iffkg+2)=ffkg(ipw,ii)*kpgx(ipw,3) 359 ffkg(ipw,iffkg+3)=ffkg(ipw,ii)*kpgx(ipw,4) 360 end do 361 parity(iffkg+1)=3-parity(ii) 362 parity(iffkg+2)=parity(iffkg+1) 363 parity(iffkg+3)=parity(iffkg+1) 364 iffkg=iffkg+3 365 end do 366 end if 367 !Note that the additional number of projectors for forces is 3*nffkge 368 369 !Third, treat first-derivative of the non-local operator 370 !with respect to an atomic displacement in one direction : 371 if(choice==2 .and. ndgxdt==1)then 372 do ii=1,nffkge 373 do ipw=1,nincpw 374 ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,idir+1) 375 end do 376 parity(iffkg+1)=3-parity(ii) 377 iffkg=iffkg+1 378 end do 379 end if 380 !Note that the additional number of projectors for this case is nffkge 381 382 383 !Fourth, treat dynamical matrices : like forces, this part could be rationalized. 384 if(choice==4)then 385 do ii=1,nffkge 386 do ipw=1,nincpw 387 kpg_x=kpgx(ipw,2) ; kpg_y=kpgx(ipw,3) ; kpg_z=kpgx(ipw,4) 388 ffkg_now=ffkg(ipw,ii) 389 ffkg(ipw,iffkg+1)=ffkg_now*kpg_x 390 ffkg(ipw,iffkg+2)=ffkg_now*kpg_y 391 ffkg(ipw,iffkg+3)=ffkg_now*kpg_z 392 ffkg(ipw,iffkg+4)=ffkg_now*kpg_x*kpg_x 393 ffkg(ipw,iffkg+5)=ffkg_now*kpg_y*kpg_y 394 ffkg(ipw,iffkg+6)=ffkg_now*kpg_z*kpg_z 395 ffkg(ipw,iffkg+7)=ffkg_now*kpg_z*kpg_y 396 ffkg(ipw,iffkg+8)=ffkg_now*kpg_z*kpg_x 397 ffkg(ipw,iffkg+9)=ffkg_now*kpg_y*kpg_x 398 end do 399 parity(iffkg+1:iffkg+3)=3-parity(ii) 400 parity(iffkg+4:iffkg+9)=parity(ii) 401 iffkg=iffkg+9 402 end do 403 end if 404 !Note that the additional number of projectors for dynamical matrices is 9*nffkge 405 406 !Treat composite projectors for the stress or 1st derivative contribution 407 !to frozen-wavefunction part of elastic tensor 408 !as well as, for ddk perturbation, the part that depend on ffnl(:,2,..) 409 if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then 410 411 iln0=0 412 do ilmn=1,lmnmax 413 if (ispinor==indlmn(6,ilmn,itypat)) then 414 iln=indlmn(5,ilmn,itypat) 415 if (iln>iln0) then 416 iln0=iln 417 ilang=1+indlmn(1,ilmn,itypat) 418 iproj=indlmn(3,ilmn,itypat) 419 if(iproj>0)then 420 ! number of unique tensor components 421 if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2 422 if(choice==5)ilangx=(ilang*(ilang+1))/2 423 424 do ii=1,ilangx 425 ! Get the starting address for the relevant tensor 426 if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 427 if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6 428 ig=ipw1 429 iffkg=iffkg+1 430 if(choice==3 .or. choice==6 .or. choice==23)then 431 do ipw=1,nincpw 432 ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 433 ig=ig+1 434 end do 435 else 436 do ipw=1,nincpw 437 ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)*& 438 & (kpgx(ipw,2)*gmet(1,idir)+ & 439 & kpgx(ipw,3)*gmet(2,idir)+ & 440 & kpgx(ipw,4)*gmet(3,idir) ) 441 ig=ig+1 442 end do 443 end if 444 if(ilang==1 .or. ilang==3)parity(iffkg)=2 445 if(ilang==2 .or. ilang==4)parity(iffkg)=1 446 if(choice==5)parity(iffkg)=3-parity(iffkg) 447 end do 448 449 ! End condition if(iproj>0) 450 end if 451 452 ! End loop on ilang (ilmn) 453 end if 454 end if 455 end do 456 457 ! End condition of stress 458 end if 459 460 !Treat composite projectors for the 2nd derivative wrt 2 strains 461 !and wrt one strain and one atomic displacement (internal strain) 462 !contributions to frozen-wavefunction part of (generalized) elastic tensor. 463 !There are 3 sets on terms (in historical order): 464 !first, terms with ffnl(:,3,...) and rank+4 tensors. 465 !second, terms with ffnl(:,1,...) and rank+1 tensors. 466 !third, terms with ffnl(:,2,...) and rank+3 tensors. 467 468 if(choice==6)then 469 470 iln0=0 471 do ilmn=1,lmnmax 472 if (ispinor==indlmn(6,ilmn,itypat)) then 473 iln=indlmn(5,ilmn,itypat) 474 if (iln>iln0) then 475 iln0=iln 476 ilang=1+indlmn(1,ilmn,itypat) 477 iproj=indlmn(3,ilmn,itypat) 478 if(iproj>0)then 479 ! First set of terms 480 ! number of unique tensor components 481 ilangx=((ilang+4)*(ilang+5))/2 482 483 do ii=1,ilangx 484 ! Get the starting address for the relevant tensor 485 jj=ii+((ilang+3)*(ilang+4)*(ilang+5))/6 486 ig=ipw1 487 iffkg=iffkg+1 488 do ipw=1,nincpw 489 ffkg(ipw,iffkg)=ffnl(ig,3,ilmn,itypat)*kpgx(ipw,jj) 490 ig=ig+1 491 end do 492 if(ilang==1 .or. ilang==3)parity(iffkg)=2 493 if(ilang==2 .or. ilang==4)parity(iffkg)=1 494 end do 495 496 ! Second set of terms 497 ! number of unique tensor components 498 ilangx=((ilang+1)*(ilang+2))/2 499 500 do ii=1,ilangx 501 ! Get the starting address for the relevant tensor 502 jj=ii+((ilang)*(ilang+1)*(ilang+2))/6 503 ig=ipw1 504 iffkg=iffkg+1 505 do ipw=1,nincpw 506 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 507 ig=ig+1 508 end do 509 if(ilang==1 .or. ilang==3)parity(iffkg)=1 510 if(ilang==2 .or. ilang==4)parity(iffkg)=2 511 end do 512 513 ! Third set of terms 514 ! number of unique tensor components 515 ilangx=((ilang+3)*(ilang+4))/2 516 517 do ii=1,ilangx 518 ! Get the starting address for the relevant tensor 519 jj=ii+((ilang+2)*(ilang+3)*(ilang+4))/6 520 ig=ipw1 521 iffkg=iffkg+1 522 do ipw=1,nincpw 523 ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 524 ig=ig+1 525 end do 526 if(ilang==1 .or. ilang==3)parity(iffkg)=1 527 if(ilang==2 .or. ilang==4)parity(iffkg)=2 528 end do 529 530 ! End condition if(iproj>0) 531 end if 532 ! End loop on ilang (ilmn) 533 end if 534 end if 535 end do 536 537 ! End condition of 2nd strain derivatives 538 end if 539 540 !For ddk perturbation, treat the part that depend on ffnl(:,1,..) 541 !no contribution from s state 542 if(nlang>=2 .and. choice==5)then 543 iln0=0 544 do ilmn=1,lmnmax 545 if (ispinor==indlmn(6,ilmn,itypat)) then 546 iln=indlmn(5,ilmn,itypat) 547 if (iln>iln0) then 548 iln0=iln 549 ilang=1+indlmn(1,ilmn,itypat) 550 if (ilang>=2) then 551 iproj=indlmn(3,ilmn,itypat) 552 if(iproj>0)then 553 ilang2=(ilang*(ilang-1))/2 554 555 do ii=1,ilang2 556 ! Get the starting address for the relevant tensor 557 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 558 ig=ipw1 559 iffkg=iffkg+1 560 do ipw=1,nincpw 561 ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 562 ig=ig+1 563 end do 564 if(ilang==2 .or. ilang==4)parity(iffkg)=2 565 if(ilang==3)parity(iffkg)=1 566 end do 567 568 ! End condition if(iproj>0) 569 end if 570 ! End loop on ilang>=2 571 end if 572 end if 573 end if 574 end do 575 ! End condition of p,d or f state 576 end if 577 578 end subroutine mkffkg