TABLE OF CONTENTS


ABINIT/m_vkbr [ Modules ]

[ Top ] [ 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.

COPYRIGHT

 Copyright (C) 2008-2018 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

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 MODULE m_vkbr
26 
27  use defs_basis
28  use defs_datatypes
29  use m_hide_blas
30  use m_errors
31  use m_abicore
32 
33  use m_gwdefs,        only : czero_gw
34  use m_fstrings,      only : sjoin, itoa
35  use m_paw_sphharm,   only : ylmc, ylmcd
36  use m_geometry,      only : normv
37  use m_crystal,       only : crystal_t
38  use m_kg,            only : mkkin
39  use m_mkffnl,        only : mkffnl
40 
41  implicit none
42 
43  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.

PARENTS

      m_vkbr

CHILDREN

      ylmcd

SOURCE

340 subroutine add_vnlr_commutator(vkbr,cryst,psps,npw,nspinor,ug1,ug2,rhotwx)
341 
342 
343 !This section has been created automatically by the script Abilint (TD).
344 !Do not modify the following lines by hand.
345 #undef ABI_FUNC
346 #define ABI_FUNC 'add_vnlr_commutator'
347 !End of the abilint section
348 
349  implicit none
350 
351 !Arguments ------------------------------------
352 !scalars
353  integer,intent(in) :: npw,nspinor
354  type(vkbr_t),intent(in) :: vkbr
355  type(crystal_t),intent(in) :: cryst
356  type(pseudopotential_type),intent(in) :: psps
357 !arrays
358  complex(gwpc),target,intent(in) :: ug1(npw*nspinor),ug2(npw*nspinor)
359  complex(gwpc),intent(inout) :: rhotwx(3,nspinor**2)
360 
361 !Local variables ------------------------------
362 !scalars
363  integer :: iat,ig,ilm,itypat,nlmn,ilmn,iln0,iln,il,in,im
364  complex(gwpc) :: cta1,cta4
365 !arrays
366  complex(gwpc) :: dum(3),cta2(3),cta3(3),gamma_term(3)
367 
368 !************************************************************************
369 
370  ABI_CHECK(nspinor == 1, "nspinor/=1 not coded")
371 
372  ! Adding term i <c,k|[Vnl,r]|v,k> ===
373  select case (vkbr%inclvkb)
374  case (2)
375   ! Complex spherical harmonics (much faster!).
376   dum=czero_gw; gamma_term=czero
377 
378   do iat=1,vkbr%natom
379     itypat = cryst%typat(iat)
380     nlmn = count(psps%indlmn(3,:,itypat) > 0)
381     iln0 = 0
382     do ilmn=1,nlmn
383       il = 1 + psps%indlmn(1,ilmn,itypat)
384       in = psps%indlmn(3,ilmn,itypat)
385       iln = psps%indlmn(5,ilmn,itypat)
386       if (iln <= iln0) cycle
387       iln0 = iln
388       !if (indlmn(6,ilmn,itypat) /= 1 .or. vkbsign(iln,itypat) == zero) cycle
389       !in = 1
390       do im=1,2*(il-1)+1
391         ! Index of im and il
392         ilm = im + (il-1)*(il-1)
393         cta1 = czero_gw; cta2(:) = czero_gw
394         cta4 = czero_gw; cta3(:) = czero_gw
395         do ig=1,npw
396           ! Here we take advantage of the property Y_(l-m)= (-i)^m Y_lm^*.
397           cta1   = cta1    + ug1(ig) * vkbr%fnl (ig,ilm,in,iat)
398           cta2(:)= cta2(:) + ug2(ig) * vkbr%fnld(:,ig,ilm,in,iat)
399           cta3(:)= cta3(:) + ug1(ig) * vkbr%fnld(:,ig,ilm,in,iat)
400           cta4   = cta4    + ug2(ig) * vkbr%fnl (ig,ilm,in,iat)
401           if (ig==1) gamma_term = gamma_term + CONJG(cta1)*cta2(:) +CONJG(cta3(:))*cta4
402         end do
403         dum(:)= dum(:) + CONJG(cta1)*cta2(:) + CONJG(cta3(:))*cta4
404       end do
405 
406     end do
407   end do
408 
409   if (vkbr%istwfk>1) then
410     dum = two * j_dpc * AIMAG(dum); if (vkbr%istwfk==2) dum = dum - j_dpc * AIMAG(gamma_term)
411   end if
412   rhotwx(:,1) = rhotwx(:,1) + dum(:)
413 
414  case default
415    MSG_ERROR(sjoin("Wrong inclvkb:", itoa(vkbr%inclvkb)))
416  end select
417 
418 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.

PARENTS

      m_vkbr

CHILDREN

      ylmcd

SOURCE

455 subroutine calc_vkb(cryst,psps,kpoint,npw_k,mpw,kg_k,vkbsign,vkb,vkbd)
456 
457 
458 !This section has been created automatically by the script Abilint (TD).
459 !Do not modify the following lines by hand.
460 #undef ABI_FUNC
461 #define ABI_FUNC 'calc_vkb'
462 !End of the abilint section
463 
464  implicit none
465 
466 !Arguments ------------------------------------
467 !scalars
468  integer,intent(in) :: npw_k, mpw
469  type(crystal_t),intent(in) :: cryst
470  type(pseudopotential_type),intent(in) :: psps
471 !arrays
472  integer,intent(in) :: kg_k(3,npw_k)
473  real(dp),intent(in) :: kpoint(3)
474  real(dp),intent(out) :: vkb (mpw,psps%lnmax,psps%ntypat)
475  real(dp),intent(out) :: vkbd(mpw,psps%lnmax,psps%ntypat)
476  real(dp),intent(out) :: vkbsign(psps%lnmax,psps%ntypat)
477 
478 !Local variables ------------------------------
479 !scalars
480  integer :: dimffnl,ider,idir,itypat,nkpg,in,il,ilmn,ig,iln,iln0,nlmn
481  real(dp) :: effmass_free,ecutsm,ecut
482 !arrays
483  real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:),ylm_gr(:,:,:),ylm_k(:,:)
484 
485 ! *************************************************************************
486 
487  DBG_ENTER("COLL")
488  ABI_CHECK(psps%usepaw==0, "You should not be here!")
489  ABI_CHECK(psps%useylm==0, "useylm/=0 not considered!")
490 
491  ! Compute KB dyadic sign.
492  vkbsign=zero
493  do itypat=1,psps%ntypat
494    iln0 = 0
495    nlmn = count(psps%indlmn(3,:,itypat) > 0)
496    do ilmn=1,nlmn
497      iln = psps%indlmn(5,ilmn,itypat)
498      if (iln <= iln0) cycle
499      iln0 = iln
500      if (abs(psps%ekb(iln,itypat)) > 1.0d-10) vkbsign(iln,itypat) = dsign(one, psps%ekb(iln,itypat))
501    end do
502  end do
503 
504  ! Allocate KB form factor and derivative wrt k+G
505  ! Here we do not use correct ordering for dimensions
506  idir=0; nkpg=0; ider=1; dimffnl=2 ! To retrieve the first derivative.
507 
508  ! Quantities used only if useylm==1
509  ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2*Psps%useylm))
510  ABI_MALLOC(ylm_gr, (npw_k, 3+6*(ider/2),psps%mpsang**2*Psps%useylm))
511  ABI_MALLOC(kpg_dum, (npw_k, nkpg))
512  ABI_MALLOC(ffnl, (npw_k,dimffnl, psps%lmnmax, psps%ntypat))
513 
514  call mkffnl(psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,cryst%gmet,cryst%gprimd,ider,idir,Psps%indlmn,&
515    kg_k,kpg_dum,kpoint,psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,&
516    psps%ntypat,Psps%pspso,Psps%qgrid_ff,cryst%rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr)
517 
518  ABI_FREE(ylm_k)
519  ABI_FREE(ylm_gr)
520  ABI_FREE(kpg_dum)
521 
522  ABI_MALLOC(modkplusg, (npw_k))
523  effmass_free = one; ecutsm = zero; ecut = huge(one)
524  call mkkin(ecut,ecutsm,effmass_free,cryst%gmet,kg_k,modkplusg,kpoint,npw_k,0,0)
525  modkplusg(:) = SQRT(half/pi**2*modkplusg(:))
526  modkplusg(:) = MAX(modkplusg(:),tol10)
527 
528  ! Calculate matrix elements.
529  vkb=zero; vkbd=zero
530 
531  do itypat=1,psps%ntypat
532    iln0 = 0
533    nlmn = count(psps%indlmn(3,:,itypat) > 0)
534    do ilmn=1,nlmn
535      il = 1 + psps%indlmn(1,ilmn,itypat)
536      in = psps%indlmn(3,ilmn,itypat)
537      iln = psps%indlmn(5,ilmn,itypat)
538      !write(*,*)ilmn, iln, il, in
539      if (iln <= iln0) cycle
540      iln0 = iln
541      !if (vkbsign(iln,itypat) == zero) cycle
542      if (ABS(psps%ekb(iln,itypat)) > 1.0d-10) then
543        ABI_CHECK(iln == ilmn, "iln != ilmn")
544        !ABI_CHECK(il == iln, "il != iln")
545        if (il==1) then
546          vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)
547          vkbd(1:npw_k,iln,itypat) = ffnl(:,2,iln,itypat)*modkplusg(:)/two_pi
548        else if (il==2) then
549          vkb(1:npw_k,iln,itypat)  = ffnl(:,1,iln,itypat)*modkplusg(:)
550          do ig=1,npw_k
551            vkbd(ig,iln,itypat) = ((ffnl(ig,2,iln,itypat)*modkplusg(ig)*modkplusg(ig))+&
552             ffnl(ig,1,iln,itypat) )/two_pi
553          end do
554        else if (il==3) then
555          vkb (1:npw_k,iln,itypat) =  ffnl(:,1,iln,itypat)*modkplusg(:)**2
556          vkbd(1:npw_k,iln,itypat) = (ffnl(:,2,iln,itypat)*modkplusg(:)**3+&
557           2*ffnl(:,1,iln,itypat)*modkplusg(:) )/two_pi
558        else if (il==4) then
559          vkb (1:npw_k,iln,itypat) =  ffnl(:,1,iln,itypat)*modkplusg(:)**3
560          vkbd(1:npw_k,iln,itypat) = (ffnl(:,2,iln,itypat)*modkplusg(:)**4+&
561           3*ffnl(:,1,iln,itypat)*modkplusg(:)**2 )/two_pi
562        end if
563        vkb (:,iln,itypat) = SQRT(4*pi/cryst%ucvol*(2*il-1)*ABS(psps%ekb(iln,itypat)))*vkb (:,iln,itypat)
564        vkbd(:,iln,itypat) = SQRT(4*pi/cryst%ucvol*(2*il-1)*ABS(psps%ekb(iln,itypat)))*vkbd(:,iln,itypat)
565      end if
566    end do
567  end do
568 
569  ABI_FREE(ffnl)
570  ABI_FREE(modkplusg)
571 
572  DBG_EXIT("COLL")
573 
574 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

PARENTS

      m_vkbr

CHILDREN

      ylmcd

SOURCE

729 subroutine ccgradvnl_ylm(cryst,psps,npw,gvec,kpoint,vkbsign,vkb,vkbd,fnl,fnld)
730 
731 
732 !This section has been created automatically by the script Abilint (TD).
733 !Do not modify the following lines by hand.
734 #undef ABI_FUNC
735 #define ABI_FUNC 'ccgradvnl_ylm'
736 !End of the abilint section
737 
738  implicit none
739 
740 !Arguments ------------------------------------
741 !scalars
742  integer,intent(in) :: npw
743  type(crystal_t),intent(in) :: cryst
744  type(pseudopotential_type),intent(in) :: psps
745 !arrays
746  integer,intent(in) :: gvec(3,npw)
747  real(dp),intent(in) :: kpoint(3)
748  real(dp),intent(in) :: vkb(npw,psps%lnmax,cryst%ntypat)
749  real(dp),intent(in) :: vkbd(npw,psps%lnmax,cryst%ntypat)
750  real(dp),intent(in) :: vkbsign(psps%lnmax,cryst%ntypat)
751  complex(gwpc),intent(out) :: fnl(npw,psps%mpsang**2,psps%mproj,cryst%natom)
752  complex(gwpc),intent(out) :: fnld(3,npw,psps%mpsang**2,psps%mproj,cryst%natom)
753 
754 !Local variables-------------------------------
755 !scalars
756  integer :: ii,iat,ig,il,im,ilm,itypat,nlmn,iln0,iln,ilmn,in
757  real(dp),parameter :: ppad=tol8
758  real(dp) :: cosphi,costh,factor,mkg,mkg2,sinphi,sinth,sq,xdotg
759  complex(dpc) :: dphi,dth,sfac
760  character(len=500) :: msg
761 !arrays
762  real(dp) :: gcart(3),kcart(3),kg(3)
763  real(dp) :: b1(3),b2(3),b3(3),a1(3),a2(3),a3(3)
764  complex(dpc) :: dylmcart(3),dylmcrys(3),gradphi(3),gradth(3)
765 !************************************************************************
766 
767  DBG_ENTER("COLL")
768 
769  if (psps%mpsang > 4) then
770    write(msg,'(3a)')&
771     'Number of angular momentum components bigger than programmed.',ch10,&
772     'Taking into account only s p d f '
773    MSG_ERROR(msg)
774  end if
775 
776  a1=cryst%rprimd(:,1); b1=two_pi*Cryst%gprimd(:,1)
777  a2=cryst%rprimd(:,2); b2=two_pi*Cryst%gprimd(:,2)
778  a3=cryst%rprimd(:,3); b3=two_pi*Cryst%gprimd(:,3)
779 
780  ! Calculate Kleiman-Bylander factor and first derivative.
781  fnl=czero_gw; fnld=czero_gw
782 
783  do ig=1,npw
784    ! Get kcart = k+G in Cartesian coordinates.
785    kg(:)= kpoint(:) + REAL(gvec(:,ig))
786    kcart(:) = kg(1)*b1(:) + kg(2)*b2(:) + kg(3)*b3(:)
787    ! Solve the problem with sinth=0. or sinphi=0
788    if (ABS(kcart(2))<ppad) kcart(2) = kcart(2) + ppad
789 
790    mkg2 = kcart(1)**2+kcart(2)**2+kcart(3)**2
791    mkg = SQRT(mkg2)
792    ! The next to solve the problem with k=Gamma.
793    !if (mkg < 0.0001) cycle
794 
795    sq=SQRT(kcart(1)**2+kcart(2)**2)
796 
797    gcart(:)=  REAL(gvec(1,ig))*b1(:)&
798 &            +REAL(gvec(2,ig))*b2(:)&
799 &            +REAL(gvec(3,ig))*b3(:)
800 
801    ! Calculate spherical coordinates (th, phi).
802    costh = kcart(3)/mkg
803    sinth = sq/mkg
804    cosphi= kcart(1)/sq
805    sinphi= kcart(2)/sq
806 
807    gradth(1)  = kcart(1)*kcart(3)/mkg**3/sinth
808    gradth(2)  = kcart(2)*kcart(3)/mkg**3/sinth
809    gradth(3)  = -(one/mkg-kcart(3)**2/mkg**3)/sinth
810    gradphi(1) = -(one/sq - kcart(1)**2/sq**3)/sinphi
811    gradphi(2) = kcart(2)*kcart(1)/sq**3/sinphi
812    gradphi(3) = czero
813 
814    do iat=1,cryst%natom
815      itypat = cryst%typat(iat)
816      xdotg = gcart(1)*cryst%xcart(1,iat)+gcart(2)*Cryst%xcart(2,iat)+gcart(3)*Cryst%xcart(3,iat)
817      ! Remember that in the GW code the reciprocal vectors
818      ! are defined such as a_i*b_j = 2pi delta_ij, no need to introduce 2pi
819      sfac=CMPLX(COS(xdotg), SIN(xdotg))
820 
821      iln0 = 0
822      nlmn = count(psps%indlmn(3,:,itypat) > 0)
823      do ilmn=1,nlmn
824        il = 1 + psps%indlmn(1,ilmn,itypat)
825        in = psps%indlmn(3,ilmn,itypat)
826        iln = psps%indlmn(5,ilmn,itypat)
827        ! spin = 1 if scalar term (spin diagonal), 2 if SOC term.
828        !spin = psps%indlmn(6, ilmn, itypat)
829        if (iln <= iln0) cycle
830        iln0 = iln
831        if (vkbsign(iln,itypat) == zero) cycle
832        !if (spin /= 1 .or. vkbsign(iln,itypat) == zero) cycle
833        factor = SQRT(four_pi/REAL(2*(il-1)+1))
834        do im=1,2*(il-1)+1
835          ! Index of im and il
836          ilm = im + (il-1)*(il-1)
837 
838          ! Calculate the first KB factor, note that fnl is simple precision complex
839          fnl(ig,ilm,in,iat) = factor*sfac*ylmc(il-1,im-il,kcart) * vkb(ig,iln,itypat) * vkbsign(iln,itypat)
840 
841          ! Calculate the second KB factor (involving first derivatives)
842          ! dYlm/dK = dYlm/dth * grad_K th + dYlm/dphi + grad_K phi
843          call ylmcd(il-1,im-il,kcart,dth,dphi)
844          dylmcart(:) = dth*gradth(:) + dphi*gradphi(:)
845 
846          ! Cartesian to crystallographic axis
847          ! Notice: a bug was discovered by Marco Cazzaniga, december 2009
848          ! the transformation matrix A=(a1,a2,a3) must act *on its left* on the
849          ! covariant vector dylmcart (a *row* vector). The previous implementation assumed A
850          ! acting on its right on a column vector, yielding wrong results for the (small)
851          ! non local contributions to the spectra, such as a spurious anisotropy in isotropic systems.
852          ! This is the correct version:
853          dylmcrys(1) = (a1(1)*dylmcart(1)+a1(2)*dylmcart(2)+a1(3)*dylmcart(3))/(two_pi)
854          dylmcrys(2) = (a2(1)*dylmcart(1)+a2(2)*dylmcart(2)+a2(3)*dylmcart(3))/(two_pi)
855          dylmcrys(3) = (a3(1)*dylmcart(1)+a3(2)*dylmcart(2)+a3(3)*dylmcart(3))/(two_pi)
856 
857          ! Note that fnld is simple precision complex, it could be possible to use double precision
858          do ii=1,3
859            fnld(ii,ig,ilm,in,iat) = factor*sfac* &
860             ( kg(ii)/mkg*ylmc(il-1,im-il,kcart)*vkbd(ig,iln,itypat) + dylmcrys(ii)*vkb(ig,iln,itypat) )
861          end do
862 
863        end do !im
864      end do !il
865    end do !iat
866  end do !ig
867 
868  DBG_EXIT("COLL")
869 
870 end subroutine ccgradvnl_ylm

m_vkbr/nc_ihr_comm [ Functions ]

[ Top ] [ m_vkbr ] [ Functions ]

NAME

  nc_pwihr_comm

FUNCTION

  Calculate the matrix elements of the commutator i[H,r]
  For norm conserving potentials 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.

PARENTS

CHILDREN

SOURCE

617 function nc_ihr_comm(vkbr,cryst,psps,npw,nspinor,istwfk,inclvkb,kpoint,ug1,ug2,gvec) result(ihr_comm)
618 
619 
620 !This section has been created automatically by the script Abilint (TD).
621 !Do not modify the following lines by hand.
622 #undef ABI_FUNC
623 #define ABI_FUNC 'nc_ihr_comm'
624 !End of the abilint section
625 
626  implicit none
627 
628 !Arguments ------------------------------------
629 !scalars
630  integer,intent(in) :: npw,nspinor,inclvkb,istwfk
631  type(vkbr_t),intent(in) :: vkbr
632  type(crystal_t),intent(in) :: cryst
633  type(Pseudopotential_type),intent(in) :: psps
634 !arrays
635  integer,intent(in) :: gvec(3,npw)
636  real(dp),intent(in) :: kpoint(3)
637  complex(gwpc),intent(in) :: ug1(npw*nspinor),ug2(npw*nspinor)
638  complex(gwpc) :: ihr_comm(3,nspinor**2)
639 
640 !Local variables ------------------------------
641 !scalars
642  integer :: ig,iab,spad1,spad2
643  complex(dpc) :: c_tmp
644 !arrays
645  integer :: spinorwf_pad(2,4)
646 
647 !************************************************************************
648 
649  ! [H, r] = -\nabla + [V_{nl}, r]
650  ! V_nl is present only in the case of NC pseudos but
651  ! not in PAW unless even the AE Hamiltonian in non-local e.g. LDA+U or LEXX.
652 
653  ! -i <c,k|\nabla_r|v,k> in reduced coordinates is always included.
654  ! -i <c,k|\nabla_r|v,k> = \sum_G u_{ck}^*(G) [k+G] u_{vk}(G)
655  ! Note that here we assume c/=v, moreover the ug are supposed to be orthonormal and
656  ! hence k+G can be replaced by G.
657  ! HM 03/08/2018: we need band velocities so we don't assume c/=v anymore and we use
658  ! k+G.
659 
660  spinorwf_pad = RESHAPE([0, 0, npw, npw, 0, npw, npw, 0], [2, 4])
661  ihr_comm = czero
662 
663  ! -i <c,k|\nabla_r|v,k> in reduced coordinates.
664  ! This term is spin diagonal if nspinor == 2
665  if (istwfk == 1) then
666    do iab=1,nspinor
667      spad1 = spinorwf_pad(1,iab); spad2 = spinorwf_pad(2,iab)
668      do ig=1,npw
669        c_tmp = CONJG(ug1(ig+spad1)) * ug2(ig+spad2)
670        ihr_comm(:,iab) = ihr_comm(:,iab) + c_tmp * (kpoint + gvec(:,ig))
671      end do
672    end do
673  else
674    ! Symmetrized expression: \sum_G  (k+G) 2i Im [ u_a^*(G) u_b(G) ]. (k0,G0) term is null.
675    ABI_CHECK(nspinor == 1, "nspinor != 1")
676    do ig=1,npw
677      c_tmp = CONJG(ug1(ig)) * ug2(ig)
678      ihr_comm(:,1) = ihr_comm(:,1) + two*j_dpc * AIMAG(c_tmp) * (kpoint + gvec(:,ig))
679    end do
680  end if
681 
682  ! Add second term $i <c,k|[Vnl,r]|v,k> $ in reduced cordinates.
683  if (inclvkb /= 0) then
684    ABI_CHECK(istwfk == vkbr%istwfk, "input istwfk /= vkbr%istwfk")
685    call add_vnlr_commutator(vkbr,cryst,psps,npw,nspinor,ug1,ug2,ihr_comm)
686  end if
687 
688 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

PARENTS

      m_vkbr

CHILDREN

      ylmcd

SOURCE

225 subroutine vkbr_free_0D(vkbr)
226 
227 
228 !This section has been created automatically by the script Abilint (TD).
229 !Do not modify the following lines by hand.
230 #undef ABI_FUNC
231 #define ABI_FUNC 'vkbr_free_0D'
232 !End of the abilint section
233 
234  implicit none
235 
236 !Arguments ------------------------------------
237 !scalars
238  type(vkbr_t),intent(inout) :: vkbr
239 
240 !************************************************************************
241 
242 !complex
243  if (allocated(vkbr%fnl)) then
244    ABI_FREE(vkbr%fnl)
245  end if
246  if (allocated(vkbr%fnld)) then
247    ABI_FREE(vkbr%fnld)
248  end if
249 
250 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

PARENTS

CHILDREN

      ylmcd

SOURCE

269 subroutine vkbr_free_1D(vkbr)
270 
271 
272 !This section has been created automatically by the script Abilint (TD).
273 !Do not modify the following lines by hand.
274 #undef ABI_FUNC
275 #define ABI_FUNC 'vkbr_free_1D'
276 !End of the abilint section
277 
278  implicit none
279 
280 !Arguments ------------------------------------
281 !arrays
282  type(vkbr_t),intent(inout) :: vkbr(:)
283 
284 !Local variables ------------------------------
285 !scalars
286  integer :: ii
287 
288 !************************************************************************
289 
290  do ii=1,SIZE(vkbr)
291    call vkbr_free_0D(vkbr(ii))
292  end do
293 
294 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.

PARENTS

      calc_optical_mels,cchi0q0,cchi0q0_intraband

CHILDREN

      ylmcd

SOURCE

133 subroutine vkbr_init(vkbr,cryst,psps,inclvkb,istwfk,npw,kpoint,gvec)
134 
135 
136 !This section has been created automatically by the script Abilint (TD).
137 !Do not modify the following lines by hand.
138 #undef ABI_FUNC
139 #define ABI_FUNC 'vkbr_init'
140 !End of the abilint section
141 
142  implicit none
143 
144 !Arguments ------------------------------------
145 !scalars
146  integer,intent(in) :: npw,inclvkb,istwfk
147  type(crystal_t),intent(in) :: cryst
148  type(vkbr_t),intent(inout) :: vkbr
149  type(pseudopotential_type),intent(in) :: psps
150 !arrays
151  integer,intent(in) :: gvec(3,npw)
152  real(dp),intent(in) :: kpoint(3)
153 
154 !Local variables-------------------------------
155 !scalars
156  integer :: ierr
157  character(len=500) :: msg
158 !arrays
159  real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:)
160 
161 !************************************************************************
162 
163  !@vkbr_t
164  vkbr%istwfk = istwfk
165  vkbr%ntypat = cryst%ntypat
166  vkbr%natom = cryst%natom
167  vkbr%mpsang  = psps%mpsang
168  vkbr%npw = npw
169  vkbr%inclvkb = inclvkb
170  vkbr%kpoint = kpoint
171 
172  ! Calculate KB form factors and derivatives.
173  ! The arrays are allocated with lnmax to support pseudos with more than projector.
174  ! Note that lnmax takes into account lloc hence arrays are in packed form and one should be
175  ! accessed with the indices provided by indlmn.
176  ! TODO: they should be calculated on-the-fly using calc_vkb
177  !       For the moment, we opt for a quick an dirty implementation.
178 
179  ABI_MALLOC(vkbsign, (psps%lnmax, cryst%ntypat))
180  ABI_MALLOC(vkb, (npw, psps%lnmax, cryst%ntypat))
181  ABI_MALLOC(vkbd, (npw, psps%lnmax, cryst%ntypat))
182  call calc_vkb(cryst,psps,kpoint,npw,npw,gvec,vkbsign,vkb,vkbd)
183 
184  select case (inclvkb)
185  case (2)
186    ! Complex spherical harmonics (CPU and mem \propto npw).
187    write(msg,'(a,f12.1)')'out-of-memory in fnl; Mb= ',one*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwpc*b2Mb
188    ABI_STAT_MALLOC(vkbr%fnl,(npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr)
189    ABI_CHECK(ierr==0, msg)
190 
191    write(msg,'(a,f12.1)')'out-of-memory in fnld; Mb= ',three*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwpc*b2Mb
192    ABI_STAT_MALLOC(vkbr%fnld,(3,npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr)
193    ABI_CHECK(ierr==0, msg)
194 
195    call ccgradvnl_ylm(cryst,psps,npw,gvec,kpoint,vkbsign,vkb,vkbd,vkbr%fnl,vkbr%fnld)
196 
197  case default
198    MSG_ERROR(sjoin("Wrong inclvkb= ",itoa(inclvkb)))
199  end select
200 
201  ABI_FREE(vkbsign)
202  ABI_FREE(vkb)
203  ABI_FREE(vkbd)
204 
205 end subroutine vkbr_init

m_vkbr/vkbr_t [ Types ]

[ Top ] [ m_vkbr ] [ 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

58  type,public :: vkbr_t
59 
60   integer :: istwfk
61   ! Storage mode of the G vectors for this k-point.
62 
63   integer :: ntypat
64   ! Number of type of atoms
65 
66   integer :: natom
67   ! Number of atoms
68 
69   integer :: mpsang
70   ! Max l+1 over atoms
71 
72   integer :: npw
73   ! Number of G-vectors.
74 
75   integer :: inclvkb
76   ! Option for calculating the matrix elements of [Vnl,r].
77 
78   real(dp) :: kpoint(3)
79   ! The k-point in reduced coordinates.
80 
81   complex(gwpc),allocatable :: fnl(:,:,:,:)
82   ! fnl(npw,mpsang**2,mproj,natom)
83 
84   complex(gwpc),allocatable :: fnld(:,:,:,:,:)
85   ! fnld(3,npw,mpsang**2,mproj,natom)
86 
87  end type vkbr_t
88 
89  public :: vkbr_init       ! vkbr_t Constructor
90  public :: vkbr_free       ! Free memory
91  public :: nc_ihr_comm     ! Compute matrix elements of the commutator i[H,r] for NC pseudos
92  public :: calc_vkb        ! Kleynman-Bylander form factors and derivatives.