TABLE OF CONTENTS
ABINIT/m_mkffnl [ 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 ]
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 ]
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