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.
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(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.
PARENTS
opernl3,opernl4a,opernl4b
CHILDREN
SOURCE
60 #if defined HAVE_CONFIG_H 61 #include "config.h" 62 #endif 63 64 #include "abi_common.h" 65 66 67 subroutine dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 68 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 69 & npw,ntens,ntypat,parity) 70 71 use defs_basis 72 use m_profiling_abi 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'dfpt_mkffkg' 78 !End of the abilint section 79 80 implicit none 81 82 !Arguments ------------------------------------ 83 !scalars 84 integer,intent(in) :: choice,idir,ipw1,ispinor,itypat,lmnmax,mblkpw,ndgxdt 85 integer,intent(in) :: nffkg,nffnl,nincpw,nkpg,nlang,npw,ntens,ntypat 86 !arrays 87 integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw) 88 integer,intent(out) :: parity(nffkg) 89 real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg) 90 real(dp),intent(in) :: kpt(3) 91 real(dp),intent(out) :: ffkg(nffkg,mblkpw),kpgx(mblkpw,ntens) 92 93 !Local variables------------------------------- 94 !scalars 95 integer :: iffkg,ig,ii,ilang,ilang2,ilangx,ilmn,iln,iln0,iproj,ipw,jj 96 integer :: nffkge 97 real(dp) :: ffkg_now,kpg_x,kpg_y,kpg_z 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(iffkg,ipw)=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(iffkg,ipw)=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 345 ! End loop on ilang (ilmn) 346 end if 347 end do 348 349 !This is the number of composite projectors for the energy 350 nffkge=iffkg 351 352 !Second, treat forces : actually, this part could be rationalized, 353 !since the outcome is a multiplication by three of the number 354 !of composite projectors for the energy, while less should be needed 355 if((choice==2.or.choice==23) .and. ndgxdt/=1)then 356 do ii=1,nffkge 357 do ipw=1,nincpw 358 ffkg(iffkg+1,ipw)=ffkg(ii,ipw)*kpgx(ipw,2) 359 ffkg(iffkg+2,ipw)=ffkg(ii,ipw)*kpgx(ipw,3) 360 ffkg(iffkg+3,ipw)=ffkg(ii,ipw)*kpgx(ipw,4) 361 end do 362 parity(iffkg+1)=3-parity(ii) 363 parity(iffkg+2)=parity(iffkg+1) 364 parity(iffkg+3)=parity(iffkg+1) 365 iffkg=iffkg+3 366 end do 367 end if 368 !Note that the additional number of projectors for forces is 3*nffkge 369 370 !Third, treat first-derivative of the non-local operator 371 !with respect to an atomic displacement in one direction : 372 if(choice==2 .and. ndgxdt==1)then 373 do ii=1,nffkge 374 do ipw=1,nincpw 375 ffkg(iffkg+1,ipw)=ffkg(ii,ipw)*kpgx(ipw,idir+1) 376 end do 377 parity(iffkg+1)=3-parity(ii) 378 iffkg=iffkg+1 379 end do 380 end if 381 !Note that the additional number of projectors for this case is nffkge 382 383 384 !Fourth, treat dynamical matrices : like forces, this part could be rationalized. 385 if(choice==4)then 386 do ii=1,nffkge 387 do ipw=1,nincpw 388 kpg_x=kpgx(ipw,2) ; kpg_y=kpgx(ipw,3) ; kpg_z=kpgx(ipw,4) 389 ffkg_now=ffkg(ii,ipw) 390 ffkg(iffkg+1,ipw)=ffkg_now*kpg_x 391 ffkg(iffkg+2,ipw)=ffkg_now*kpg_y 392 ffkg(iffkg+3,ipw)=ffkg_now*kpg_z 393 ffkg(iffkg+4,ipw)=ffkg_now*kpg_x*kpg_x 394 ffkg(iffkg+5,ipw)=ffkg_now*kpg_y*kpg_y 395 ffkg(iffkg+6,ipw)=ffkg_now*kpg_z*kpg_z 396 ffkg(iffkg+7,ipw)=ffkg_now*kpg_z*kpg_y 397 ffkg(iffkg+8,ipw)=ffkg_now*kpg_z*kpg_x 398 ffkg(iffkg+9,ipw)=ffkg_now*kpg_y*kpg_x 399 end do 400 parity(iffkg+1:iffkg+3)=3-parity(ii) 401 parity(iffkg+4:iffkg+9)=parity(ii) 402 iffkg=iffkg+9 403 end do 404 end if 405 !Note that the additional number of projectors for dynamical matrices is 9*nffkge 406 407 !Treat composite projectors for the stress or 1st derivative contribution 408 !to frozen-wavefunction part of elastic tensor 409 !as well as, for ddk perturbation, the part that depend on ffnl(:,2,..) 410 if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then 411 412 iln0=0 413 do ilmn=1,lmnmax 414 if (ispinor==indlmn(6,ilmn,itypat)) then 415 iln=indlmn(5,ilmn,itypat) 416 if (iln>iln0) then 417 iln0=iln 418 ilang=1+indlmn(1,ilmn,itypat) 419 iproj=indlmn(3,ilmn,itypat) 420 if(iproj>0)then 421 ! number of unique tensor components 422 if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2 423 if(choice==5)ilangx=(ilang*(ilang+1))/2 424 425 do ii=1,ilangx 426 ! Get the starting address for the relevant tensor 427 if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 428 if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6 429 ig=ipw1 430 iffkg=iffkg+1 431 if(choice==3 .or. choice==6 .or. choice==23)then 432 do ipw=1,nincpw 433 ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 434 ig=ig+1 435 end do 436 else 437 do ipw=1,nincpw 438 ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)*& 439 & (kpgx(ipw,2)*gmet(1,idir)+ & 440 & kpgx(ipw,3)*gmet(2,idir)+ & 441 & kpgx(ipw,4)*gmet(3,idir) ) 442 ig=ig+1 443 end do 444 end if 445 if(ilang==1 .or. ilang==3)parity(iffkg)=2 446 if(ilang==2 .or. ilang==4)parity(iffkg)=1 447 if(choice==5)parity(iffkg)=3-parity(iffkg) 448 end do 449 450 ! End condition if(iproj>0) 451 end if 452 ! End condition on iln 453 end if 454 ! End condition if(ispinor=indlmn(6,...)) 455 end if 456 ! End loop on ilmn 457 end do 458 459 ! End condition of stress 460 end if 461 462 !Treat composite projectors for the 2nd derivative wrt 2 strains 463 !and wrt one strain and one atomic displacement (internal strain) 464 !contributions to frozen-wavefunction part of (generalized) elastic tensor. 465 !There are 3 sets on terms (in historical order): 466 !first, terms with ffnl(:,3,...) and rank+4 tensors. 467 !second, terms with ffnl(:,1,...) and rank+1 tensors. 468 !third, terms with ffnl(:,2,...) and rank+3 tensors. 469 470 if(choice==6)then 471 472 iln0=0 473 do ilmn=1,lmnmax 474 if (ispinor==indlmn(6,ilmn,itypat)) then 475 iln=indlmn(5,ilmn,itypat) 476 if (iln>iln0) then 477 iln0=iln 478 ilang=1+indlmn(1,ilmn,itypat) 479 iproj=indlmn(3,ilmn,itypat) 480 481 if(iproj>0)then 482 ! First set of terms 483 ! number of unique tensor components 484 ilangx=((ilang+4)*(ilang+5))/2 485 486 do ii=1,ilangx 487 ! Get the starting address for the relevant tensor 488 jj=ii+((ilang+3)*(ilang+4)*(ilang+5))/6 489 ig=ipw1 490 iffkg=iffkg+1 491 do ipw=1,nincpw 492 ffkg(iffkg,ipw)=ffnl(ig,3,ilmn,itypat)*kpgx(ipw,jj) 493 ig=ig+1 494 end do 495 if(ilang==1 .or. ilang==3)parity(iffkg)=2 496 if(ilang==2 .or. ilang==4)parity(iffkg)=1 497 end do 498 499 ! Second set of terms 500 ! number of unique tensor components 501 ilangx=((ilang+1)*(ilang+2))/2 502 503 do ii=1,ilangx 504 ! Get the starting address for the relevant tensor 505 jj=ii+((ilang)*(ilang+1)*(ilang+2))/6 506 ig=ipw1 507 iffkg=iffkg+1 508 do ipw=1,nincpw 509 ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 510 ig=ig+1 511 end do 512 if(ilang==1 .or. ilang==3)parity(iffkg)=1 513 if(ilang==2 .or. ilang==4)parity(iffkg)=2 514 end do 515 516 ! Third set of terms 517 ! number of unique tensor components 518 ilangx=((ilang+3)*(ilang+4))/2 519 520 do ii=1,ilangx 521 ! Get the starting address for the relevant tensor 522 jj=ii+((ilang+2)*(ilang+3)*(ilang+4))/6 523 ig=ipw1 524 iffkg=iffkg+1 525 do ipw=1,nincpw 526 ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj) 527 ig=ig+1 528 end do 529 if(ilang==1 .or. ilang==3)parity(iffkg)=1 530 if(ilang==2 .or. ilang==4)parity(iffkg)=2 531 end do 532 533 ! End condition if(iproj>0) 534 end if 535 ! End condition on iln 536 end if 537 ! End condition if(ispinor=indlmn(6,...)) 538 end if 539 ! End loop on ilmn 540 end do 541 542 ! End condition of 2nd strain derivatives 543 end if 544 545 !For ddk perturbation, treat the part that depend on ffnl(:,1,..) 546 !no contribution from s state 547 if(nlang>=2 .and. choice==5)then 548 iln0=0 549 do ilmn=1,lmnmax 550 if (ispinor==indlmn(6,ilmn,itypat)) then 551 iln=indlmn(5,ilmn,itypat) 552 if (iln>iln0) then 553 iln0=iln 554 ilang=1+indlmn(1,ilmn,itypat) 555 if (ilang>=2) then 556 iproj=indlmn(3,ilmn,itypat) 557 if(iproj>0)then 558 ilang2=(ilang*(ilang-1))/2 559 560 do ii=1,ilang2 561 ! Get the starting address for the relevant tensor 562 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 563 ig=ipw1 564 iffkg=iffkg+1 565 do ipw=1,nincpw 566 ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj) 567 ig=ig+1 568 end do 569 if(ilang==2 .or. ilang==4)parity(iffkg)=2 570 if(ilang==3)parity(iffkg)=1 571 end do 572 573 ! End condition if(iproj>0) 574 end if 575 ! End condition if(ilang>=2) 576 end if 577 ! End condition if(iln>iln0) 578 end if 579 ! End condition if(ispinor=indlmn(6,...)) 580 end if 581 ! End loop on ilmn 582 end do 583 ! End condition of p,d or f state 584 end if 585 586 !DEBUG 587 !write(std_out,*)' dfpt_mkffkg : exit ' 588 !ENDDEBUG 589 590 end subroutine dfpt_mkffkg