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