TABLE OF CONTENTS
- ABINIT/m_vkbr
- m_vkbr/add_vnlr_commutator
- m_vkbr/calc_vkb
- m_vkbr/ccgradvnl_ylm
- m_vkbr/nc_ihr_comm
- m_vkbr/vkbr_free_0D
- m_vkbr/vkbr_free_1D
- m_vkbr/vkbr_init
- m_vkbr/vkbr_t
ABINIT/m_vkbr [ Modules ]
NAME
m_vkbr
FUNCTION
This module provides objects and methods used to calculate the matrix elements
of the commutator [H,r] needed for the correct treatment of the optical limit q-->0
in the matrix elements <k-q,b1|e^{-iqr}|k,b2> when non-local pseudopotentials are used.
NOTES
This module is deprecated. Use ddkop_t in m_ddk.F90
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, FB) 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
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 MODULE m_vkbr 28 29 use defs_basis 30 use m_hide_blas 31 use m_errors 32 use m_abicore 33 34 use defs_datatypes, only : pseudopotential_type 35 use m_gwdefs, only : czero_gw 36 use m_fstrings, only : sjoin, itoa 37 use m_paw_sphharm, only : ylmc, ylmcd 38 use m_geometry, only : normv 39 use m_crystal, only : crystal_t 40 use m_kg, only : mkkin 41 use m_mkffnl, only : mkffnl 42 43 implicit none 44 45 private
m_vkbr/add_vnlr_commutator [ Functions ]
[ Top ] [ m_vkbr ] [ Functions ]
NAME
add_vnlr_commutator
FUNCTION
Calculate the matrix elements of the dipole operator <phi1|r|phi2>. For norm conserving potentials the commutator [Vnl,r] is included according to inclvkb.
INPUTS
vkbr<vkbr_t> cryst<crystal_t>=Datatype gathering info on the crystal structure. psps<pseudopotential_type>Structure gathering info on the pseudopotentials. npw=Number of G for wavefunctions. nspinor=Number of spinorial components. ug1(npw*nspinor)=Left wavefunction. ug2(npw*nspinor)=Right wavefunction
SIDE EFFECTS
rhotwx(3,nspinor**2)= Updated. Matrix elements in reduced coordinates, see NOTES below.
NOTES
1) <k b1|e^{-iq.r}|k b2> = \delta_{b1 b2} -iq <k b1|r|k b2> = \delta_{b1 b2} -iq ( <k b1| [H,r] |k b2> / (e1-e2) ).
This routine calculates the matrix elements of ir*(e1-e2)
Remember that [H,r] = -\nabla + [V_nl,r]
2) The Fourier transform of a two-point real function f(r1,r2) satisfies:
a) f_{\Gamma}(G1,G2) = f_{\Gamma}(-G1,-G2)^*
b) f_{G0/2} (G1,G2) = f_{G0/2}(-G1-G0,-G2-G0)^*
TODO
*) Spinorial case is not implemented.
SOURCE
289 subroutine add_vnlr_commutator(vkbr,cryst,psps,npw,nspinor,ug1,ug2,rhotwx) 290 291 !Arguments ------------------------------------ 292 !scalars 293 integer,intent(in) :: npw,nspinor 294 type(vkbr_t),intent(in) :: vkbr 295 type(crystal_t),intent(in) :: cryst 296 type(pseudopotential_type),intent(in) :: psps 297 !arrays 298 complex(gwpc),target,intent(in) :: ug1(npw*nspinor),ug2(npw*nspinor) 299 complex(gwpc),intent(inout) :: rhotwx(3,nspinor**2) 300 301 !Local variables ------------------------------ 302 !scalars 303 integer :: iat,ig,ilm,itypat,nlmn,ilmn,iln0,iln,il,in,im 304 complex(gwpc) :: cta1,cta4 305 !arrays 306 complex(gwpc) :: dum(3),cta2(3),cta3(3),gamma_term(3) 307 308 !************************************************************************ 309 310 ABI_CHECK(nspinor == 1, "inclvkb > 0 with nspinor == 2 is not coded") 311 312 ! Adding term i <c,k|[Vnl,r]|v,k> === 313 select case (vkbr%inclvkb) 314 case (2) 315 ! Complex spherical harmonics (much faster!). 316 dum=czero_gw; gamma_term=czero 317 318 do iat=1,vkbr%natom 319 itypat = cryst%typat(iat) 320 nlmn = count(psps%indlmn(3,:,itypat) > 0) 321 iln0 = 0 322 do ilmn=1,nlmn 323 il = 1 + psps%indlmn(1,ilmn,itypat) 324 in = psps%indlmn(3,ilmn,itypat) 325 iln = psps%indlmn(5,ilmn,itypat) 326 if (iln <= iln0) cycle 327 iln0 = iln 328 !if (indlmn(6,ilmn,itypat) /= 1 .or. vkbsign(iln,itypat) == zero) cycle 329 !in = 1 330 do im=1,2*(il-1)+1 331 ! Index of im and il 332 ilm = im + (il-1)*(il-1) 333 cta1 = czero_gw; cta2(:) = czero_gw 334 cta4 = czero_gw; cta3(:) = czero_gw 335 do ig=1,npw 336 ! Here we take advantage of the property Y_(l-m)= (-i)^m Y_lm^*. 337 cta1 = cta1 + ug1(ig) * vkbr%fnl (ig,ilm,in,iat) 338 cta2(:)= cta2(:) + ug2(ig) * vkbr%fnld(:,ig,ilm,in,iat) 339 cta3(:)= cta3(:) + ug1(ig) * vkbr%fnld(:,ig,ilm,in,iat) 340 cta4 = cta4 + ug2(ig) * vkbr%fnl (ig,ilm,in,iat) 341 if (ig==1) gamma_term = gamma_term + CONJG(cta1)*cta2(:) +CONJG(cta3(:))*cta4 342 end do 343 dum(:)= dum(:) + CONJG(cta1)*cta2(:) + CONJG(cta3(:))*cta4 344 end do 345 346 end do 347 end do 348 349 if (vkbr%istwfk>1) then 350 dum = two * j_dpc * AIMAG(dum); if (vkbr%istwfk==2) dum = dum - j_dpc * AIMAG(gamma_term) 351 end if 352 rhotwx(:,1) = rhotwx(:,1) + dum(:) 353 354 case default 355 ABI_ERROR(sjoin("Wrong inclvkb:", itoa(vkbr%inclvkb))) 356 end select 357 358 end subroutine add_vnlr_commutator
m_vkbr/calc_vkb [ Functions ]
[ Top ] [ m_vkbr ] [ Functions ]
NAME
calc_vkb
FUNCTION
This routine calculates the Kleynman-Bylander form factors and its derivatives needed for the evaluation of the matrix elements of the dipole operator <phi1|r|phi2>.
INPUTS
cryst<crystal_t>=Crystalline structure psps<pseudopotential_type>=Structured datatype gathering information on the pseudopotentials. kpoint(3)=The k-point in reduced coordinates. npw_k=Number of plane waves for this k-point. kg_k(3,npw_k)=Reduced coordinates of the G-vectors. rprimd(3,3)=dimensional primitive translations for real space (bohr)
OUTPUT
vkb (npw_k, %lnmax, %ntypat)=KB form factors. vkbd(npw_k, %lnmax, %ntypat)=KB form factor derivatives. vkbsign(%lnmax, %ntypat) =KS dyadic sign.
TODO
SOC not implemented.
SOURCE
389 subroutine calc_vkb(cryst,psps,kpoint,npw_k,mpw,kg_k,vkbsign,vkb,vkbd) 390 391 !Arguments ------------------------------------ 392 !scalars 393 integer,intent(in) :: npw_k, mpw 394 type(crystal_t),intent(in) :: cryst 395 type(pseudopotential_type),intent(in) :: psps 396 !arrays 397 integer,intent(in) :: kg_k(3,npw_k) 398 real(dp),intent(in) :: kpoint(3) 399 real(dp),intent(out) :: vkb (mpw,psps%lnmax,psps%ntypat) 400 real(dp),intent(out) :: vkbd(mpw,psps%lnmax,psps%ntypat) 401 real(dp),intent(out) :: vkbsign(psps%lnmax,psps%ntypat) 402 403 !Local variables ------------------------------ 404 !scalars 405 integer :: dimffnl,ider,idir,itypat,nkpg,in,il,ilmn,ig,iln,iln0,nlmn 406 real(dp) :: effmass_free,ecutsm,ecut 407 !arrays 408 real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:),ylm_gr(:,:,:),ylm_k(:,:) 409 410 ! ************************************************************************* 411 412 DBG_ENTER("COLL") 413 ABI_CHECK(psps%usepaw==0, "You should not be here!") 414 ABI_CHECK(psps%useylm==0, "useylm/=0 not considered!") 415 416 ! Compute KB dyadic sign. 417 vkbsign=zero 418 do itypat=1,psps%ntypat 419 iln0 = 0 420 nlmn = count(psps%indlmn(3,:,itypat) > 0) 421 do ilmn=1,nlmn 422 iln = psps%indlmn(5,ilmn,itypat) 423 if (iln <= iln0) cycle 424 iln0 = iln 425 if (abs(psps%ekb(iln,itypat)) > 1.0d-10) vkbsign(iln,itypat) = dsign(one, psps%ekb(iln,itypat)) 426 end do 427 end do 428 429 ! Allocate KB form factor and derivative wrt k+G 430 ! Here we do not use correct ordering for dimensions 431 idir=0; nkpg=0; ider=1; dimffnl=2 ! To retrieve the first derivative. 432 433 ! Quantities used only if useylm==1 434 ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2*Psps%useylm)) 435 ABI_MALLOC(ylm_gr, (npw_k, 3+6*(ider/2),psps%mpsang**2*Psps%useylm)) 436 ABI_MALLOC(kpg_dum, (npw_k, nkpg)) 437 ABI_MALLOC(ffnl, (npw_k,dimffnl, psps%lmnmax, psps%ntypat)) 438 439 call mkffnl(psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,cryst%gmet,cryst%gprimd,ider,idir,Psps%indlmn,& 440 kg_k,kpg_dum,kpoint,psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,& 441 psps%ntypat,Psps%pspso,Psps%qgrid_ff,cryst%rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr) 442 443 ABI_FREE(ylm_k) 444 ABI_FREE(ylm_gr) 445 ABI_FREE(kpg_dum) 446 447 ABI_MALLOC(modkplusg, (npw_k)) 448 effmass_free = one; ecutsm = zero; ecut = huge(one) 449 call mkkin(ecut,ecutsm,effmass_free,cryst%gmet,kg_k,modkplusg,kpoint,npw_k,0,0) 450 modkplusg(:) = SQRT(half/pi**2*modkplusg(:)) 451 modkplusg(:) = MAX(modkplusg(:),tol10) 452 453 ! Calculate matrix elements. 454 vkb=zero; vkbd=zero 455 456 do itypat=1,psps%ntypat 457 iln0 = 0 458 nlmn = count(psps%indlmn(3,:,itypat) > 0) 459 do ilmn=1,nlmn 460 il = 1 + psps%indlmn(1,ilmn,itypat) 461 in = psps%indlmn(3,ilmn,itypat) 462 iln = psps%indlmn(5,ilmn,itypat) 463 !write(*,*)ilmn, iln, il, in 464 if (iln <= iln0) cycle 465 iln0 = iln 466 !if (vkbsign(iln,itypat) == zero) cycle 467 if (ABS(psps%ekb(iln,itypat)) > 1.0d-10) then 468 ABI_CHECK(iln == ilmn, "iln != ilmn") 469 !ABI_CHECK(il == iln, "il != iln") 470 if (il==1) then 471 vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat) 472 vkbd(1:npw_k,iln,itypat) = ffnl(:,2,iln,itypat)*modkplusg(:)/two_pi 473 else if (il==2) then 474 vkb(1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)*modkplusg(:) 475 do ig=1,npw_k 476 vkbd(ig,iln,itypat) = ((ffnl(ig,2,iln,itypat)*modkplusg(ig)*modkplusg(ig))+& 477 ffnl(ig,1,iln,itypat) )/two_pi 478 end do 479 else if (il==3) then 480 vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)*modkplusg(:)**2 481 vkbd(1:npw_k,iln,itypat) = (ffnl(:,2,iln,itypat)*modkplusg(:)**3+& 482 2*ffnl(:,1,iln,itypat)*modkplusg(:) )/two_pi 483 else if (il==4) then 484 vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)*modkplusg(:)**3 485 vkbd(1:npw_k,iln,itypat) = (ffnl(:,2,iln,itypat)*modkplusg(:)**4+& 486 3*ffnl(:,1,iln,itypat)*modkplusg(:)**2 )/two_pi 487 end if 488 vkb (:,iln,itypat) = SQRT(4*pi/cryst%ucvol*(2*il-1)*ABS(psps%ekb(iln,itypat)))*vkb (:,iln,itypat) 489 vkbd(:,iln,itypat) = SQRT(4*pi/cryst%ucvol*(2*il-1)*ABS(psps%ekb(iln,itypat)))*vkbd(:,iln,itypat) 490 end if 491 end do 492 end do 493 494 ABI_FREE(ffnl) 495 ABI_FREE(modkplusg) 496 497 DBG_EXIT("COLL") 498 499 end subroutine calc_vkb
m_vkbr/ccgradvnl_ylm [ Functions ]
[ Top ] [ m_vkbr ] [ Functions ]
NAME
ccgradvnl_ylm
FUNCTION
Compute Vnl(K) and grad_K Vnl(K) three reciprocal lattice units components using spherical harmonics instead of Legendre polynomials Needed for chi0(q=0)
INPUTS
cryst<crystal_t>=Unit cell and symmetries psps<pseudopotential_type>Structure gathering info on the pseudopotentials. npw=number of planewaves for wavefunctions gvec(3,npw)=integer coordinates of each plane wave in reciprocal space kpoint(3)=K-point in reduced coordinates. vkbsign(lnmax,ntypat)=sign of each KB dyadic product vkb(npw,lnmax,ntypat)=KB projector function vkbd(npw,lnmax,ntypat)=derivative of the KB projector function in reciprocal space
OUTPUT
fnl(npw,mpsang*2,natom), fnld(3,npw,mpsang*2,natom)
NOTES
Subroutine taken from the EXC code All the calculations are done in double precision, but the output arrays fnl and fnld are in single precision, should use double precision after modification of the other subroutines
SOURCE
634 subroutine ccgradvnl_ylm(cryst,psps,npw,gvec,kpoint,vkbsign,vkb,vkbd,fnl,fnld) 635 636 !Arguments ------------------------------------ 637 !scalars 638 integer,intent(in) :: npw 639 type(crystal_t),intent(in) :: cryst 640 type(pseudopotential_type),intent(in) :: psps 641 !arrays 642 integer,intent(in) :: gvec(3,npw) 643 real(dp),intent(in) :: kpoint(3) 644 real(dp),intent(in) :: vkb(npw,psps%lnmax,cryst%ntypat) 645 real(dp),intent(in) :: vkbd(npw,psps%lnmax,cryst%ntypat) 646 real(dp),intent(in) :: vkbsign(psps%lnmax,cryst%ntypat) 647 complex(gwpc),intent(out) :: fnl(npw,psps%mpsang**2,psps%mproj,cryst%natom) 648 complex(gwpc),intent(out) :: fnld(3,npw,psps%mpsang**2,psps%mproj,cryst%natom) 649 650 !Local variables------------------------------- 651 !scalars 652 integer :: ii,iat,ig,il,im,ilm,itypat,nlmn,iln0,iln,ilmn,in 653 real(dp),parameter :: ppad=tol6 654 real(dp) :: cosphi,costh,factor,mkg,mkg2,sinphi,sinth,sq,xdotg 655 complex(dpc) :: dphi,dth,sfac 656 character(len=500) :: msg 657 !arrays 658 real(dp) :: gcart(3),kcart(3),kg(3) 659 real(dp) :: b1(3),b2(3),b3(3),a1(3),a2(3),a3(3) 660 complex(dpc) :: dylmcart(3),dylmcrys(3),gradphi(3),gradth(3) 661 !************************************************************************ 662 663 DBG_ENTER("COLL") 664 665 if (psps%mpsang > 4) then 666 write(msg,'(3a)')& 667 'Number of angular momentum components bigger than programmed.',ch10,& 668 'Taking into account only s p d f ' 669 ABI_ERROR(msg) 670 end if 671 672 a1=cryst%rprimd(:,1); b1=two_pi*Cryst%gprimd(:,1) 673 a2=cryst%rprimd(:,2); b2=two_pi*Cryst%gprimd(:,2) 674 a3=cryst%rprimd(:,3); b3=two_pi*Cryst%gprimd(:,3) 675 676 ! Calculate Kleiman-Bylander factor and first derivative. 677 fnl=czero_gw; fnld=czero_gw 678 679 do ig=1,npw 680 ! Get kcart = k+G in Cartesian coordinates. 681 kg(:)= kpoint(:) + REAL(gvec(:,ig)) 682 kcart(:) = kg(1)*b1(:) + kg(2)*b2(:) + kg(3)*b3(:) 683 ! Solve the problem with sinth=0. or sinphi=0 684 if (ABS(kcart(2))<ppad) kcart(2) = kcart(2) + ppad 685 686 mkg2 = kcart(1)**2+kcart(2)**2+kcart(3)**2 687 mkg = SQRT(mkg2) 688 ! The next to solve the problem with k=Gamma. 689 !if (mkg < 0.0001) cycle 690 691 sq=SQRT(kcart(1)**2+kcart(2)**2) 692 693 gcart(:)= REAL(gvec(1,ig))*b1(:)& 694 & +REAL(gvec(2,ig))*b2(:)& 695 & +REAL(gvec(3,ig))*b3(:) 696 697 ! Calculate spherical coordinates (th, phi). 698 costh = kcart(3)/mkg 699 sinth = sq/mkg 700 cosphi= kcart(1)/sq 701 sinphi= kcart(2)/sq 702 703 gradth(1) = kcart(1)*kcart(3)/mkg**3/sinth 704 gradth(2) = kcart(2)*kcart(3)/mkg**3/sinth 705 gradth(3) = -(one/mkg-kcart(3)**2/mkg**3)/sinth 706 gradphi(1) = -(one/sq - kcart(1)**2/sq**3)/sinphi 707 gradphi(2) = kcart(2)*kcart(1)/sq**3/sinphi 708 gradphi(3) = czero 709 710 do iat=1,cryst%natom 711 itypat = cryst%typat(iat) 712 xdotg = gcart(1)*cryst%xcart(1,iat)+gcart(2)*Cryst%xcart(2,iat)+gcart(3)*Cryst%xcart(3,iat) 713 ! Remember that in the GW code the reciprocal vectors 714 ! are defined such as a_i*b_j = 2pi delta_ij, no need to introduce 2pi 715 sfac=CMPLX(COS(xdotg), SIN(xdotg), kind=dpc) 716 717 iln0 = 0 718 nlmn = count(psps%indlmn(3,:,itypat) > 0) 719 do ilmn=1,nlmn 720 il = 1 + psps%indlmn(1,ilmn,itypat) 721 in = psps%indlmn(3,ilmn,itypat) 722 iln = psps%indlmn(5,ilmn,itypat) 723 ! spin = 1 if scalar term (spin diagonal), 2 if SOC term. 724 !spin = psps%indlmn(6, ilmn, itypat) 725 if (iln <= iln0) cycle 726 iln0 = iln 727 if (vkbsign(iln,itypat) == zero) cycle 728 !if (spin /= 1 .or. vkbsign(iln,itypat) == zero) cycle 729 factor = SQRT(four_pi/REAL(2*(il-1)+1)) 730 do im=1,2*(il-1)+1 731 ! Index of im and il 732 ilm = im + (il-1)*(il-1) 733 734 ! Calculate the first KB factor, note that fnl is simple precision complex 735 fnl(ig,ilm,in,iat) = factor*sfac*ylmc(il-1,im-il,kcart) * vkb(ig,iln,itypat) * vkbsign(iln,itypat) 736 737 ! Calculate the second KB factor (involving first derivatives) 738 ! dYlm/dK = dYlm/dth * grad_K th + dYlm/dphi + grad_K phi 739 call ylmcd(il-1,im-il,kcart,dth,dphi) 740 dylmcart(:) = dth*gradth(:) + dphi*gradphi(:) 741 742 ! Cartesian to crystallographic axis 743 ! Notice: a bug was discovered by Marco Cazzaniga, december 2009 744 ! the transformation matrix A=(a1,a2,a3) must act *on its left* on the 745 ! covariant vector dylmcart (a *row* vector). The previous implementation assumed A 746 ! acting on its right on a column vector, yielding wrong results for the (small) 747 ! non local contributions to the spectra, such as a spurious anisotropy in isotropic systems. 748 ! This is the correct version: 749 dylmcrys(1) = (a1(1)*dylmcart(1)+a1(2)*dylmcart(2)+a1(3)*dylmcart(3))/(two_pi) 750 dylmcrys(2) = (a2(1)*dylmcart(1)+a2(2)*dylmcart(2)+a2(3)*dylmcart(3))/(two_pi) 751 dylmcrys(3) = (a3(1)*dylmcart(1)+a3(2)*dylmcart(2)+a3(3)*dylmcart(3))/(two_pi) 752 753 ! Note that fnld is simple precision complex, it could be possible to use double precision 754 do ii=1,3 755 fnld(ii,ig,ilm,in,iat) = factor*sfac* & 756 ( kg(ii)/mkg*ylmc(il-1,im-il,kcart)*vkbd(ig,iln,itypat) + dylmcrys(ii)*vkb(ig,iln,itypat) ) 757 end do 758 759 end do !im 760 end do !il 761 end do !iat 762 end do !ig 763 764 DBG_EXIT("COLL") 765 766 end subroutine ccgradvnl_ylm
m_vkbr/nc_ihr_comm [ Functions ]
[ Top ] [ m_vkbr ] [ Functions ]
NAME
nc_ihr_comm
FUNCTION
Calculate the matrix elements of the commutator i[H,r] For NC pseudppotentials, the commutator i[Vnl,r] is included depending on inclvkb.
INPUTS
vkbr<vkbr_t> cryst<crystal_t>=Unit cell and symmetries psps<pseudopotential_type>Structure gathering info on the pseudopotentials. nspinor=Number of spinorial components. npw=Number of G for wavefunctions. istwfk=Storage mode for wavefunctions. inclvkb=Option defining whether [Vnl,r] is added or not. kpoint(3)=k-point in reduced coordinates. ug1(npw*nspinor)=Left wavefunction. ug2(npw*nspinor)=Right wavefunction gvec(3,npw)=Planes waves for wavefunctions.
OUTPUT
ihr_comm(3,nspinor**2)= Matrix elements of the commutator i[H,r] between the input states. Result is in reduced coordinates. ug1 and ug2 are supposed to be orthogonal.
NOTES
<k b1|e^{-iq.r}|k b2> = \delta_{b1 b2} -iq <k b1|r|k b2> = \delta_{b1 b2} -iq ( <k b1| [H,r] |k b2> / (e1-e2) ).
Remember that [H,r] = -\nabla + [V_nl,r]
TODO
*) Spinorial case is not implemented.
SOURCE
538 function nc_ihr_comm(vkbr,cryst,psps,npw,nspinor,istwfk,inclvkb,kpoint,ug1,ug2,gvec) result(ihr_comm) 539 540 !Arguments ------------------------------------ 541 !scalars 542 integer,intent(in) :: npw,nspinor,inclvkb,istwfk 543 type(vkbr_t),intent(in) :: vkbr 544 type(crystal_t),intent(in) :: cryst 545 type(Pseudopotential_type),intent(in) :: psps 546 !arrays 547 integer,intent(in) :: gvec(3,npw) 548 real(dp),intent(in) :: kpoint(3) 549 complex(gwpc),intent(in) :: ug1(npw*nspinor),ug2(npw*nspinor) 550 complex(gwpc) :: ihr_comm(3,nspinor**2) 551 552 !Local variables ------------------------------ 553 !scalars 554 integer :: ig,iab,spad1,spad2 555 complex(dpc) :: c_tmp 556 !arrays 557 integer :: spinorwf_pad(2,4) 558 559 !************************************************************************ 560 561 ! [H, r] = -\nabla + [V_{nl}, r] 562 ! V_nl is present only in the case of NC pseudos but 563 ! not in PAW unless even the AE Hamiltonian in non-local e.g. DFT+U or LEXX. 564 565 ! -i <c,k|\nabla_r|v,k> in reduced coordinates is always included. 566 ! -i <c,k|\nabla_r|v,k> = \sum_G u_{ck}^*(G) [k+G] u_{vk}(G) 567 ! Note that here we assume c/=v, moreover the ug are supposed to be orthonormal and 568 ! hence k+G can be replaced by G. 569 ! HM 03/08/2018: we need band velocities so we don't assume c/=v anymore and we use k+G. 570 571 spinorwf_pad = RESHAPE([0, 0, npw, npw, 0, npw, npw, 0], [2, 4]) 572 ihr_comm = czero 573 574 ! -i <c,k|\nabla_r|v,k> in reduced coordinates. 575 ! This term is spin diagonal if nspinor == 2 576 if (istwfk == 1) then 577 do iab=1,nspinor 578 spad1 = spinorwf_pad(1,iab); spad2 = spinorwf_pad(2,iab) 579 do ig=1,npw 580 c_tmp = GWPC_CONJG(ug1(ig+spad1)) * ug2(ig+spad2) 581 ihr_comm(:,iab) = ihr_comm(:,iab) + c_tmp * (kpoint + gvec(:,ig)) 582 end do 583 end do 584 else 585 ! Symmetrized expression: \sum_G (k+G) 2i Im [ u_a^*(G) u_b(G) ]. (k0, G0) term is null. 586 ABI_CHECK(nspinor == 1, "nspinor != 1") 587 do ig=1,npw 588 c_tmp = GWPC_CONJG(ug1(ig)) * ug2(ig) 589 ihr_comm(:,1) = ihr_comm(:,1) + two*j_dpc * AIMAG(c_tmp) * (kpoint + gvec(:,ig)) 590 end do 591 end if 592 593 ! Add second term $i <c,k|[Vnl,r]|v,k> $ in reduced cordinates. 594 if (inclvkb /= 0) then 595 ABI_CHECK(istwfk == vkbr%istwfk, "input istwfk /= vkbr%istwfk") 596 call add_vnlr_commutator(vkbr,cryst,psps,npw,nspinor,ug1,ug2,ihr_comm) 597 end if 598 599 end function nc_ihr_comm
m_vkbr/vkbr_free_0D [ Functions ]
[ Top ] [ m_vkbr ] [ Functions ]
NAME
vkbr_free_0D
FUNCTION
Free all memory allocated in a structure of type vkbr_t
SOURCE
207 subroutine vkbr_free_0D(vkbr) 208 209 !Arguments ------------------------------------ 210 !scalars 211 type(vkbr_t),intent(inout) :: vkbr 212 213 !************************************************************************ 214 215 !complex 216 ABI_SFREE(vkbr%fnl) 217 ABI_SFREE(vkbr%fnld) 218 219 end subroutine vkbr_free_0D
m_vkbr/vkbr_free_1D [ Functions ]
[ Top ] [ m_vkbr ] [ Functions ]
NAME
vkbr_free_1D
FUNCTION
Free all memory allocated in a structure of type vkbr_t
SOURCE
233 subroutine vkbr_free_1D(vkbr) 234 235 !Arguments ------------------------------------ 236 !arrays 237 type(vkbr_t),intent(inout) :: vkbr(:) 238 239 !Local variables ------------------------------ 240 !scalars 241 integer :: ii 242 243 !************************************************************************ 244 245 do ii=1,SIZE(vkbr) 246 call vkbr_free_0D(vkbr(ii)) 247 end do 248 249 end subroutine vkbr_free_1D
m_vkbr/vkbr_init [ Functions ]
[ Top ] [ m_vkbr ] [ Functions ]
NAME
vkbr_init
FUNCTION
Creation method the the vkbr_t structures datatype.
INPUTS
cryst<crystal_t>=Datatype gathering info on the crystal structure. psps<pseudopotential_type>Structure gathering info on the pseudopotentials. inclvkb=Option defining the algorithm used for the application of [Vnl,r]. 2 for Spherical harmonics istwfk=Storage mode for the wavefunctions at this k-point. npw=Number of planewaves in <k+G1|[Vnl,r]|k+G2> kpoint(3)=K-point of interest in reduced coordinates. gvec(3,npw)=Reduced coordinates of the G-vectors.
OUTPUT
vkbr<vkbr_t>=Structure containing arrays needed for calculating <\psi_1|[Vnl,r]\psi_2>. Completely initialized in output.
SOURCE
130 subroutine vkbr_init(vkbr,cryst,psps,inclvkb,istwfk,npw,kpoint,gvec) 131 132 !Arguments ------------------------------------ 133 !scalars 134 integer,intent(in) :: npw,inclvkb,istwfk 135 type(crystal_t),intent(in) :: cryst 136 type(vkbr_t),intent(inout) :: vkbr 137 type(pseudopotential_type),intent(in) :: psps 138 !arrays 139 integer,intent(in) :: gvec(3,npw) 140 real(dp),intent(in) :: kpoint(3) 141 142 !Local variables------------------------------- 143 !scalars 144 integer :: ierr 145 character(len=500) :: msg 146 !arrays 147 real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:) 148 149 !************************************************************************ 150 151 !@vkbr_t 152 vkbr%istwfk = istwfk 153 vkbr%ntypat = cryst%ntypat 154 vkbr%natom = cryst%natom 155 vkbr%mpsang = psps%mpsang 156 vkbr%npw = npw 157 vkbr%inclvkb = inclvkb 158 vkbr%kpoint = kpoint 159 160 ! Calculate KB form factors and derivatives. 161 ! The arrays are allocated with lnmax to support pseudos with more than projector. 162 ! Note that lnmax takes into account lloc hence arrays are in packed form and one should be 163 ! accessed with the indices provided by indlmn. 164 ! TODO: they should be calculated on-the-fly using calc_vkb 165 ! For the moment, we opt for a quick an dirty implementation. 166 167 ABI_MALLOC(vkbsign, (psps%lnmax, cryst%ntypat)) 168 ABI_MALLOC(vkb, (npw, psps%lnmax, cryst%ntypat)) 169 ABI_MALLOC(vkbd, (npw, psps%lnmax, cryst%ntypat)) 170 call calc_vkb(cryst,psps,kpoint,npw,npw,gvec,vkbsign,vkb,vkbd) 171 172 select case (inclvkb) 173 case (2) 174 ! Complex spherical harmonics (CPU and mem \propto npw). 175 write(msg,'(a,f12.1)')'out-of-memory in fnl; Mb= ',one*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwpc*b2Mb 176 ABI_STAT_MALLOC(vkbr%fnl,(npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr) 177 ABI_CHECK(ierr==0, msg) 178 179 write(msg,'(a,f12.1)')'out-of-memory in fnld; Mb= ',three*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwpc*b2Mb 180 ABI_STAT_MALLOC(vkbr%fnld,(3,npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr) 181 ABI_CHECK(ierr==0, msg) 182 183 call ccgradvnl_ylm(cryst,psps,npw,gvec,kpoint,vkbsign,vkb,vkbd,vkbr%fnl,vkbr%fnld) 184 185 case default 186 ABI_ERROR(sjoin("Wrong inclvkb= ",itoa(inclvkb))) 187 end select 188 189 ABI_FREE(vkbsign) 190 ABI_FREE(vkb) 191 ABI_FREE(vkbd) 192 193 end subroutine vkbr_init
m_vkbr/vkbr_t [ Types ]
NAME
FUNCTION
Matrix elements in |k+G> space needed for the evaluation of the matrix elements of the commutator [Vnl,r] for the optical limit in <kb1|e^{-iqr}|kb2>.
SOURCE
60 type,public :: vkbr_t 61 62 integer :: istwfk 63 ! Storage mode of the G vectors for this k-point. 64 65 integer :: ntypat 66 ! Number of type of atoms 67 68 integer :: natom 69 ! Number of atoms 70 71 integer :: mpsang 72 ! Max l+1 over atoms 73 74 integer :: npw 75 ! Number of G-vectors. 76 77 integer :: inclvkb 78 ! Option for calculating the matrix elements of [Vnl,r]. 79 ! 0 to exclude commutator, 2 to include it 80 81 real(dp) :: kpoint(3) 82 ! The k-point in reduced coordinates. 83 84 complex(gwpc),allocatable :: fnl(:,:,:,:) 85 ! fnl(npw,mpsang**2,mproj,natom) 86 87 complex(gwpc),allocatable :: fnld(:,:,:,:,:) 88 ! fnld(3,npw,mpsang**2,mproj,natom) 89 90 end type vkbr_t 91 92 public :: vkbr_init ! vkbr_t Constructor 93 public :: vkbr_free ! Free memory 94 public :: nc_ihr_comm ! Compute matrix elements of the commutator i[H,r] for NC pseudos 95 public :: calc_vkb ! Kleynman-Bylander form factors and derivatives.