## ABINIT/dfpt_mkffkg [ Functions ]

[ Top ] [ 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.
```

PARENTS

```      opernl3,opernl4a,opernl4b
```

CHILDREN

SOURCE

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

## ABINIT/m_mkffkg [ Modules ]

[ Top ] [ Modules ]

NAME

```  m_mkffkg
```

FUNCTION

```  Copyright (C) 1998-2018 ABINIT group (DCA, XG, MT, DRH)
GNU General Public License, see ~abinit/COPYING
or http://www.gnu.org/copyleft/gpl.txt .
```

PARENTS

CHILDREN

SOURCE

```21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24
25 #include "abi_common.h"
26
27 module m_mkffkg
28
29  use defs_basis
30  use m_abicore
31
32  implicit none
33
34  private
```

## ABINIT/mkffkg [ Functions ]

[ Top ] [ 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.
```

PARENTS

```      opernl2
```

CHILDREN

SOURCE

``` 670 subroutine mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
671 &                  kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
672 &                  npw,ntens,ntypat,parity)
673
674
675 !This section has been created automatically by the script Abilint (TD).
676 !Do not modify the following lines by hand.
677 #undef ABI_FUNC
678 #define ABI_FUNC 'mkffkg'
679 !End of the abilint section
680
681  implicit none
682
683 !Arguments ------------------------------------
684 !scalars
685  integer,intent(in) :: choice,idir,ipw1,ispinor,itypat,lmnmax,mblkpw,ndgxdt
686  integer,intent(in) :: nffkg,nffnl,nincpw,nkpg,nlang,npw,ntens,ntypat
687 !arrays
688  integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw)
689  integer,intent(out) :: parity(nffkg)
690  real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg)
691  real(dp),intent(in) :: kpt(3)
692  real(dp),intent(out) :: ffkg(mblkpw,nffkg),kpgx(mblkpw,ntens)
693
694 !Local variables-------------------------------
695 !scalars
696  integer :: iffkg,ig,ii,ilang,ilang2,ilangx,ilmn,iln,iln0,iproj,ipw,jj
697  integer :: nffkge
698  real(dp) :: ffkg_now,kpg_x,kpg_y,kpg_z
699 !arrays
700
701 ! *************************************************************************
702
703  jj=0;ilangx=0
704
705 !This will be useless after all the modifications have been done
706  do ipw=1,nincpw
707    kpgx(ipw,1)=1.0d0
708  end do
709
710 !Initialize kpgx array related to tensors defined below
711  if ( nlang>=2 .or. choice==2 .or. choice==3 .or. choice==4 .or. choice==5&
712 & .or. choice==6 .or. choice==23) then
713    if (nkpg>=3) then
714      kpgx(1:nincpw,2)=kpg_k(ipw1+1:ipw1+nincpw,1)
715      kpgx(1:nincpw,3)=kpg_k(ipw1+1:ipw1+nincpw,2)
716      kpgx(1:nincpw,4)=kpg_k(ipw1+1:ipw1+nincpw,3)
717    else
718      ig=ipw1
719      do ipw=1,nincpw
720        kpgx(ipw,2)=kpt(1)+dble(kg_k(1,ig))
721        kpgx(ipw,3)=kpt(2)+dble(kg_k(2,ig))
722        kpgx(ipw,4)=kpt(3)+dble(kg_k(3,ig))
723        ig=ig+1
724      end do
725    end if
726  end if
727  if (nlang>=3 .or. choice==3 .or. choice==6 .or. choice==23) then
728 !  Define (k+G) part of rank 2 symmetric tensor (6 components), l=2
729 !  Compressed storage is 11 22 33 32 31 21
730    if (nkpg>=9) then
731      kpgx(1:nincpw,5) =kpg_k(ipw1+1:ipw1+nincpw,4)
732      kpgx(1:nincpw,6) =kpg_k(ipw1+1:ipw1+nincpw,5)
733      kpgx(1:nincpw,7) =kpg_k(ipw1+1:ipw1+nincpw,6)
734      kpgx(1:nincpw,8) =kpg_k(ipw1+1:ipw1+nincpw,7)
735      kpgx(1:nincpw,9) =kpg_k(ipw1+1:ipw1+nincpw,8)
736      kpgx(1:nincpw,10)=kpg_k(ipw1+1:ipw1+nincpw,9)
737    else
738      do ipw=1,nincpw
739        kpgx(ipw, 5) =      kpgx(ipw, 2)*kpgx(ipw, 2)
740        kpgx(ipw, 6) =      kpgx(ipw, 3)*kpgx(ipw, 3)
741        kpgx(ipw, 7) =      kpgx(ipw, 4)*kpgx(ipw, 4)
742        kpgx(ipw, 8) =      kpgx(ipw, 4)*kpgx(ipw, 3)
743        kpgx(ipw, 9) =      kpgx(ipw, 4)*kpgx(ipw, 2)
744        kpgx(ipw,10) =      kpgx(ipw, 3)*kpgx(ipw, 2)
745      end do
746    end if
747  end if
748  if (nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) .or. choice==6) then
749 !  Define (k+G) part of rank 3 symmetric tensor (10 components), l=3
750 !  Compressed storage is 111 221 331 321 311 211 222 332 322 333
751    do ipw=1,nincpw
752      kpgx(ipw,11) =     kpgx(ipw, 5)*kpgx(ipw, 2)
753      kpgx(ipw,12) =     kpgx(ipw, 6)*kpgx(ipw, 2)
754      kpgx(ipw,13) =     kpgx(ipw, 7)*kpgx(ipw, 2)
755      kpgx(ipw,14) =     kpgx(ipw, 8)*kpgx(ipw, 2)
756      kpgx(ipw,15) =     kpgx(ipw, 9)*kpgx(ipw, 2)
757      kpgx(ipw,16) =     kpgx(ipw,10)*kpgx(ipw, 2)
758      kpgx(ipw,17) =     kpgx(ipw, 6)*kpgx(ipw, 3)
759      kpgx(ipw,18) =     kpgx(ipw, 7)*kpgx(ipw, 3)
760      kpgx(ipw,19) =     kpgx(ipw, 8)*kpgx(ipw, 3)
761      kpgx(ipw,20) =     kpgx(ipw, 7)*kpgx(ipw, 4)
762    end do
763  end if
764  if (((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6) then
766 !  Define (k+G) part of rank 4 symmetric tensor (15 components), l=2
767 !  Compressed storage is 1111 2211 3311 3211 3111 2111 2221 3321 3221
768 !  3331 2222 3322 3222 3332 3333
769    do ipw=1,nincpw
770      kpgx(ipw,21) =     kpgx(ipw, 5)*kpgx(ipw, 5)
771      kpgx(ipw,22) =     kpgx(ipw, 6)*kpgx(ipw, 5)
772      kpgx(ipw,23) =     kpgx(ipw, 7)*kpgx(ipw, 5)
773      kpgx(ipw,24) =     kpgx(ipw, 8)*kpgx(ipw, 5)
774      kpgx(ipw,25) =     kpgx(ipw, 9)*kpgx(ipw, 5)
775      kpgx(ipw,26) =     kpgx(ipw,10)*kpgx(ipw, 5)
776      kpgx(ipw,27) =     kpgx(ipw, 6)*kpgx(ipw,10)
777      kpgx(ipw,28) =     kpgx(ipw, 7)*kpgx(ipw,10)
778      kpgx(ipw,29) =     kpgx(ipw, 8)*kpgx(ipw,10)
779      kpgx(ipw,30) =     kpgx(ipw, 7)*kpgx(ipw, 9)
780      kpgx(ipw,31) =     kpgx(ipw, 6)*kpgx(ipw, 6)
781      kpgx(ipw,32) =     kpgx(ipw, 7)*kpgx(ipw, 6)
782      kpgx(ipw,33) =     kpgx(ipw, 8)*kpgx(ipw, 6)
783      kpgx(ipw,34) =     kpgx(ipw, 7)*kpgx(ipw, 8)
784      kpgx(ipw,35) =     kpgx(ipw, 7)*kpgx(ipw, 7)
785    end do
786  end if
787  if (((choice==3.or.choice==23) .and. nlang>=4) .or. (choice==6 .and. nlang>=2)) then
788 !  Define (k+G) part of rank 5 symmetric tensor (21 components), l=3
789 !  Compressed storage is 11111 22111 33111 32111 31111 21111
790 !  22211 33211 32211     33311 22221 33221 32221 33321 33331
791 !  22222 33222 32222     33322 33332 33333
792    do ipw=1,nincpw
793      kpgx(ipw,36) =     kpgx(ipw,21)*kpgx(ipw, 2)
794      kpgx(ipw,37) =     kpgx(ipw,22)*kpgx(ipw, 2)
795      kpgx(ipw,38) =     kpgx(ipw,23)*kpgx(ipw, 2)
796      kpgx(ipw,39) =     kpgx(ipw,24)*kpgx(ipw, 2)
797      kpgx(ipw,40) =     kpgx(ipw,25)*kpgx(ipw, 2)
798      kpgx(ipw,41) =     kpgx(ipw,26)*kpgx(ipw, 2)
799      kpgx(ipw,42) =     kpgx(ipw,27)*kpgx(ipw, 2)
800      kpgx(ipw,43) =     kpgx(ipw,28)*kpgx(ipw, 2)
801      kpgx(ipw,44) =     kpgx(ipw,29)*kpgx(ipw, 2)
802      kpgx(ipw,45) =     kpgx(ipw,30)*kpgx(ipw, 2)
803      kpgx(ipw,46) =     kpgx(ipw,31)*kpgx(ipw, 2)
804      kpgx(ipw,47) =     kpgx(ipw,32)*kpgx(ipw, 2)
805      kpgx(ipw,48) =     kpgx(ipw,33)*kpgx(ipw, 2)
806      kpgx(ipw,49) =     kpgx(ipw,34)*kpgx(ipw, 2)
807      kpgx(ipw,50) =     kpgx(ipw,35)*kpgx(ipw, 2)
808      kpgx(ipw,51) =     kpgx(ipw,31)*kpgx(ipw, 3)
809      kpgx(ipw,52) =     kpgx(ipw,32)*kpgx(ipw, 3)
810      kpgx(ipw,53) =     kpgx(ipw,33)*kpgx(ipw, 3)
811      kpgx(ipw,54) =     kpgx(ipw,34)*kpgx(ipw, 3)
812      kpgx(ipw,55) =     kpgx(ipw,35)*kpgx(ipw, 3)
813      kpgx(ipw,56) =     kpgx(ipw,35)*kpgx(ipw, 4)
814    end do
815  end if
816  if (choice==6 .and. nlang>=3) then
817 !  Define (k+G) part of rank 6 symmetric tensor (28 components)
818 !  Compressed storage is
819 !  111111 221111 331111 321111 311111 211111 222111 332111 322111
820 !  333111 222211 332211 322211 333211 333311 222221 332221 322221
821 !  333221 333321 333331 222222 332222 322222 333222 333322 333332
822 !  333333
823    do ipw=1,nincpw
824      kpgx(ipw,57) =     kpgx(ipw,36)*kpgx(ipw, 2)
825      kpgx(ipw,58) =     kpgx(ipw,37)*kpgx(ipw, 2)
826      kpgx(ipw,59) =     kpgx(ipw,38)*kpgx(ipw, 2)
827      kpgx(ipw,60) =     kpgx(ipw,39)*kpgx(ipw, 2)
828      kpgx(ipw,61) =     kpgx(ipw,40)*kpgx(ipw, 2)
829      kpgx(ipw,62) =     kpgx(ipw,41)*kpgx(ipw, 2)
830      kpgx(ipw,63) =     kpgx(ipw,42)*kpgx(ipw, 2)
831      kpgx(ipw,64) =     kpgx(ipw,43)*kpgx(ipw, 2)
832      kpgx(ipw,65) =     kpgx(ipw,44)*kpgx(ipw, 2)
833      kpgx(ipw,66) =     kpgx(ipw,45)*kpgx(ipw, 2)
834      kpgx(ipw,67) =     kpgx(ipw,46)*kpgx(ipw, 2)
835      kpgx(ipw,68) =     kpgx(ipw,47)*kpgx(ipw, 2)
836      kpgx(ipw,69) =     kpgx(ipw,48)*kpgx(ipw, 2)
837      kpgx(ipw,70) =     kpgx(ipw,49)*kpgx(ipw, 2)
838      kpgx(ipw,71) =     kpgx(ipw,50)*kpgx(ipw, 2)
839      kpgx(ipw,72) =     kpgx(ipw,51)*kpgx(ipw, 2)
840      kpgx(ipw,73) =     kpgx(ipw,52)*kpgx(ipw, 2)
841      kpgx(ipw,74) =     kpgx(ipw,53)*kpgx(ipw, 2)
842      kpgx(ipw,75) =     kpgx(ipw,54)*kpgx(ipw, 2)
843      kpgx(ipw,76) =     kpgx(ipw,55)*kpgx(ipw, 2)
844      kpgx(ipw,77) =     kpgx(ipw,56)*kpgx(ipw, 2)
845      kpgx(ipw,78) =     kpgx(ipw,51)*kpgx(ipw, 3)
846      kpgx(ipw,79) =     kpgx(ipw,52)*kpgx(ipw, 3)
847      kpgx(ipw,80) =     kpgx(ipw,53)*kpgx(ipw, 3)
848      kpgx(ipw,81) =     kpgx(ipw,54)*kpgx(ipw, 3)
849      kpgx(ipw,82) =     kpgx(ipw,55)*kpgx(ipw, 3)
850      kpgx(ipw,83) =     kpgx(ipw,56)*kpgx(ipw, 3)
851      kpgx(ipw,84) =     kpgx(ipw,56)*kpgx(ipw, 4)
852    end do
853  end if
854  if (choice==6 .and. nlang==4) then
855 !  Define (k+G) part of rank 7 symmetric tensor (36 components)
856 !  Compressed storage is
857 !  1111111 2211111 3311111 3211111 3111111 2111111 2221111 3321111 3221111
858 !  3331111 2222111 3322111 3222111 3332111 3333111 2222211 3322211 3222211
859 !  3332211 3333211 3333311 2222221 3322221 3222221 3332221 3333221 3333321
860 !  3333331 2222222 3322222 3222222 3332222 3333222 3333322 3333332 3333333
861    do ipw=1,nincpw
862      kpgx(ipw,85) =     kpgx(ipw,57)*kpgx(ipw, 2)
863      kpgx(ipw,86) =     kpgx(ipw,58)*kpgx(ipw, 2)
864      kpgx(ipw,87) =     kpgx(ipw,59)*kpgx(ipw, 2)
865      kpgx(ipw,88) =     kpgx(ipw,60)*kpgx(ipw, 2)
866      kpgx(ipw,89) =     kpgx(ipw,61)*kpgx(ipw, 2)
867      kpgx(ipw,90) =     kpgx(ipw,62)*kpgx(ipw, 2)
868      kpgx(ipw,91) =     kpgx(ipw,63)*kpgx(ipw, 2)
869      kpgx(ipw,92) =     kpgx(ipw,64)*kpgx(ipw, 2)
870      kpgx(ipw,93) =     kpgx(ipw,65)*kpgx(ipw, 2)
871      kpgx(ipw,94) =     kpgx(ipw,66)*kpgx(ipw, 2)
872      kpgx(ipw,95) =     kpgx(ipw,67)*kpgx(ipw, 2)
873      kpgx(ipw,96) =     kpgx(ipw,68)*kpgx(ipw, 2)
874      kpgx(ipw,97) =     kpgx(ipw,69)*kpgx(ipw, 2)
875      kpgx(ipw,98) =     kpgx(ipw,70)*kpgx(ipw, 2)
876      kpgx(ipw,99) =     kpgx(ipw,71)*kpgx(ipw, 2)
877      kpgx(ipw,100) =    kpgx(ipw,72)*kpgx(ipw, 2)
878      kpgx(ipw,101) =    kpgx(ipw,73)*kpgx(ipw, 2)
879      kpgx(ipw,102) =    kpgx(ipw,74)*kpgx(ipw, 2)
880      kpgx(ipw,103) =    kpgx(ipw,75)*kpgx(ipw, 2)
881      kpgx(ipw,104) =    kpgx(ipw,76)*kpgx(ipw, 2)
882      kpgx(ipw,105) =    kpgx(ipw,77)*kpgx(ipw, 2)
883      kpgx(ipw,106) =    kpgx(ipw,78)*kpgx(ipw, 2)
884      kpgx(ipw,107) =    kpgx(ipw,79)*kpgx(ipw, 2)
885      kpgx(ipw,108) =    kpgx(ipw,80)*kpgx(ipw, 2)
886      kpgx(ipw,109) =    kpgx(ipw,81)*kpgx(ipw, 2)
887      kpgx(ipw,110) =    kpgx(ipw,82)*kpgx(ipw, 2)
888      kpgx(ipw,111) =    kpgx(ipw,83)*kpgx(ipw, 2)
889      kpgx(ipw,112) =    kpgx(ipw,84)*kpgx(ipw, 2)
890      kpgx(ipw,113) =    kpgx(ipw,78)*kpgx(ipw, 3)
891      kpgx(ipw,114) =    kpgx(ipw,79)*kpgx(ipw, 3)
892      kpgx(ipw,115) =    kpgx(ipw,80)*kpgx(ipw, 3)
893      kpgx(ipw,116) =    kpgx(ipw,81)*kpgx(ipw, 3)
894      kpgx(ipw,117) =    kpgx(ipw,82)*kpgx(ipw, 3)
895      kpgx(ipw,118) =    kpgx(ipw,83)*kpgx(ipw, 3)
896      kpgx(ipw,119) =    kpgx(ipw,84)*kpgx(ipw, 3)
897      kpgx(ipw,120) =    kpgx(ipw,84)*kpgx(ipw, 4)
898    end do
899  end if
900
901 !*****************************************************************************
902 !
903 !Packing of composite projectors in ffkg
904
905  iffkg=0
906
907 !Treat composite projectors for the energy
908  iln0=0
909  do ilmn=1,lmnmax
910    iln=indlmn(5,ilmn,itypat)
911    if (iln>iln0) then
912      iln0=iln
913      ilang=1+indlmn(1,ilmn,itypat)
914      iproj=indlmn(3,ilmn,itypat)
915      if(iproj>0)then
916        ilang2=(ilang*(ilang+1))/2
917
918        if(ilang==1)then
919 !        Treat s-component separately
920          ig=ipw1
921          iffkg=iffkg+1
922          do ipw=1,nincpw
923            ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)
924            ig=ig+1
925          end do
926          parity(iffkg)=2
927        else
928 !        Treat other components (could be made faster by treating explicitely
929 !        each angular momentum)
930          do ii=1,ilang2
931 !          Get the starting address for the relevant tensor
932            jj=ii+((ilang-1)*ilang*(ilang+1))/6
933            ig=ipw1
934            iffkg=iffkg+1
935            do ipw=1,nincpw
936              ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
937              ig=ig+1
938            end do
939            if(ilang==2 .or. ilang==4)parity(iffkg)=1
940            if(ilang==3)parity(iffkg)=2
941          end do
942        end if
943
944 !      End condition if(iproj>0)
945      end if
946 !    End loop on ilang (ilmn)
947    end if
948  end do
949
950 !This is the number of composite projectors for the energy
951  nffkge=iffkg
952
953 !Second, treat forces : actually, this part could be rationalized,
954 !since the outcome is a multiplication by three of the number
955 !of composite projectors for the energy, while less should be needed
956  if((choice==2.or.choice==23) .and. ndgxdt/=1)then
957    do ii=1,nffkge
958      do ipw=1,nincpw
959        ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,2)
960        ffkg(ipw,iffkg+2)=ffkg(ipw,ii)*kpgx(ipw,3)
961        ffkg(ipw,iffkg+3)=ffkg(ipw,ii)*kpgx(ipw,4)
962      end do
963      parity(iffkg+1)=3-parity(ii)
964      parity(iffkg+2)=parity(iffkg+1)
965      parity(iffkg+3)=parity(iffkg+1)
966      iffkg=iffkg+3
967    end do
968  end if
969 !Note that the additional number of projectors for forces is 3*nffkge
970
971 !Third, treat first-derivative of the non-local operator
972 !with respect to an atomic displacement in one direction :
973  if(choice==2 .and. ndgxdt==1)then
974    do ii=1,nffkge
975      do ipw=1,nincpw
976        ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,idir+1)
977      end do
978      parity(iffkg+1)=3-parity(ii)
979      iffkg=iffkg+1
980    end do
981  end if
982 !Note that the additional number of projectors for this case is nffkge
983
984
985 !Fourth, treat dynamical matrices : like forces, this part could be rationalized.
986  if(choice==4)then
987    do ii=1,nffkge
988      do ipw=1,nincpw
989        kpg_x=kpgx(ipw,2) ; kpg_y=kpgx(ipw,3) ; kpg_z=kpgx(ipw,4)
990        ffkg_now=ffkg(ipw,ii)
991        ffkg(ipw,iffkg+1)=ffkg_now*kpg_x
992        ffkg(ipw,iffkg+2)=ffkg_now*kpg_y
993        ffkg(ipw,iffkg+3)=ffkg_now*kpg_z
994        ffkg(ipw,iffkg+4)=ffkg_now*kpg_x*kpg_x
995        ffkg(ipw,iffkg+5)=ffkg_now*kpg_y*kpg_y
996        ffkg(ipw,iffkg+6)=ffkg_now*kpg_z*kpg_z
997        ffkg(ipw,iffkg+7)=ffkg_now*kpg_z*kpg_y
998        ffkg(ipw,iffkg+8)=ffkg_now*kpg_z*kpg_x
999        ffkg(ipw,iffkg+9)=ffkg_now*kpg_y*kpg_x
1000      end do
1001      parity(iffkg+1:iffkg+3)=3-parity(ii)
1002      parity(iffkg+4:iffkg+9)=parity(ii)
1003      iffkg=iffkg+9
1004    end do
1005  end if
1006 !Note that the additional number of projectors for dynamical matrices is 9*nffkge
1007
1008 !Treat composite projectors for the stress or 1st derivative contribution
1009 !to frozen-wavefunction part of elastic tensor
1010 !as well as, for ddk perturbation, the part that depend on ffnl(:,2,..)
1011  if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then
1012
1013    iln0=0
1014    do ilmn=1,lmnmax
1015      if (ispinor==indlmn(6,ilmn,itypat)) then
1016        iln=indlmn(5,ilmn,itypat)
1017        if (iln>iln0) then
1018          iln0=iln
1019          ilang=1+indlmn(1,ilmn,itypat)
1020          iproj=indlmn(3,ilmn,itypat)
1021          if(iproj>0)then
1022 !          number of unique tensor components
1023            if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2
1024            if(choice==5)ilangx=(ilang*(ilang+1))/2
1025
1026            do ii=1,ilangx
1027 !            Get the starting address for the relevant tensor
1028              if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
1029              if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6
1030              ig=ipw1
1031              iffkg=iffkg+1
1032              if(choice==3 .or. choice==6 .or. choice==23)then
1033                do ipw=1,nincpw
1034                  ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)
1035                  ig=ig+1
1036                end do
1037              else
1038                do ipw=1,nincpw
1039                  ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)*&
1040 &                 (kpgx(ipw,2)*gmet(1,idir)+ &
1041 &                 kpgx(ipw,3)*gmet(2,idir)+ &
1042 &                 kpgx(ipw,4)*gmet(3,idir) )
1043                  ig=ig+1
1044                end do
1045              end if
1046              if(ilang==1 .or. ilang==3)parity(iffkg)=2
1047              if(ilang==2 .or. ilang==4)parity(iffkg)=1
1048              if(choice==5)parity(iffkg)=3-parity(iffkg)
1049            end do
1050
1051 !          End condition if(iproj>0)
1052          end if
1053
1054 !        End loop on ilang (ilmn)
1055        end if
1056      end if
1057    end do
1058
1059 !  End condition of stress
1060  end if
1061
1062 !Treat composite projectors for the 2nd derivative wrt 2 strains
1063 !and wrt one strain and one atomic displacement (internal strain)
1064 !contributions to frozen-wavefunction part of (generalized) elastic tensor.
1065 !There are 3 sets on terms (in historical order):
1066 !first,  terms with ffnl(:,3,...) and rank+4 tensors.
1067 !second, terms with ffnl(:,1,...) and rank+1 tensors.
1068 !third,  terms with ffnl(:,2,...) and rank+3 tensors.
1069
1070  if(choice==6)then
1071
1072    iln0=0
1073    do ilmn=1,lmnmax
1074      if (ispinor==indlmn(6,ilmn,itypat)) then
1075        iln=indlmn(5,ilmn,itypat)
1076        if (iln>iln0) then
1077          iln0=iln
1078          ilang=1+indlmn(1,ilmn,itypat)
1079          iproj=indlmn(3,ilmn,itypat)
1080          if(iproj>0)then
1081 !          First set of terms
1082 !          number of unique tensor components
1083            ilangx=((ilang+4)*(ilang+5))/2
1084
1085            do ii=1,ilangx
1086 !            Get the starting address for the relevant tensor
1087              jj=ii+((ilang+3)*(ilang+4)*(ilang+5))/6
1088              ig=ipw1
1089              iffkg=iffkg+1
1090              do ipw=1,nincpw
1091                ffkg(ipw,iffkg)=ffnl(ig,3,ilmn,itypat)*kpgx(ipw,jj)
1092                ig=ig+1
1093              end do
1094              if(ilang==1 .or. ilang==3)parity(iffkg)=2
1095              if(ilang==2 .or. ilang==4)parity(iffkg)=1
1096            end do
1097
1098 !          Second set of terms
1099 !          number of unique tensor components
1100            ilangx=((ilang+1)*(ilang+2))/2
1101
1102            do ii=1,ilangx
1103 !            Get the starting address for the relevant tensor
1104              jj=ii+((ilang)*(ilang+1)*(ilang+2))/6
1105              ig=ipw1
1106              iffkg=iffkg+1
1107              do ipw=1,nincpw
1108                ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
1109                ig=ig+1
1110              end do
1111              if(ilang==1 .or. ilang==3)parity(iffkg)=1
1112              if(ilang==2 .or. ilang==4)parity(iffkg)=2
1113            end do
1114
1115 !          Third set of terms
1116 !          number of unique tensor components
1117            ilangx=((ilang+3)*(ilang+4))/2
1118
1119            do ii=1,ilangx
1120 !            Get the starting address for the relevant tensor
1121              jj=ii+((ilang+2)*(ilang+3)*(ilang+4))/6
1122              ig=ipw1
1123              iffkg=iffkg+1
1124              do ipw=1,nincpw
1125                ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)
1126                ig=ig+1
1127              end do
1128              if(ilang==1 .or. ilang==3)parity(iffkg)=1
1129              if(ilang==2 .or. ilang==4)parity(iffkg)=2
1130            end do
1131
1132 !          End condition if(iproj>0)
1133          end if
1134 !        End loop on ilang (ilmn)
1135        end if
1136      end if
1137    end do
1138
1139 !  End condition of 2nd strain derivatives
1140  end if
1141
1142 !For ddk perturbation, treat the part that depend on ffnl(:,1,..)
1143 !no contribution from s state
1144  if(nlang>=2 .and. choice==5)then
1145    iln0=0
1146    do ilmn=1,lmnmax
1147      if (ispinor==indlmn(6,ilmn,itypat)) then
1148        iln=indlmn(5,ilmn,itypat)
1149        if (iln>iln0) then
1150          iln0=iln
1151          ilang=1+indlmn(1,ilmn,itypat)
1152          if (ilang>=2) then
1153            iproj=indlmn(3,ilmn,itypat)
1154            if(iproj>0)then
1155              ilang2=(ilang*(ilang-1))/2
1156
1157              do ii=1,ilang2
1158 !              Get the starting address for the relevant tensor
1159                jj=ii+((ilang-2)*(ilang-1)*ilang)/6
1160                ig=ipw1
1161                iffkg=iffkg+1
1162                do ipw=1,nincpw
1163                  ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
1164                  ig=ig+1
1165                end do
1166                if(ilang==2 .or. ilang==4)parity(iffkg)=2
1167                if(ilang==3)parity(iffkg)=1
1168              end do
1169
1170 !            End condition if(iproj>0)
1171            end if
1172 !          End loop on ilang>=2
1173          end if
1174        end if
1175      end if
1176    end do
1177 !  End condition of p,d or f state
1178  end if
1179
1180 end subroutine mkffkg
```