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-2022 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 .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_mkffnl
24 
25  use defs_basis
26  use m_abicore
27  use m_errors
28  use m_splines
29  use m_xmpi
30 
31  use m_time,     only : timab
32  use m_kg,       only : mkkin
33  use m_sort,     only : sort_dp
34  use defs_datatypes, only : pseudopotential_type
35  use m_crystal,  only : crystal_t
36 
37  implicit none
38 
39  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)
 [comm]=MPI communicator. Default: xmpi_comm_self.

OUTPUT

  ffnl(npw,dimffnl,lmnmax,ntypat)=described below
 [request]=Used in conjunction with [comm] to perform non-blocking xmpi_isum_ip. Client code must
  wait on request before using ffnl. If not present, blocking API is used.

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.

SOURCE

238 subroutine mkffnl(dimekb, dimffnl, ekb, ffnl, ffspl, gmet, gprimd, ider, idir, indlmn, &
239                   kg, kpg, kpt, lmnmax, lnmax, mpsang, mqgrid, nkpg, npw, ntypat, pspso, &
240                   qgrid, rmet, usepaw, useylm, ylm, ylm_gr, &
241                   comm, request) ! optional
242 
243 !Arguments ------------------------------------
244 !scalars
245  integer,intent(in) :: dimekb,dimffnl,ider,idir,lmnmax,lnmax,mpsang,mqgrid,nkpg
246  integer,intent(in) :: npw,ntypat,usepaw,useylm
247  integer,optional,intent(in) :: comm
248  integer ABI_ASYNC, optional,intent(out):: request
249 !arrays
250  integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg(3,npw),pspso(ntypat)
251  real(dp),intent(in) :: ekb(dimekb,ntypat*(1-usepaw))
252  real(dp),intent(in) :: ffspl(mqgrid,2,lnmax,ntypat),gmet(3,3),gprimd(3,3)
253  real(dp),intent(in) :: kpg(npw,nkpg),kpt(3),qgrid(mqgrid),rmet(3,3)
254  real(dp),intent(in) :: ylm(:,:),ylm_gr(:,:,:)
255  real(dp),intent(out) :: ffnl(npw,dimffnl,lmnmax,ntypat)
256  ! MG: Should be ABI_ASYNC due to optional non-Blocking API but NAG complains
257  ! Error: m_d2frnl.F90, line 600: Array section FFNL_STR(:,:,:,:,MU) supplied for dummy FFNL (no. 4) of MKFFNL,
258  ! the dummy is ASYNCHRONOUS but not assumed-shape
259  ! so we declare request as ASYNCHRONOUS
260 
261 !Local variables-------------------------------
262 !scalars
263  integer :: ider_tmp,iffnl,ig,ig0,il,ilm,ilmn,iln,iln0,im,itypat,mu,mua,mub,nlmn,nu,nua,nub
264  integer :: nprocs, my_rank, cnt, ierr
265  real(dp),parameter :: renorm_factor=0.5d0/pi**2,tol_norm=tol10
266  real(dp) :: ecut,ecutsm,effmass_free,fact,kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3,rmetab,yp1
267  logical :: testnl=.false.
268  character(len=500) :: msg
269 !arrays
270  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
271  integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
272  real(dp) :: rprimd(3,3),tsec(2)
273  real(dp),allocatable :: dffnl_cart(:,:),dffnl_red(:,:),dffnl_tmp(:)
274  real(dp),allocatable :: d2ffnl_cart(:,:),d2ffnl_red(:,:),d2ffnl_tmp(:)
275  real(dp),allocatable :: kpgc(:,:),kpgn(:,:),kpgnorm(:),kpgnorm_inv(:)
276  real(dp),allocatable :: wk_ffnl1(:),wk_ffnl2(:),wk_ffnl3(:),wk_ffspl(:,:)
277 
278 ! *************************************************************************
279 
280  ! Keep track of time spent in mkffnl
281  call timab(16, 1, tsec)
282 
283  nprocs = 1; my_rank = 0
284  if (present(comm)) then
285    nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
286  end if
287 
288  ! Compatibility tests
289  if (mpsang>4) then
290    write(msg,'(a,i0,a,a)')&
291    'Called with mpsang > 4, =',mpsang,ch10,&
292    'This subroutine will not accept lmax+1 > 4.'
293    ABI_BUG(msg)
294  end if
295  if (idir<-7.or.idir>4) then
296    ABI_BUG('Called with idir<-6 or idir>4 !')
297  end if
298  if (useylm==0) then
299    iffnl=1+ider
300  else
301    iffnl=1
302    if (ider>=1) then
303      if (idir==0) iffnl=iffnl+3
304      if (idir/=0) iffnl=iffnl+1
305      if (idir==4) iffnl=iffnl+2
306      if (idir==-7) iffnl=iffnl+5
307    end if
308    if (ider==2) then
309      if (idir==0) iffnl=iffnl+6
310      if (idir==4) iffnl=iffnl+6
311    end if
312  end if
313  if (iffnl/=dimffnl) then
314    write(msg,'(2(a,i1),a,i2)') 'Incompatibility between ider, idir and dimffnl : ider = ',ider,&
315                                ' idir = ',idir,' dimffnl = ',dimffnl
316    ABI_BUG(msg)
317  end if
318  if (useylm==1) then
319    ABI_CHECK(size(ylm,1)==npw,'BUG: wrong ylm size (1)')
320    ABI_CHECK(size(ylm,2)==mpsang**2,'BUG: wrong ylm size (2)')
321    if(ider>0)then
322      ABI_CHECK(size(ylm_gr,1)==npw,'BUG: wrong ylm_gr size (1)')
323      ABI_CHECK(size(ylm_gr,2)>=3+6*(ider/2),'BUG: wrong ylm_gr size (2)')
324      ABI_CHECK(size(ylm_gr,3)==mpsang**2,'BUG: wrong ylm_gr size (3)')
325    end if
326  end if
327 
328  ! Get (k+G) and |k+G|
329  ABI_MALLOC(kpgnorm,(npw))
330  ABI_MALLOC(kpgnorm_inv,(npw))
331 
332  ig0=-1 ! index of |k+g|=0 vector
333 
334  if (useylm==1) then
335    ABI_MALLOC(kpgc,(npw,3))
336    if (ider>=1) then
337      ABI_MALLOC(kpgn,(npw,3))
338    end if
339    if (nkpg<3) then
340 !$OMP PARALLEL DO PRIVATE(ig,kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3)
341      do ig=1,npw
342        kpg1=kpt(1)+dble(kg(1,ig))
343        kpg2=kpt(2)+dble(kg(2,ig))
344        kpg3=kpt(3)+dble(kg(3,ig))
345        kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3)
346        kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3)
347        kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3)
348        kpgc(ig,1)=kpgc1
349        kpgc(ig,2)=kpgc2
350        kpgc(ig,3)=kpgc3
351        kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
352        if (kpgnorm(ig)<=tol_norm) ig0=ig
353        if (ider>=1) then
354          kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
355          kpgn(ig,1)=kpg1*kpgnorm_inv(ig)
356          kpgn(ig,2)=kpg2*kpgnorm_inv(ig)
357          kpgn(ig,3)=kpg3*kpgnorm_inv(ig)
358        end if
359      end do
360    else
361 !$OMP PARALLEL DO PRIVATE(ig,kpgc1,kpgc2,kpgc3)
362      do ig=1,npw
363        kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
364        kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
365        kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
366        kpgc(ig,1)=kpgc1
367        kpgc(ig,2)=kpgc2
368        kpgc(ig,3)=kpgc3
369        kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
370        if (kpgnorm(ig)<=tol_norm) ig0=ig
371        if (ider>=1) then
372          kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
373          kpgn(ig,1:3)=kpg(ig,1:3)*kpgnorm_inv(ig)
374        end if
375      end do
376    end if
377  else
378    if (nkpg<3) then
379      ecut=huge(0.0d0)*0.1d0;ecutsm=zero;effmass_free=one
380      ! Note that with ecutsm=0, the right kinetic energy is computed
381      call mkkin(ecut,ecutsm,effmass_free,gmet,kg,kpgnorm,kpt,npw,0,0)
382 !$OMP PARALLEL DO
383      do ig=1,npw
384        kpgnorm(ig)=sqrt(renorm_factor*kpgnorm(ig))
385        kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
386        if (kpgnorm(ig)<=tol_norm) ig0=ig
387      end do
388    else
389 !$OMP PARALLEL DO PRIVATE(ig,kpgc1,kpgc2,kpgc3)
390      do ig=1,npw
391        kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
392        kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
393        kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
394        kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
395        kpgnorm_inv(ig)=1.d0/max(kpgnorm(ig),tol_norm)
396        if (kpgnorm(ig)<=tol_norm) ig0=ig
397      end do
398    end if
399  end if
400 
401  ! Need rprimd in some cases
402  if (ider>=1.and.useylm==1.and.ig0>0) then
403    do mu=1,3
404      do nu=1,3
405        rprimd(mu,nu)=gprimd(mu,1)*rmet(1,nu)+gprimd(mu,2)*rmet(2,nu)+gprimd(mu,3)*rmet(3,nu)
406      end do
407    end do
408  end if
409 
410  ! Allocate several temporary arrays
411  ABI_MALLOC(wk_ffnl1,(npw))
412  ABI_MALLOC(wk_ffnl2,(npw))
413  ABI_MALLOC(wk_ffnl3,(npw))
414  ABI_MALLOC(wk_ffspl,(mqgrid,2))
415 
416  if (ider>=1.and.useylm==1) then
417    ABI_MALLOC(dffnl_red,(npw,3))
418    if (idir/=0) then
419      ABI_MALLOC(dffnl_cart,(npw,3))
420    end if
421    if (idir>0) then
422      ABI_MALLOC(dffnl_tmp,(npw))
423    end if
424  end if
425  if (ider>=2 .and. useylm==1) then
426    ABI_MALLOC(d2ffnl_red,(npw,6))
427    if (idir==4) then
428      ABI_MALLOC(d2ffnl_cart,(npw,6))
429      ABI_MALLOC(d2ffnl_tmp,(npw))
430    end if
431  end if
432 
433  ! Loop over types of atoms
434  ffnl = zero; cnt = 0
435  do itypat=1,ntypat
436 
437    ! Loop over (l,m,n) values
438    iln0=0; nlmn=count(indlmn(3,:,itypat)>0)
439 
440    do ilmn=1,nlmn
441      il=indlmn(1,ilmn,itypat)
442      im=indlmn(2,ilmn,itypat)
443      ilm =indlmn(4,ilmn,itypat)
444      iln =indlmn(5,ilmn,itypat)
445      iffnl=ilmn;if (useylm==0) iffnl=iln
446 
447      ! Special case: we enter the loop in case of spin-orbit calculation
448      ! even if the psp has no spin-orbit component.
449      if (indlmn(6,ilmn,itypat) ==1 .or. pspso(itypat) /=0) then
450 
451        ! Compute FFNL only if ekb>0 or paw
452        if (usepaw==1) testnl=.true.
453        if (usepaw==0) testnl=(abs(ekb(iln,itypat))>tol_norm)
454 
455        if (testnl) then
456          cnt = cnt + 1
457          if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism (optional)
458          !
459          ! Store form factors (from ffspl)
460          ! -------------------------------
461          ! MG: This part is an hotspot of in the EPH code due to the large number of k-points used
462          ! To improve memory locality, I tried to:
463          !
464          !      1) call a new version of splfit that operates on wk_ffspl with shape: (2,mqgrid)
465          !      2) pass a sorted kpgnorm array and then rearrange the output spline
466          !
467          ! but I didn't manage to make it significantly faster.
468          ! For the time being, we rely on MPI-parallelism via the optional MPI communicator.
469 
470          if (iln > iln0) then
471            wk_ffspl(:,:)=ffspl(:,:,iln,itypat)
472            ider_tmp = min(ider, 1)
473            call splfit(qgrid,wk_ffnl2,wk_ffspl,ider_tmp,kpgnorm,wk_ffnl1,mqgrid,npw)
474            if (ider == 2) then
475              call splfit(qgrid,wk_ffnl3,wk_ffspl,ider,kpgnorm,wk_ffnl1,mqgrid,npw)
476            end if
477          end if
478 
479          ! Store FFNL and FFNL derivatives
480          ! -------------------------------
481 
482          ! =========================================================================
483          ! A-USE OF SPHER. HARMONICS IN APPLICATION OF NL OPERATOR:
484          ! ffnl(K,l,m,n)=fnl(K).Ylm(K)
485          ! --if (idir==0)
486          ! ffnl_prime(K,1:3,l,m,n)=3 reduced coordinates of d(ffnl)/dK^cart
487          ! =fnl_prime(K).Ylm(K).K^red_i/|K|+fnl(K).(dYlm/dK^cart)^red_i
488          ! --if (0<idir<4)
489          ! ffnl_prime(K,l,m,n)=cart. coordinate idir of d(ffnl)/dK^red
490          ! --if (idir==4)
491          ! ffnl_prime(K,l,m,n)=3 cart. coordinates of d(ffnl)/dK^red
492          ! --if (-7<=idir<0) - |idir|=(mu,nu) (1->11,2->22,3->33,4->32,5->31,6->21)
493          ! ffnl_prime(K,l,m,n)=1/2 [d(ffnl)/dK^cart_mu K^cart_nu + d(ffnl)/dK^cart_nu K^cart_mu]
494          ! ffnl_prime_prime(K,l,m,n)=6 reduced coordinates of d2(ffnl)/dK^cart.dK^cart
495 
496          if (useylm==1) then
497 !$OMP PARALLEL DO
498            do ig=1,npw
499              ffnl(ig,1,iffnl,itypat)=ylm(ig,ilm)*wk_ffnl1(ig)
500            end do
501 
502            if (ider>=1) then
503 !$OMP PARALLEL DO COLLAPSE(2)
504              do mu=1,3
505                do ig=1,npw
506                  dffnl_red(ig,mu)=ylm(ig,ilm)*wk_ffnl2(ig)*kpgn(ig,mu)+ylm_gr(ig,mu,ilm)*wk_ffnl1(ig)
507                end do
508              end do
509              ! Special cases |k+g|=0
510              if (ig0>0) then
511                do mu=1,3
512                  dffnl_red(ig0,mu)=zero
513                  if (il==1) then
514                    !Retrieve 1st-deriv. of ffnl at q=zero according to spline routine
515                    yp1=(wk_ffspl(2,1)-wk_ffspl(1,1))/qgrid(2)-sixth*qgrid(2)*(two*wk_ffspl(1,2)+wk_ffspl(2,2))
516                    fact=yp1*sqrt(three/four_pi)
517                    if (im==-1) dffnl_red(ig0,mu)=fact*rprimd(2,mu)
518                    if (im== 0) dffnl_red(ig0,mu)=fact*rprimd(3,mu)
519                    if (im==+1) dffnl_red(ig0,mu)=fact*rprimd(1,mu)
520                  end if
521                end do
522              end if
523              if (idir==0) then
524 !$OMP PARALLEL DO COLLAPSE(2)
525                do mu=1,3
526                  do ig=1,npw
527                    ffnl(ig,1+mu,iffnl,itypat)=dffnl_red(ig,mu)
528                  end do
529                end do
530              else
531                dffnl_cart=zero
532 !$OMP PARALLEL DO COLLAPSE(2)
533                do mu=1,3
534                  do ig=1,npw
535                    do nu=1,3
536                      dffnl_cart(ig,mu)=dffnl_cart(ig,mu)+dffnl_red(ig,nu)*gprimd(mu,nu)
537                    end do
538                  end do
539                end do
540                if (idir>=1.and.idir<=3) then
541                  dffnl_tmp=zero
542 !$OMP PARALLEL PRIVATE(nu,ig)
543 !$OMP DO
544                  do ig=1,npw
545                    do nu=1,3
546                      dffnl_tmp(ig)=dffnl_tmp(ig) + dffnl_cart(ig,nu)*gprimd(nu,idir)
547                    end do
548                  end do
549 !$OMP END DO
550 !$OMP WORKSHARE
551                  ffnl(:,2,iffnl,itypat)=dffnl_tmp(:)
552 !$OMP END WORKSHARE
553 !$OMP END PARALLEL
554                else if (idir==4) then
555                  do mu=1,3
556 !$OMP PARALLEL PRIVATE(nu,ig)
557 !$OMP WORKSHARE
558                    dffnl_tmp=zero
559 !$OMP END WORKSHARE
560 !$OMP DO
561                    do ig=1,npw
562                      do nu=1,3
563                        dffnl_tmp(ig)=dffnl_tmp(ig) + dffnl_cart(ig,nu)*gprimd(nu,mu)
564                      end do
565                    end do
566 !$OMP END DO
567 !$OMP WORKSHARE
568                    ffnl(:,1+mu,iffnl,itypat)=dffnl_tmp(:)
569 !$OMP END WORKSHARE
570 !$OMP END PARALLEL
571                  end do
572                else if (idir/=-7) then
573                  mu=abs(idir);mua=alpha(mu);mub=beta(mu)
574 !$OMP PARALLEL DO
575                  do ig=1,npw
576                    ffnl(ig,2,iffnl,itypat)=0.5d0* (dffnl_cart(ig,mua)*kpgc(ig,mub) + dffnl_cart(ig,mub)*kpgc(ig,mua))
577                  end do
578                else if (idir==-7) then
579 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(mua, mub)
580                  do mu=1,6
581                    do ig=1,npw
582                      mua=alpha(mu);mub=beta(mu)
583                      ffnl(ig,1+mu,iffnl,itypat)=0.5d0 * (dffnl_cart(ig,mua)*kpgc(ig,mub) + dffnl_cart(ig,mub)*kpgc(ig,mua))
584                    end do
585                  end do
586                end if
587              end if
588            end if
589 
590            if (ider==2) then
591              do mu=1,6
592                mua=alpha(mu);mub=beta(mu)
593                rmetab=rmet(mua,mub)
594 !$OMP PARALLEL DO
595                do ig=1,npw
596                  d2ffnl_red(ig,mu)= &
597                  ylm_gr(ig,3+mu,ilm)*wk_ffnl1(ig) &
598                  + (rmetab-kpgn(ig,mua)*kpgn(ig,mub))*ylm(ig,ilm)*wk_ffnl2(ig)*kpgnorm_inv(ig) &
599                  + ylm(ig,ilm)*kpgn(ig,mua)*kpgn(ig,mub)*wk_ffnl3(ig) &
600                  + (ylm_gr(ig,mua,ilm)*kpgn(ig,mub)+ylm_gr(ig,mub,ilm)*kpgn(ig,mua))*wk_ffnl2(ig)
601                end do
602                ! Special cases |k+g|=0
603                if (ig0>0) then
604                  d2ffnl_red(ig0,mu)=zero
605                  if (il==0) then
606                    d2ffnl_red(ig0,mu)=wk_ffspl(1,2)*rmetab/sqrt(four_pi)
607                  end if
608                  if (il==2) then
609                    fact=wk_ffspl(1,2)*quarter*sqrt(15._dp/pi)
610                    if (im==-2) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(2,mub)+rprimd(2,mua)*rprimd(1,mub))
611                    if (im==-1) d2ffnl_red(ig0,mu)=fact*(rprimd(2,mua)*rprimd(3,mub)+rprimd(3,mua)*rprimd(2,mub))
612                    if (im==+1) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(3,mub)+rprimd(3,mua)*rprimd(1,mub))
613                    if (im==+2) d2ffnl_red(ig0,mu)=fact*(rprimd(1,mua)*rprimd(1,mub)-rprimd(2,mua)*rprimd(2,mub))
614                    if (im== 0) d2ffnl_red(ig0,mu)=(fact/sqrt3)*(two*rprimd(3,mua)*rprimd(3,mub) &
615                                                   -rprimd(1,mua)*rprimd(1,mub)-rprimd(2,mua)*rprimd(2,mub))
616                  end if
617                end if
618              end do
619              if (idir==0) then
620 !$OMP PARALLEL DO COLLAPSE(2)
621                do mu=1,6
622                  do ig=1,npw
623                    ffnl(ig,4+mu,iffnl,itypat)=d2ffnl_red(ig,mu)
624                  end do
625                end do
626              else if (idir==4) then
627                d2ffnl_cart=zero
628 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(mu,mua,mub,ig,nu,nua,nub)
629                do mu=1,6
630                  do ig=1,npw
631                    mua=alpha(mu);mub=beta(mu)
632                    do nua=1,3
633                      do nub=1,3
634                        nu=gamma(nua,nub)
635                        d2ffnl_cart(ig,mu)=d2ffnl_cart(ig,mu)+d2ffnl_red(ig,nu)*gprimd(mua,nua)*gprimd(mub,nub)
636                      end do
637                    end do
638                  end do
639                end do
640                do mu=1,6
641                  mua=alpha(mu);mub=beta(mu)
642 !$OMP PARALLEL PRIVATE(nu,nua,nub,ig)
643 !$OMP WORKSHARE
644                  d2ffnl_tmp=zero
645 !$OMP END WORKSHARE
646 !$OMP DO
647                  do ig=1,npw
648                    do nua=1,3
649                      do nub=1,3
650                        nu=gamma(nua,nub)
651                        d2ffnl_tmp(ig)=d2ffnl_tmp(ig)+d2ffnl_cart(ig,nu)*gprimd(nua,mua)*gprimd(nub,mub)
652                      end do
653                    end do
654                  end do
655 !$OMP END DO
656 !$OMP WORKSHARE
657                  ffnl(:,4+mu,iffnl,itypat)=d2ffnl_tmp(:)
658 !$OMP END WORKSHARE
659 !$OMP END PARALLEL
660                end do
661              end if
662            end if
663 
664            ! =========================================================================
665            ! B-USE OF LEGENDRE POLYNOMIAL IN APPLICATION OF NL OPERATOR:
666            ! ffnl(K,l,n)=fnl(K)/|K|^l
667            ! ffnl_prime(K,l,n)=(fnl_prime(K)-l*fnl(K)/|K|)/|K|^(l+1)
668            ! ffnl_prime_prime(K,l,n)=(fnl_prime_prime(K)-(2*l+1)*fnl_prime(K)/|K|
669            ! +l*(l+2)*fnl(K)/|K|^2)/|K|^(l+2)
670          else if (iln>iln0) then
671 
672            if (il==0) then
673 !$OMP PARALLEL DO
674              do ig=1,npw
675                ffnl(ig,1,iffnl,itypat)=wk_ffnl1(ig)
676              end do
677            else
678 !$OMP PARALLEL DO
679              do ig=1,npw
680                ffnl(ig,1,iffnl,itypat)=wk_ffnl1(ig)*kpgnorm_inv(ig)**il
681              end do
682            end if
683            if (ider>=1) then
684 !$OMP PARALLEL DO
685              do ig=1,npw
686                ffnl(ig,2,iffnl,itypat)= (wk_ffnl2(ig)-dble(il)*wk_ffnl1(ig)*kpgnorm_inv(ig))*kpgnorm_inv(ig)**(il+1)
687              end do
688              if (ider==2) then
689 !$OMP PARALLEL DO
690                do ig=1,npw
691                  ffnl(ig,3,iffnl,itypat)= (wk_ffnl3(ig)-       &
692                    dble(2*il+1)*wk_ffnl2(ig)*kpgnorm_inv(ig)+   &
693                    dble(il*(il+2))*wk_ffnl1(ig)*kpgnorm_inv(ig)**2)*kpgnorm_inv(ig)**(il+2)
694                end do
695              end if
696            end if
697 
698          end if  ! Use of Ylm or not
699 
700        else
701          ! No NL part
702 !$OMP PARALLEL DO COLLAPSE(2)
703          do mu=1,dimffnl
704            do ig=1,npw
705              ffnl(ig,mu,iffnl,itypat)=zero
706            end do
707          end do
708 
709        end if ! testnl (a nonlocal part exists)
710      end if ! special case: spin orbit calc. & no spin-orbit psp
711 
712      if (iln > iln0) iln0 = iln
713 
714    end do ! loop over (l,m,n) values
715  end do ! loop over atom types
716 
717  ABI_FREE(kpgnorm_inv)
718  ABI_FREE(kpgnorm)
719  ABI_FREE(wk_ffnl1)
720  ABI_FREE(wk_ffnl2)
721  ABI_FREE(wk_ffnl3)
722  ABI_FREE(wk_ffspl)
723 
724  ! Optional deallocations.
725  ABI_SFREE(kpgc)
726  ABI_SFREE(kpgn)
727  ABI_SFREE(dffnl_red)
728  ABI_SFREE(d2ffnl_red)
729  ABI_SFREE(dffnl_cart)
730  ABI_SFREE(d2ffnl_cart)
731  ABI_SFREE(dffnl_tmp)
732  ABI_SFREE(d2ffnl_tmp)
733 
734  if (nprocs > 1) then
735    ! Blocking/non-blocking depending on the presence of request.
736    if (present(request)) then
737      call xmpi_isum_ip(ffnl, comm, request, ierr)
738    else
739      call xmpi_sum(ffnl, comm, ierr)
740    end if
741  else
742    if (present(request)) request = xmpi_request_null
743  end if
744 
745  call timab(16, 2, tsec)
746 
747 end subroutine mkffnl

ABINIT/mkffnl_objs [ Functions ]

[ Top ] [ Functions ]

NAME

 mkffnl_objs

FUNCTION

  Simplified wrapper around mkffnl in which input parameters are passed via crystal_t and pseudopotential_type.

INPUTS

  See mkffnl

OUTPUT

  ffnl(npw,dimffnl,lmnmax,ntypat)=described below
 [request]=Used in conjunction with [comm] to perform non-blocking xmpi_isum_ip. Client code must
  wait on request before using ffnl. If not present, blocking API is used.

SOURCE

 66 subroutine mkffnl_objs(cryst, psps, dimffnl, ffnl, ider, idir, kg, kpg, kpt, nkpg, npw, ylm, ylm_gr, &
 67                        comm, request) ! optional
 68 
 69 !Arguments ------------------------------------
 70 !scalars
 71  type(crystal_t),intent(in) :: cryst
 72  type(pseudopotential_type),intent(in) :: psps
 73  integer,intent(in) :: dimffnl, ider, idir, npw, nkpg
 74  integer,optional,intent(in) :: comm
 75  integer ABI_ASYNC, optional,intent(out):: request
 76 !arrays
 77  integer,intent(in) :: kg(3,npw)
 78  real(dp),intent(in) :: kpg(npw, nkpg), kpt(3)
 79  real(dp),intent(in) :: ylm(:,:), ylm_gr(:,:,:)
 80  real(dp),intent(out) :: ffnl(npw, dimffnl, psps%lmnmax, psps%ntypat)
 81 !
 82 !!Local variables-------------------------------
 83  integer :: my_comm
 84 
 85 ! *************************************************************************
 86  my_comm = xmpi_comm_self; if (present(comm)) my_comm = comm
 87 
 88  if (present(request)) then
 89     call mkffnl(psps%dimekb, dimffnl, psps%ekb, ffnl, psps%ffspl, &
 90                 cryst%gmet, cryst%gprimd, ider, idir, psps%indlmn, kg, kpg, kpt, psps%lmnmax, &
 91                 psps%lnmax, psps%mpsang, psps%mqgrid_ff, nkpg, npw, psps%ntypat, &
 92                 psps%pspso, psps%qgrid_ff, cryst%rmet, psps%usepaw, psps%useylm, ylm, ylm_gr, &
 93                 comm=comm, request=request)
 94 
 95  else
 96     call mkffnl(psps%dimekb, dimffnl, psps%ekb, ffnl, psps%ffspl, &
 97                 cryst%gmet, cryst%gprimd, ider, idir, psps%indlmn, kg, kpg, kpt, psps%lmnmax, &
 98                 psps%lnmax, psps%mpsang, psps%mqgrid_ff, nkpg, npw, psps%ntypat, &
 99                 psps%pspso, psps%qgrid_ff, cryst%rmet, psps%usepaw, psps%useylm, ylm, ylm_gr, &
100                 comm=comm)
101  end if
102 
103 end subroutine mkffnl_objs