TABLE OF CONTENTS


ABINIT/calc_optical_mels [ Functions ]

[ Top ] [ Functions ]

NAME

  calc_optical_mels

FUNCTION

  Calculate all optical matrix elements in the BZ.

INPUTS

 lomo_spin(Wfd%nsppol)=Index of the lomo band for the different spins.
 lomo_min,max_band=minimum and max band index to be calculated.
 nkbz=Number of points in the full Brillouin zone.
 inclvkb=if different from 0, [Vnl,r] is included in the calculation of the
   matrix element of the velocity operator. No meaning for PAW (except for DFT+U)
 qpt(3)
 Kmesh<kmesh_t>=Info on the k-point sampling for wave functions.
 Cryst<crystal_t>=Structure defining the crystalline structure.
 KS_Bst<ebands_t>
 Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data
 Psps <pseudopotential_type>=variables related to pseudopotentials.
 Hur(Cryst%natom*usepaw)<pawhur_t>=Only for PAW and DFT+U, quantities used to evaluate the commutator [H_u,r].
 Wfd<wfdgw_t>=Handler for the wavefunctions.

OUTPUT

 opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol)=Matrix elements <c k|e^{+iqr}|v k>

SOURCE

 77 subroutine calc_optical_mels(Wfd,Kmesh,KS_Bst,Cryst,Psps,Pawtab,Hur,&
 78 &  inclvkb,lomo_spin,lomo_min,max_band,nkbz,qpoint,opt_cvk)
 79 
 80 !Arguments ------------------------------------
 81 !scalars
 82  integer,intent(in) :: nkbz,inclvkb,lomo_min,max_band
 83  type(kmesh_t),intent(in) :: Kmesh
 84  type(crystal_t),intent(in) :: Cryst
 85  type(pseudopotential_type),intent(in) :: Psps
 86  type(ebands_t),intent(in) :: KS_Bst
 87  type(wfdgw_t),target,intent(inout) :: Wfd
 88 !arrays
 89  integer,intent(in) :: lomo_spin(Wfd%nsppol)
 90  real(dp),intent(in) :: qpoint(3)
 91  complex(dpc),intent(out) :: opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,Wfd%nsppol)
 92  type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
 93  type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
 94 
 95 !Local variables ------------------------------
 96 !scalars
 97  integer :: nsppol,usepaw,nspinor,comm,spin,npw_k,istwf_k,my_nbbp
 98  integer :: ik_bz,ik_ibz,itim_k,isym_k,ib_c,ib_v,ierr,my_rank
 99  real(dp) :: ediff
100  complex(dpc) :: emcvk
101  character(len=500) :: msg
102  type(vkbr_t) :: vkbr
103  type(wave_t),pointer :: wave_v, wave_c
104 !arrays
105  integer,allocatable :: bbp_distrb(:,:)
106  integer,ABI_CONTIGUOUS pointer :: kg_k(:,:)
107  real(dp) :: mat_dp(3,3),qrot(3),b1(3),b2(3),b3(3),kbz(3)
108  complex(dpc),allocatable :: ir_kibz(:,:,:,:,:)
109  complex(gwpc), ABI_CONTIGUOUS pointer :: ug_c(:),ug_v(:)
110  complex(gwpc) :: ihrc(3,Wfd%nspinor**2)
111  logical :: bbp_mask(Wfd%mband,Wfd%mband)
112  type(pawcprj_type),allocatable :: Cp_v(:,:),Cp_c(:,:)
113 
114 !************************************************************************
115 
116  call wrtout(std_out," Calculating optical matrix elements in the IBZ","COLL")
117  ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
118 
119  comm = Wfd%comm
120  my_rank = Wfd%my_rank
121 
122  nsppol  = Wfd%nsppol
123  nspinor = Wfd%nspinor
124  usepaw  = Wfd%usepaw
125 
126  if (usepaw==1) then
127    ABI_MALLOC(Cp_v,(Wfd%natom,nspinor))
128    call pawcprj_alloc(Cp_v,0,Wfd%nlmn_atm)
129    ABI_MALLOC(Cp_c,(Wfd%natom,nspinor))
130    call pawcprj_alloc(Cp_c,0,Wfd%nlmn_atm)
131  end if
132 
133  if (inclvkb==1.and.usepaw==0) then
134    ABI_ERROR("inclvkb==1 not coded,using inclvkb==2")
135  end if
136  !
137  ! Calculate the matrix elements of ir in the IBZ.
138  ABI_MALLOC(ir_kibz,(3,lomo_min:max_band,lomo_min:max_band,Wfd%nkibz,nsppol))
139  ir_kibz=czero
140 
141  ABI_MALLOC(bbp_distrb, (Wfd%mband,Wfd%mband))
142 
143  do spin=1,nsppol
144    do ik_ibz=1,Wfd%nkibz
145     !
146     ! Distribute the (b,b') entries.
147     bbp_mask=.FALSE.; bbp_mask(lomo_spin(spin):max_band,lomo_spin(spin):max_band)=.TRUE.
148     call wfd%distribute_bbp(ik_ibz,spin,"All",my_nbbp,bbp_distrb,bbp_mask=bbp_mask)
149     if (ALL(bbp_distrb/=my_rank)) CYCLE
150 
151     istwf_k = Wfd%istwfk(ik_ibz)
152     ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded") ! KB stuff is missing.
153     npw_k = Wfd%npwarr(ik_ibz)
154     kg_k  => Wfd%Kdata(ik_ibz)%kg_k
155 
156     if (inclvkb/=0.and.usepaw==0) then
157       ! Prepare term i <n,k|[Vnl,r]|n"k>
158       call vkbr_init(vkbr,Cryst,Psps,inclvkb,istwf_k,npw_k,Kmesh%ibz(:,ik_ibz),kg_k)
159     end if
160 
161     ! Note: spinorial case is not coded therefore we work with ihrc(:,1).
162     ! TODO: The lower triangle can be Reconstructed by symmetry.
163     do ib_v=lomo_spin(spin),max_band ! Loop over bands
164       if ( ALL(bbp_distrb(ib_v,:)/=my_rank) ) CYCLE
165 
166       ABI_CHECK(wfd%get_wave_ptr(ib_v, ik_ibz, spin, wave_v, msg) == 0, msg)
167       ug_v => wave_v%ug
168       if (usepaw==1) call wfd%get_cprj(ib_v,ik_ibz,spin,Cryst,Cp_v,sorted=.FALSE.)
169 
170       do ib_c=lomo_spin(spin),max_band
171        if (bbp_distrb(ib_v,ib_c)/=my_rank) CYCLE
172        ABI_CHECK(wfd%get_wave_ptr(ib_c, ik_ibz, spin, wave_c, msg) == 0, msg)
173        ug_c => wave_c%ug
174 
175        if (usepaw==0) then
176          ! Calculate matrix elements of i[H,r] for NC pseudopotentials.
177          ihrc = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,Kmesh%ibz(:,ik_ibz),ug_c,ug_v,kg_k)
178 
179        else
180          ! Matrix elements of i[H,r] for PAW.
181          call wfd%get_cprj(ib_c,ik_ibz,spin,Cryst,Cp_c,sorted=.FALSE.)
182 
183          ihrc = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug_c,ug_v,kg_k,Cp_c,Cp_v,HUr)
184        end if
185        !
186        ! Save matrix elements of i*r in the IBZ
187        ediff = KS_Bst%eig(ib_c,ik_ibz,spin) - KS_BSt%eig(ib_v,ik_ibz,spin)
188        if (ABS(ediff)<tol16) ediff=tol6  ! Treat a possible degeneracy between v and c.
189        ir_kibz(:,ib_c,ib_v,ik_ibz,spin) = ihrc(:,1)/ediff
190 
191       end do !ib_c
192     end do !ib_v
193 
194     call vkbr_free(vkbr)
195    end do !spin
196  end do !ik_ibz
197 
198  ! Collect results on each node.
199  call xmpi_sum(ir_kibz,comm,ierr)
200 
201  ABI_FREE(bbp_distrb)
202 
203  if (usepaw==1) then
204    call pawcprj_free(Cp_v)
205    ABI_FREE(Cp_v)
206    call pawcprj_free(Cp_c)
207    ABI_FREE(Cp_c)
208  end if
209  !
210  ! ======================================================
211  ! ==== Calculate Fcv(kBZ) in the full Brilouin zone ====
212  ! ======================================================
213  !
214  ! Symmetrization of the matrix elements of the position operator.
215  ! <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
216  !   where S is one of the symrec operations in reciprocal space, R is the
217  !   corresponding operation in real space, \tau being the associated fractional translations.
218  !
219  ! q.Mcv( Sk) =  S^{-1}q. Mcv(k)
220  ! q.Mcv(-Sk) = -S^{-1}q. CONJG(Mcv(k)) if time-reversal is used.
221 
222  b1=Cryst%gprimd(:,1)*two_pi
223  b2=Cryst%gprimd(:,2)*two_pi
224  b3=Cryst%gprimd(:,3)*two_pi
225 
226  opt_cvk = czero
227  do spin=1,nsppol
228    do ik_bz=1,nkbz
229     !
230     ! Get ik_ibz, and symmetries index from ik_bz.
231     call Kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k)
232 
233     mat_dp = DBLE(Cryst%symrec(:,:,isym_k))
234     call matrginv(mat_dp,3,3) ! Invert
235     qrot = (3-2*itim_k) * MATMUL(mat_dp,qpoint)
236 
237     do ib_v=lomo_spin(spin),max_band !  Loops over the bands C and V start
238       do ib_c=lomo_spin(spin),max_band
239         !if (ib_c==ib_v) CYCLE
240         emcvk = pdtqrc(qrot,ir_kibz(:,ib_c,ib_v,ik_ibz,spin),b1,b2,b3)
241         if (itim_k==2) emcvk = CONJG(emcvk)
242         opt_cvk(ib_c,ib_v,ik_bz,spin) = emcvk
243       end do !ib_c
244     end do !ib_v
245 
246    end do !ik_bz
247  end do !spin
248 
249  ABI_FREE(ir_kibz)
250 
251  call xmpi_barrier(comm)
252 
253 contains

ABINIT/m_wfd_optic [ Modules ]

[ Top ] [ Modules ]

NAME

  m_wfd_optic

FUNCTION

  Functions to compute optical matrix elements using the wavefunction descriptor.

COPYRIGHT

  Copyright (C) 2008-2022 ABINIT group (MG)
  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

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_wfd_optic
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27  use m_xmpi
28 
29  use defs_datatypes,      only : ebands_t, pseudopotential_type
30  use m_hide_lapack,       only : matrginv
31  use m_bz_mesh,           only : kmesh_t
32  use m_crystal,           only : crystal_t
33  use m_vkbr,              only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm
34  use m_wfd,               only : wfdgw_t, wave_t
35  use m_pawtab,            only : pawtab_type
36  use m_pawcprj,           only : pawcprj_type, pawcprj_alloc, pawcprj_free
37  use m_paw_hr,            only : pawhur_t, paw_ihr
38 
39  implicit none
40 
41  private

ABINIT/pdtqrc [ Functions ]

[ Top ] [ Functions ]

NAME

  pdtqrc

FUNCTION

  Calculate the dot product of a real vector with a complex vector, where each is in terms of b1-b3

INPUTS

OUTPUT

SOURCE

269 pure function pdtqrc(R,C,b1,b2,b3)
270 
271 !Arguments ------------------------------------
272 !arrays
273  real(dp),intent(in) :: R(3),b1(3),b2(3),b3(3)
274  complex(dpc),intent(in) :: C(3)
275  complex(dpc) :: pdtqrc
276 
277 !Local variables ------------------------------
278 !scalars
279  integer :: ii
280 
281 !************************************************************************
282 
283  pdtqrc=czero
284  do ii=1,3
285    pdtqrc = pdtqrc + (R(1)*b1(ii)+R(2)*b2(ii)+R(3)*b3(ii)) * &
286 &                    (C(1)*b1(ii)+C(2)*b2(ii)+C(3)*b3(ii))
287  end do
288 
289 end function pdtqrc