TABLE OF CONTENTS


ABINIT/m_mkffnl [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mkffnl

FUNCTION

 Make FFNL, nonlocal form factors, for each type of atom up to ntypat
 and for each angular momentum.

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, 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

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_mkffnl
29 
30  use defs_basis
31  use m_abicore
32  use m_errors
33  use m_splines
34 
35  use m_time,     only : timab
36  use m_kg,       only : mkkin
37 
38  implicit none
39 
40  private

ABINIT/mkffnl [ Functions ]

[ Top ] [ Functions ]

NAME

 mkffnl

FUNCTION

 Make FFNL, nonlocal form factors, for each type of atom up to ntypat
 and for each angular momentum.
 When Legendre polynomials are used in the application of the
   nonlocal operator, FFNLs depend on (l,n) components; in this
   case, form factors are real and divided by |k+G|^l;
 When spherical harmonics are used, FFNLs depend on (l,m,n)
   components; in this case, form factors are multiplied by Ylm(k+G).

INPUTS

  dimekb=second dimension of ekb (see ekb)
  dimffnl=second dimension of ffnl (1+number of derivatives)
  ekb(dimekb,ntypat*(1-usepaw))=(Real) Kleinman-Bylander energies (hartree)
                                ->NORM-CONSERVING PSPS ONLY
  ffspl(mqgrid,2,lnmax,ntypat)=form factors and spline fit to 2nd derivative
  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  gprimd(3,3)=dimensional reciprocal space primitive translations
  ider=0=>no derivative wanted; 1=>1st derivative wanted; 2=>1st and 2nd derivatives wanted
  idir=ONLY WHEN YLMs ARE USED:
       When 1st derivative has to be computed:  (see more info below)
       - Determine the direction(s) of the derivatives(s)
       - Determine the set of coordinates (reduced or cartesians)
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,spin for i=ln  (if useylm=0)
                                                     or i=lmn (if useylm=1)
  kg(3,npw)=integer coordinates of planewaves in basis sphere for this k point.
  kpg(npw,nkpg)= (k+G) components (only if useylm=1)
  kpt(3)=reduced coordinates of k point
  lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
        =if useylm=0, max number of (l,n)   comp. over all type of psps
  lnmax=max. number of (l,n) components over all type of psps
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mqgrid=size of q (or |G|) grid for f(q)
  nkpg=second dimension of kpg_k (0 if useylm=0)
  npw=number of planewaves in basis sphere
  ntypat=number of types of atoms
  usepaw= 0 for non paw calculation; =1 for paw calculation
  pspso(ntypat)=spin-orbit characteristics for each atom type (1, 2, or 3)
  qgrid(mqgrid)=uniform grid of q values from 0 to qmax
  rmet(3,3)=real space metric (bohr**2)
  useylm=governs the way the nonlocal operator is to be applied:
         1=using Ylm, 0=using Legendre polynomials
  ylm   (npw,mpsang*mpsang*useylm)=real spherical harmonics for each G and k point
  ylm_gr(npw,3,mpsang*mpsang*useylm)=gradients of real spherical harmonics wrt (k+G)

OUTPUT

  ffnl(npw,dimffnl,lmnmax,ntypat)=described below

NOTES

  Uses spline fit ffspl provided by Numerical Recipes spline subroutine.
  Form factor $f_l(q)$ is defined by
   \begin{equation}
  \textrm{f}_l(q)=\frac{1}{dvrms} \int_0^\infty [j_l(2 \pi r q) u_l(r) dV(r) r dr]
   \end{equation}
   where u_l(r)=reference state wavefunction, dV(r)=nonlocal psp
   correction, j_l(arg)=spherical Bessel function for angular momentum l,
   and
   \begin{equation}
    \textrm{dvrms} =  \int_0^\infty [(u_l(r) dV(r))^2 dr])^{1/2}
   \end{equation}
   which is square root of mean square dV, i.e.
     $ (\langle (dV)^2 \rangle)^{1/2} $ .
   This routine is passed f_l(q) in spline form in the array ffspl and then
   constructs the values of $f_l(q)$ on the relevant (k+G) in array ffnl.
   The evaluation of the integrals defining ffspl was done in mkkbff.

  Delivers the following (for each atom type t, or itypat):
   --------------------------
   Using Legendre polynomials in the application of nl operator:
     ffnl are real.
     ffnl(ig,1,(l,0,n),itypat) $= f_ln(k+G)/|k+G|^l $
     === if ider>=1
       ffnl(ig,2,(l,0,n),itypat) $=(fprime_ln(k+G)-l*f_ln(k+G)/|k+G|)/|k+G|^(l+1) $
     === if ider=2
       ffnl(ig,3,(l,0,n),itypat) $=(fprimeprime_ln(k+G)-(2l+1)*fprime_ln(k+G)/|k+G|
                                   +l(l+2)*f_ln(k+G)/|k+G|**2)/|k+G|^(l+2)
   --------------------------
   Using spherical harmonics in the application of nl operator:
     ffnl are real (we use REAL spherical harmonics).
     ffnl(ig,1,(l,m,n),itypat) = ffnl_1
                              $= f_ln(k+G) * Y_lm(k+G) $
     === if ider>=1
     --if (idir==0)
       ffnl(ig,1+i,(l,m,n),itypat) = dffnl_i = 3 reduced coord. of d(ffnl_1)/dK^cart
         $= fprime_ln(k+G).Y_lm(k+G).(k+G)^red_i/|k+G|+f_ln(k+G).(dY_lm/dK^cart)^red_i $
         for i=1..3
     --if (1<=idir<=3)
       ffnl(ig,2,(l,m,n),itypat)= cart. coordinate idir of d(ffnl_1)/dK^red
                                = Sum_(mu,nu) [ Gprim(mu,idir) Gprim(mu,nu) dffnl_nu ]
     --if (idir==4)
       ffnl(ig,1+i,(l,m,n),itypat)= 3 cart. coordinates of d(ffnl_1)/dK^red
                                  = Sum_(mu,nu) [ Gprim(mu,i) Gprim(mu,nu) dffnl_nu ]
     --if (-6<idir<-1)
       ffnl(ig,2,(l,m,n),itypat)=1/2 [d(ffnl)/dK^cart_mu K^cart_nu + d(ffnl)/dK^cart_nu K^cart_mu]
                                with d(ffnl)/dK^cart_i = Sum_nu [ Gprim(nu,i) dffnl_nu ]
                                for |idir|->(mu,nu) (1->11,2->22,3->33,4->32,5->31,6->21)
     --if (idir==-7)
       ffnl(ig,2:7,(l,m,n),itypat)=1/2 [d(ffnl)/dK^cart_mu K^cart_nu + d(ffnl)/dK^cart_nu K^cart_mu]
                                with d(ffnl)/dK^cart_i = Sum_nu [ Gprim(nu,i) dffnl_nu ]
                                for all (mu,nu) (6 independant terms)
     === if ider==2
     --if (idir==0)
       ffnl(ig,4+i,(l,m,n),itypat) = d2ffnl_mu,nu = 6 reduced coord. of d2(ffnl_1)/dK^cart.dK^cart
        for all i=(mu,nu) (6 independant terms)
     --if (idir==4)
       ffnl(ig,4+i,(l,m,n),itypat) = d2ffnl_i =6 cart. coordinates of d2(ffnl_1)/dK^red.dK^red
        for all i=(mu,nu) (6 independant terms)
        = Sum_(mu1,mu2,mu3,mu4) [ Gprim(mu1,mu) Gprim(mu2,nu) Gprim(mu1,mu3) Gprim(mu2,mu4) d2ffnl_mu3,mu4 ]
   --------------------------

  1) l may be 0, 1, 2, or 3 in this version.

  2) Norm-conserving psps : only FFNL for which ekb is not zero are calculated.

  3) Each expression above approaches a constant as $|k+G| \rightarrow 0 $.
     In the cases where $|k+G|$ is in the denominator, there is always a
     factor of $(k+G)_mu$ multiplying the ffnl term where it is actually used,
     so that we may replace the ffnl term by any constant when $|k+G| = 0$.
     Below we replace 1/0 by 1/tol10, thus creating an arbitrary constant
     which will later be multiplied by 0.

TODO

  Some parts can be rewritten with BLAS1 calls.

PARENTS

      ctocprj,d2frnl,dfpt_nsteltwf,dfpt_nstpaw,dfpt_nstwf,dfpt_rhofermi
      dfptnl_resp,energy,fock2ACE,forstrnps,getgh1c,ks_ddiago,m_io_kss
      m_shirley,m_vkbr,m_wfd,nonlop_test,vtorho

CHILDREN

      mkkin,splfit,timab

SOURCE

186 subroutine mkffnl(dimekb,dimffnl,ekb,ffnl,ffspl,gmet,gprimd,ider,idir,indlmn,&
187 &                   kg,kpg,kpt,lmnmax,lnmax,mpsang,mqgrid,nkpg,npw,ntypat,pspso,&
188 &                   qgrid,rmet,usepaw,useylm,ylm,ylm_gr)
189 
190 
191 !This section has been created automatically by the script Abilint (TD).
192 !Do not modify the following lines by hand.
193 #undef ABI_FUNC
194 #define ABI_FUNC 'mkffnl'
195 !End of the abilint section
196 
197  implicit none
198 
199 !Arguments ------------------------------------
200 !scalars
201  integer,intent(in) :: dimekb,dimffnl,ider,idir,lmnmax,lnmax,mpsang,mqgrid,nkpg
202  integer,intent(in) :: npw,ntypat,usepaw,useylm
203 !arrays
204  integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg(3,npw),pspso(ntypat)
205  real(dp),intent(in) :: ekb(dimekb,ntypat*(1-usepaw))
206  real(dp),intent(in) :: ffspl(mqgrid,2,lnmax,ntypat),gmet(3,3),gprimd(3,3)
207  real(dp),intent(in) :: kpg(npw,nkpg),kpt(3),qgrid(mqgrid),rmet(3,3)
208  real(dp),intent(in) :: ylm(:,:),ylm_gr(:,:,:)
209  real(dp),intent(out) :: ffnl(npw,dimffnl,lmnmax,ntypat)
210 
211 !Local variables-------------------------------
212 !scalars
213  integer :: ider_tmp,iffnl,ig,ig0,il,ilm,ilmn,iln,iln0,im,itypat,mu,mua,mub,nlmn,nu,nua,nub
214  real(dp),parameter :: renorm_factor=0.5d0/pi**2,tol_norm=tol10
215  real(dp) :: ecut,ecutsm,effmass_free,fact,kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3,rmetab,yp1
216  logical :: testnl=.false.
217  character(len=500) :: message
218 !arrays
219  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
220  integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
221  real(dp) :: rprimd(3,3),tsec(2)
222  real(dp),allocatable :: dffnl_cart(:,:),dffnl_red(:,:),dffnl_tmp(:)
223  real(dp),allocatable :: d2ffnl_cart(:,:),d2ffnl_red(:,:),d2ffnl_tmp(:)
224  real(dp),allocatable :: kpgc(:,:),kpgn(:,:),kpgnorm(:),kpgnorm_inv(:),wk_ffnl1(:)
225  real(dp),allocatable :: wk_ffnl2(:),wk_ffnl3(:),wk_ffspl(:,:)
226 
227 ! *************************************************************************
228 
229  DBG_ENTER("COLL")
230 
231 !Keep track of time spent in mkffnl
232  call timab(16,1,tsec)
233 
234 !Compatibility tests
235  if (mpsang>4) then
236    write(message,'(a,i0,a,a)')&
237 &   'Called with mpsang > 4, =',mpsang,ch10,&
238 &   'This subroutine will not accept lmax+1 > 4.'
239    MSG_BUG(message)
240  end if
241  if (idir<-7.or.idir>4) then
242    MSG_BUG('Called with idir<-6 or idir>4 !')
243  end if
244  if (useylm==0) then
245    iffnl=1+ider
246  else
247    iffnl=1
248    if (ider>=1) then
249      if (idir==0) iffnl=iffnl+3
250      if (idir/=0) iffnl=iffnl+1
251      if (idir==4) iffnl=iffnl+2
252      if (idir==-7) iffnl=iffnl+5
253    end if
254    if (ider==2) then
255      if (idir==0) iffnl=iffnl+6
256      if (idir==4) iffnl=iffnl+6
257    end if
258  end if
259  if (iffnl/=dimffnl) then
260    write(message,'(2(a,i1),a,i2)') 'Incompatibility between ider, idir and dimffnl : ider = ',ider,&
261 &                                    ' idir = ',idir,' dimffnl = ',dimffnl
262    MSG_BUG(message)
263  end if
264  if (useylm==1) then
265    ABI_CHECK(size(ylm,1)==npw,'BUG: wrong ylm size (1)')
266    ABI_CHECK(size(ylm,2)==mpsang**2,'BUG: wrong ylm size (2)')
267    if(ider>0)then
268      ABI_CHECK(size(ylm_gr,1)==npw,'BUG: wrong ylm_gr size (1)')
269      ABI_CHECK(size(ylm_gr,2)>=3+6*(ider/2),'BUG: wrong ylm_gr size (2)')
270      ABI_CHECK(size(ylm_gr,3)==mpsang**2,'BUG: wrong ylm_gr size (3)')
271    end if
272  end if
273 !Get (k+G) and |k+G|:
274  ABI_ALLOCATE(kpgnorm,(npw))
275  ABI_ALLOCATE(kpgnorm_inv,(npw))
276  ig0=-1 ! index of |k+g|=0 vector
277 
278  if (useylm==1) then
279    ABI_ALLOCATE(kpgc,(npw,3))
280    if (ider>=1) then
281      ABI_ALLOCATE(kpgn,(npw,3))
282    end if
283    if (nkpg<3) then
284 !$OMP PARALLEL DO PRIVATE(ig,kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3)
285      do ig=1,npw
286        kpg1=kpt(1)+dble(kg(1,ig))
287        kpg2=kpt(2)+dble(kg(2,ig))
288        kpg3=kpt(3)+dble(kg(3,ig))
289        kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3)
290        kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3)
291        kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3)
292        kpgc(ig,1)=kpgc1
293        kpgc(ig,2)=kpgc2
294        kpgc(ig,3)=kpgc3
295        kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
296        if (kpgnorm(ig)<=tol_norm) ig0=ig
297        if (ider>=1) then
298          kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
299          kpgn(ig,1)=kpg1*kpgnorm_inv(ig)
300          kpgn(ig,2)=kpg2*kpgnorm_inv(ig)
301          kpgn(ig,3)=kpg3*kpgnorm_inv(ig)
302        end if
303      end do
304    else
305 !$OMP PARALLEL DO PRIVATE(ig,kpgc1,kpgc2,kpgc3)
306      do ig=1,npw
307        kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
308        kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
309        kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
310        kpgc(ig,1)=kpgc1
311        kpgc(ig,2)=kpgc2
312        kpgc(ig,3)=kpgc3
313        kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
314        if (kpgnorm(ig)<=tol_norm) ig0=ig
315        if (ider>=1) then
316          kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
317          kpgn(ig,1:3)=kpg(ig,1:3)*kpgnorm_inv(ig)
318        end if
319      end do
320    end if
321  else
322    if (nkpg<3) then
323      ecut=huge(0.0d0)*0.1d0;ecutsm=zero;effmass_free=one
324 !    Note that with ecutsm=0, the right kinetic energy is computed
325      call mkkin(ecut,ecutsm,effmass_free,gmet,kg,kpgnorm,kpt,npw,0,0)
326 !$OMP PARALLEL DO
327      do ig=1,npw
328        kpgnorm(ig)=sqrt(renorm_factor*kpgnorm(ig))
329        kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
330        if (kpgnorm(ig)<=tol_norm) ig0=ig
331      end do
332    else
333 !$OMP PARALLEL DO PRIVATE(ig,kpgc1,kpgc2,kpgc3)
334      do ig=1,npw
335        kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
336        kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
337        kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
338        kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
339        kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
340        if (kpgnorm(ig)<=tol_norm) ig0=ig
341      end do
342    end if
343  end if
344 
345 !Need rprimd in some cases
346  if (ider>=1.and.useylm==1.and.ig0>0) then
347    do mu=1,3
348      do nu=1,3
349        rprimd(mu,nu)=gprimd(mu,1)*rmet(1,nu)+gprimd(mu,2)*rmet(2,nu)+gprimd(mu,3)*rmet(3,nu)
350      end do
351    end do
352  end if
353 
354 !Allocate several temporary arrays
355  ABI_ALLOCATE(wk_ffnl1,(npw))
356  ABI_ALLOCATE(wk_ffnl2,(npw))
357  ABI_ALLOCATE(wk_ffnl3,(npw))
358  ABI_ALLOCATE(wk_ffspl,(mqgrid,2))
359  if (ider>=1.and.useylm==1) then
360    ABI_ALLOCATE(dffnl_red,(npw,3))
361    if (idir/=0)  then
362      ABI_ALLOCATE(dffnl_cart,(npw,3))
363    end if
364    if (idir>0)   then
365      ABI_ALLOCATE(dffnl_tmp,(npw))
366    end if
367  end if
368  if (ider>=2.and.useylm==1) then
369    ABI_ALLOCATE(d2ffnl_red,(npw,6))
370    if (idir==4)  then
371      ABI_ALLOCATE(d2ffnl_cart,(npw,6))
372      ABI_ALLOCATE(d2ffnl_tmp,(npw))
373    end if
374  end if
375 
376 !Loop over types of atoms
377  do itypat=1,ntypat
378 
379 !  Loop over (l,m,n) values
380    iln0=0;nlmn=count(indlmn(3,:,itypat)>0)
381    ffnl(:,:,:,itypat)=zero
382    do ilmn=1,nlmn
383      il=indlmn(1,ilmn,itypat)
384      im=indlmn(2,ilmn,itypat)
385      ilm =indlmn(4,ilmn,itypat)
386      iln =indlmn(5,ilmn,itypat)
387      iffnl=ilmn;if (useylm==0) iffnl=iln
388 !    Special case : we enter the loop in case of spin-orbit calculation
389 !       even if the psp has no spin-orbit component.
390      if ((indlmn(6,ilmn,itypat)==1).or.(pspso(itypat)/=0)) then
391 
392 !      Compute FFNL only if ekb>0 or paw
393        if (usepaw==1) testnl=.true.
394        if (usepaw==0) testnl=(abs(ekb(iln,itypat))>tol_norm)
395        if (testnl) then
396 
397 !        Store form factors (from ffspl)
398 !        -------------------------------
399          if (iln>iln0) then
400            do ig=1,mqgrid
401              wk_ffspl(ig,1)=ffspl(ig,1,iln,itypat)
402              wk_ffspl(ig,2)=ffspl(ig,2,iln,itypat)
403            end do
404            ider_tmp=min(ider,1)
405            call splfit(qgrid,wk_ffnl2,wk_ffspl,ider_tmp,kpgnorm,wk_ffnl1,mqgrid,npw)
406            if(ider==2) then
407              call splfit(qgrid,wk_ffnl3,wk_ffspl,ider,kpgnorm,wk_ffnl1,mqgrid,npw)
408            end if
409          end if
410 
411 !        Store FFNL and FFNL derivatives
412 !        -------------------------------
413 
414 !        =========================================================================
415 !        A-USE OF SPHER. HARMONICS IN APPLICATION OF NL OPERATOR:
416 !        ffnl(K,l,m,n)=fnl(K).Ylm(K)
417 !        --if (idir==0)
418 !        ffnl_prime(K,1:3,l,m,n)=3 reduced coordinates of d(ffnl)/dK^cart
419 !        =fnl_prime(K).Ylm(K).K^red_i/|K|+fnl(K).(dYlm/dK^cart)^red_i
420 !        --if (0<idir<4)
421 !        ffnl_prime(K,l,m,n)=cart. coordinate idir of d(ffnl)/dK^red
422 !        --if (idir==4)
423 !        ffnl_prime(K,l,m,n)=3 cart. coordinates of d(ffnl)/dK^red
424 !        --if (-7<=idir<0) - |idir|=(mu,nu) (1->11,2->22,3->33,4->32,5->31,6->21)
425 !        ffnl_prime(K,l,m,n)=1/2 [d(ffnl)/dK^cart_mu K^cart_nu + d(ffnl)/dK^cart_nu K^cart_mu]
426 !        ffnl_prime_prime(K,l,m,n)=6 reduced coordinates of d2(ffnl)/dK^cart.dK^cart
427 
428          if (useylm==1) then
429 !$OMP PARALLEL DO
430            do ig=1,npw
431              ffnl(ig,1,iffnl,itypat)=ylm(ig,ilm)*wk_ffnl1(ig)
432            end do
433 
434            if (ider>=1) then
435 !$OMP PARALLEL DO COLLAPSE(2)
436              do mu=1,3
437                do ig=1,npw
438                  dffnl_red(ig,mu)=ylm(ig,ilm)*wk_ffnl2(ig)*kpgn(ig,mu)+ylm_gr(ig,mu,ilm)*wk_ffnl1(ig)
439                end do
440              end do
441              !Special cases |k+g|=0
442              if (ig0>0) then
443                do mu=1,3
444                  dffnl_red(ig0,mu)=zero
445                  if (il==1) then
446                    !Retrieve 1st-deriv. of ffnl at q=zero according to spline routine
447                    yp1=(wk_ffspl(2,1)-wk_ffspl(1,1))/qgrid(2)-sixth*qgrid(2)*(two*wk_ffspl(1,2)+wk_ffspl(2,2))
448                    fact=yp1*sqrt(three/four_pi)
449                    if (im==-1) dffnl_red(ig0,mu)=fact*rprimd(2,mu)
450                    if (im== 0) dffnl_red(ig0,mu)=fact*rprimd(3,mu)
451                    if (im==+1) dffnl_red(ig0,mu)=fact*rprimd(1,mu)
452                  end if
453                end do
454              end if
455              if (idir==0) then
456 !$OMP PARALLEL DO COLLAPSE(2)
457                do mu=1,3
458                  do ig=1,npw
459                    ffnl(ig,1+mu,iffnl,itypat)=dffnl_red(ig,mu)
460                  end do
461                end do
462              else
463                dffnl_cart=zero
464 !$OMP PARALLEL DO COLLAPSE(2)
465                do mu=1,3
466                  do ig=1,npw
467                    do nu=1,3
468                      dffnl_cart(ig,mu)=dffnl_cart(ig,mu)+dffnl_red(ig,nu)*gprimd(mu,nu)
469                    end do
470                  end do
471                end do
472                if (idir>=1.and.idir<=3) then
473                  dffnl_tmp=zero
474 !$OMP PARALLEL PRIVATE(nu,ig)
475 !$OMP DO
476                  do ig=1,npw
477                    do nu=1,3
478                      dffnl_tmp(ig)=dffnl_tmp(ig) + dffnl_cart(ig,nu)*gprimd(nu,idir)
479                    end do
480                  end do
481 !$OMP END DO
482 !$OMP WORKSHARE
483                  ffnl(:,2,iffnl,itypat)=dffnl_tmp(:)
484 !$OMP END WORKSHARE
485 !$OMP END PARALLEL
486                else if (idir==4) then
487                  do mu=1,3
488 !$OMP PARALLEL PRIVATE(nu,ig)
489 !$OMP WORKSHARE
490                    dffnl_tmp=zero
491 !$OMP END WORKSHARE
492 !$OMP DO
493                    do ig=1,npw
494                      do nu=1,3
495                        dffnl_tmp(ig)=dffnl_tmp(ig) + dffnl_cart(ig,nu)*gprimd(nu,mu)
496                      end do
497                    end do
498 !$OMP END DO
499 !$OMP WORKSHARE
500                    ffnl(:,1+mu,iffnl,itypat)=dffnl_tmp(:)
501 !$OMP END WORKSHARE
502 !$OMP END PARALLEL
503                  end do
504                else if (idir/=-7) then
505                  mu=abs(idir);mua=alpha(mu);mub=beta(mu)
506 !$OMP PARALLEL DO
507                  do ig=1,npw
508                    ffnl(ig,2,iffnl,itypat)=0.5d0* &
509 &                   (dffnl_cart(ig,mua)*kpgc(ig,mub) &
510 &                   +dffnl_cart(ig,mub)*kpgc(ig,mua))
511                  end do
512                else if (idir==-7) then
513 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(mua, mub)
514                  do mu=1,6
515                    do ig=1,npw
516                      mua=alpha(mu);mub=beta(mu)
517                      ffnl(ig,1+mu,iffnl,itypat)=0.5d0* &
518 &                     (dffnl_cart(ig,mua)*kpgc(ig,mub) &
519 &                     +dffnl_cart(ig,mub)*kpgc(ig,mua))
520                    end do
521                  end do
522                end if
523              end if
524            end if
525 
526            if (ider==2) then
527              do mu=1,6
528                mua=alpha(mu);mub=beta(mu)
529                rmetab=rmet(mua,mub)
530 !$OMP PARALLEL DO
531                do ig=1,npw
532                  d2ffnl_red(ig,mu)= &
533 &                 ylm_gr(ig,3+mu,ilm)*wk_ffnl1(ig) &
534 &                 + (rmetab-kpgn(ig,mua)*kpgn(ig,mub))*ylm(ig,ilm)*wk_ffnl2(ig)*kpgnorm_inv(ig) &
535 &                 + ylm(ig,ilm)*kpgn(ig,mua)*kpgn(ig,mub)*wk_ffnl3(ig) &
536 &                 + (ylm_gr(ig,mua,ilm)*kpgn(ig,mub)+ylm_gr(ig,mub,ilm)*kpgn(ig,mua))*wk_ffnl2(ig)
537                end do
538                !Special cases |k+g|=0
539                if (ig0>0) then
540                  d2ffnl_red(ig0,mu)=zero
541                  if (il==0) then
542                    d2ffnl_red(ig0,mu)=wk_ffspl(1,2)*rmetab/sqrt(four_pi)
543                  end if
544                  if (il==2) then
545                    fact=wk_ffspl(1,2)*quarter*sqrt(15._dp/pi)
546                    if (im==-2) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(2,mub)+rprimd(2,mua)*rprimd(1,mub))
547                    if (im==-1) d2ffnl_red(ig0,mu)=fact*(rprimd(2,mua)*rprimd(3,mub)+rprimd(3,mua)*rprimd(2,mub))
548                    if (im==+1) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(3,mub)+rprimd(3,mua)*rprimd(1,mub))
549                    if (im==+2) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(1,mub)-rprimd(2,mua)*rprimd(2,mub))
550                    if (im== 0) d2ffnl_red(ig0,mu)=(fact/sqrt3)*(two*rprimd(3,mua)*rprimd(3,mub) &
551 &                   -rprimd(1,mua)*rprimd(1,mub)-rprimd(2,mua)*rprimd(2,mub))
552                  end if
553                end if
554              end do
555              if (idir==0) then
556 !$OMP PARALLEL DO COLLAPSE(2)
557                do mu=1,6
558                  do ig=1,npw
559                    ffnl(ig,4+mu,iffnl,itypat)=d2ffnl_red(ig,mu)
560                  end do
561                end do
562              else if (idir==4) then
563                d2ffnl_cart=zero
564 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(mu,mua,mub,ig,nu,nua,nub)
565                do mu=1,6
566                  do ig=1,npw
567                    mua=alpha(mu);mub=beta(mu)
568                    do nua=1,3
569                      do nub=1,3
570                        nu=gamma(nua,nub)
571                        d2ffnl_cart(ig,mu)=d2ffnl_cart(ig,mu)+d2ffnl_red(ig,nu)*gprimd(mua,nua)*gprimd(mub,nub)
572                      end do
573                    end do
574                  end do
575                end do
576                do mu=1,6
577                  mua=alpha(mu);mub=beta(mu)
578 !$OMP PARALLEL PRIVATE(nu,nua,nub,ig)
579 !$OMP WORKSHARE
580                  d2ffnl_tmp=zero
581 !$OMP END WORKSHARE
582 !$OMP DO
583                  do ig=1,npw
584                    do nua=1,3
585                      do nub=1,3
586                        nu=gamma(nua,nub)
587                        d2ffnl_tmp(ig)=d2ffnl_tmp(ig)+d2ffnl_cart(ig,nu)*gprimd(nua,mua)*gprimd(nub,mub)
588                      end do
589                    end do
590                  end do
591 !$OMP END DO
592 !$OMP WORKSHARE
593                  ffnl(:,4+mu,iffnl,itypat)=d2ffnl_tmp(:)
594 !$OMP END WORKSHARE
595 !$OMP END PARALLEL
596                end do
597              end if
598            end if
599 
600 !          =========================================================================
601 !          B-USE OF LEGENDRE POLYNOMIAL IN APPLICATION OF NL OPERATOR:
602 !          ffnl(K,l,n)=fnl(K)/|K|^l
603 !          ffnl_prime(K,l,n)=(fnl_prime(K)-l*fnl(K)/|K|)/|K|^(l+1)
604 !          ffnl_prime_prime(K,l,n)=(fnl_prime_prime(K)-(2*l+1)*fnl_prime(K)/|K|
605 !          +l*(l+2)*fnl(K)/|K|^2)/|K|^(l+2)
606          else if (iln>iln0) then
607 
608            if (il==0) then
609 !$OMP PARALLEL DO
610              do ig=1,npw
611                ffnl(ig,1,iffnl,itypat)=wk_ffnl1(ig)
612              end do
613            else
614 !$OMP PARALLEL DO
615              do ig=1,npw
616                ffnl(ig,1,iffnl,itypat)=wk_ffnl1(ig)*kpgnorm_inv(ig)**il
617              end do
618            end if
619            if (ider>=1) then
620 !$OMP PARALLEL DO
621              do ig=1,npw
622                ffnl(ig,2,iffnl,itypat)= (wk_ffnl2(ig)-dble(il)*wk_ffnl1(ig)*kpgnorm_inv(ig))*kpgnorm_inv(ig)**(il+1)
623              end do
624              if (ider==2) then
625 !$OMP PARALLEL DO
626                do ig=1,npw
627                  ffnl(ig,3,iffnl,itypat)= (wk_ffnl3(ig)-       &
628 &                 dble(2*il+1)*wk_ffnl2(ig)*kpgnorm_inv(ig)+   &
629 &                 dble(il*(il+2))*wk_ffnl1(ig)*kpgnorm_inv(ig)**2)*kpgnorm_inv(ig)**(il+2)
630                end do
631              end if
632            end if
633 
634          end if  ! End if - Use of Ylm or not
635 
636        else ! No NL part
637 !$OMP PARALLEL DO COLLAPSE(2)
638          do mu=1,dimffnl
639            do ig=1,npw
640              ffnl(ig,mu,iffnl,itypat)=zero
641            end do
642          end do
643 
644        end if ! End if - a nonlocal part exists
645      end if ! End if - special case : spin orbit calc. & no spin-orbit psp
646 
647      if (iln>iln0) iln0=iln
648    end do ! End do - loop over (l,m,n) values
649 
650  end do  ! End do - loop over atom types
651 
652  if (allocated(kpgc))  then
653    ABI_DEALLOCATE(kpgc)
654  end if
655  if (allocated(kpgn))  then
656    ABI_DEALLOCATE(kpgn)
657  end if
658  if (allocated(dffnl_red))  then
659    ABI_DEALLOCATE(dffnl_red)
660  end if
661  if (allocated(d2ffnl_red))  then
662    ABI_DEALLOCATE(d2ffnl_red)
663  end if
664  if (allocated(dffnl_cart))  then
665    ABI_DEALLOCATE(dffnl_cart)
666  end if
667  if (allocated(d2ffnl_cart))  then
668    ABI_DEALLOCATE(d2ffnl_cart)
669  end if
670  if (allocated(dffnl_tmp))  then
671    ABI_DEALLOCATE(dffnl_tmp)
672  end if
673  if (allocated(d2ffnl_tmp))  then
674    ABI_DEALLOCATE(d2ffnl_tmp)
675  end if
676 
677  ABI_DEALLOCATE(kpgnorm_inv)
678  ABI_DEALLOCATE(kpgnorm)
679  ABI_DEALLOCATE(wk_ffnl1)
680  ABI_DEALLOCATE(wk_ffnl2)
681  ABI_DEALLOCATE(wk_ffnl3)
682  ABI_DEALLOCATE(wk_ffspl)
683 
684  call timab(16,2,tsec)
685 
686  DBG_EXIT("COLL")
687 
688 end subroutine mkffnl