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 (C) 1998-2018 ABINIT group (DCA, XG, GMR, MT, DRH)
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:
- 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
$$\textrm{f}_l(q)=\frac{1}{dvrms} \int_0^\infty [j_l(2 \pi r q) u_l(r) dV(r) r dr]$$
where u_l(r)=reference state wavefunction, dV(r)=nonlocal psp
correction, j_l(arg)=spherical Bessel function for angular momentum l,
and
$$\textrm{dvrms} = \int_0^\infty [(u_l(r) dV(r))^2 dr])^{1/2}$$
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