TABLE OF CONTENTS


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
189 !  Add additional tensors for strain gradients
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

  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 .

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
 765 !  Add additional tensors for strain gradients
 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