TABLE OF CONTENTS


ABINIT/m_paw_hr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_hr

FUNCTION

  This module provides objects and methods to calculate the matrix elements 
  of the commutator PAW [H,r] needed for the correct treatment of the optical limit q-->0
  in the matrix elements <k-q,b1|e^{-iqr}|k,b2>. As PAW is a full potential method 
  the commutator reduces to the contribution given by the velocity operator. 
  However, when the all-electron Hamiltonian is non-local (e.g. DFT+U or 
  LEXX) additional on-site terms have to be considered in the calculation of the
  matrix elements of [H.r].

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

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_paw_hr
29 
30  use defs_basis
31  use m_abicore
32  use m_errors
33 
34  use m_crystal,        only : crystal_t
35  use m_pawang,         only : pawang_type
36  use m_pawrad,         only : pawrad_type, simp_gen
37  use m_pawtab,         only : pawtab_type
38  use m_paw_ij,         only : paw_ij_type
39  use m_pawfgrtab,      only : pawfgrtab_type
40  use m_pawcprj,        only : pawcprj_type
41  use m_pawdij,         only : pawpupot
42  use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t
43  
44  implicit none
45 
46  private

m_paw_hr/paw_cross_ihr_comm [ Functions ]

[ Top ] [ m_paw_hr ] [ Functions ]

NAME

 paw_cross_ihr_comm

FUNCTION

  Adds the PAW cross term contribution to the matrix elements of the  i\nabla operator.
  in cartesian coordinates. Should take also into account the contribution arising from the U 
  part of the Hamiltonian (if any)

INPUTS

  ihr_comm = the commutator [H,r] evaluated between states i and j, with only the plane-wave and 
              the onsite parts included
  isppol=Spin index.
  nspinor=Number of spinori components.
  nr=Number real-space points on the fine fft grid for the ae wavefunctions
  kpoint(3)=k-point in reduced coordinates.
  Cryst<crystal_t>=Info on the crystal structure.
    %natom=Number of atoms in unit cell
    %typat(natom)
  Pawfgrtab(ntypat)= PAW tabulated data on the fine grid
    %lmn_size Number of (l,m,n) elements for the paw basis
    %nfgr Number of points on the fine grid
    %ifftsph Indexes of the fine-grid points on the fft mesh
  Paw_onsite(ntypat)= PAW tabulated data on the fine grid points inside the sphere
    %phi_gr(3,nfgr,lmn_size) gradient of phi in cartesian coordinates
    %tphi_gr(3,nfgr,lmn_size) gradient of tphi in cartesian coordinates
  ur_ae1(nr),ur_ae2(nr)=Left and right AE wavefunction.
  ur_ae_onsite1(nr),ur_ae_onsite2(nr)=Left and right AE onsite wavefunction.
  Cprj_kb1(natom,nspinor),Cprj_kb2(natom,nspinor) <type(pawcprj_type)>=
   projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to 
   wavefunctions (k,b1,s) and (k,b2,s), respectively.

OUTPUT

SIDE EFFECTS

  The cross-term contribution is added to the commutator 

SOURCE

321 subroutine paw_cross_ihr_comm(ihr_comm,nspinor,nr,Cryst,Pawfgrtab,Paw_onsite,&
322 & ur_ae1,ur_ae2,ur_ae_onsite1,ur_ae_onsite2,Cprj_kb1,Cprj_kb2)
323 
324  implicit none
325 
326 !Arguments ------------------------------------
327 !scalars
328  integer,intent(in) :: nspinor,nr
329  type(crystal_t),intent(in) :: Cryst
330 !arrays
331  complex(gwpc),intent(inout) :: ihr_comm(3,nspinor**2)
332  complex(gwpc),intent(in) :: ur_ae1(nr),ur_ae2(nr)
333  complex(gwpc),intent(in) :: ur_ae_onsite1(nr),ur_ae_onsite2(nr)
334  type(pawfgrtab_type),intent(in) :: Pawfgrtab(Cryst%natom)
335  type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
336  type(pawcprj_type),intent(in) :: Cprj_kb1(Cryst%natom,nspinor),Cprj_kb2(Cryst%natom,nspinor)
337 
338 !Local variables-------------------------------
339  integer :: iatom,lmn_size,ilmn,ifgd,ifftsph,nfgd
340  complex(dpc) :: cp1, cp2
341  complex(dpc) :: cross1,cross2
342 !arrays
343  real(dp) :: gspace_cart2red(3,3)
344  complex(gwpc) :: ihr_comm_cart(3,nspinor**2)
345  complex(dpc) :: dphigr(3), dphigr1(3),dphigr2(3)
346 
347 ! *************************************************************************
348 
349  ABI_CHECK(nspinor==1,"nspinor + pawcross not implemented")
350 
351  ! [H, r] = -\nabla + [V_{nl}, r]
352  ! The V_nl part, present in case of DFT+U, is omitted for the cross terms contribution
353  ! Recall that delta_rho_tw_ij = (psi_i - phi_i)* (phi_j - tphi_j) + (phi_i - tphi_i)* (psi_j - phi_j)
354  ihr_comm_cart(:,1) = czero
355 
356  do iatom=1,Cryst%natom
357    lmn_size = Paw_onsite(iatom)%lmn_size
358    nfgd = Pawfgrtab(iatom)%nfgd
359 
360    do ifgd=1,nfgd
361 
362      ifftsph = Pawfgrtab(iatom)%ifftsph(ifgd)
363 
364      cross1 = ur_ae1(ifftsph) - ur_ae_onsite1(ifftsph)
365      cross2 = ur_ae2(ifftsph) - ur_ae_onsite2(ifftsph)
366 
367      do ilmn=1,lmn_size
368 
369        dphigr(1:3) = Paw_onsite(iatom)%phi_gr(1:3,ifgd,ilmn) - Paw_onsite(iatom)%tphi_gr(1:3,ifgd,ilmn)
370 
371        cp1 = CMPLX(Cprj_kb1(iatom,1)%cp(1,ilmn),Cprj_kb1(iatom,1)%cp(2,ilmn)) * sqrt(Cryst%ucvol) ! that damn magic factor
372        cp2 = CMPLX(Cprj_kb2(iatom,1)%cp(1,ilmn),Cprj_kb2(iatom,1)%cp(2,ilmn)) * sqrt(Cryst%ucvol)
373 
374        dphigr1(1:3) = cp1 * dphigr(1:3)
375        dphigr2(1:3) = cp2 * dphigr(1:3)
376 
377        ihr_comm_cart(1,1) = ihr_comm_cart(1,1) - j_dpc * (CONJG(cross1) * dphigr2(1) - CONJG(dphigr1(1)) * cross2) / nr
378        ihr_comm_cart(2,1) = ihr_comm_cart(2,1) - j_dpc * (CONJG(cross1) * dphigr2(2) - CONJG(dphigr1(2)) * cross2) / nr
379        ihr_comm_cart(3,1) = ihr_comm_cart(3,1) - j_dpc * (CONJG(cross1) * dphigr2(3) - CONJG(dphigr1(3)) * cross2) / nr
380 
381      end do
382    end do
383  end do
384 
385  ! Go to reduced coordinate
386  gspace_cart2red = TRANSPOSE(Cryst%rprimd)
387  ihr_comm(:,1) = ihr_comm(:,1) +  MATMUL(gspace_cart2red, ihr_comm_cart(:,1)) / two_pi
388 
389 end subroutine paw_cross_ihr_comm

m_paw_hr/paw_ihr [ Functions ]

[ Top ] [ m_paw_hr ] [ Functions ]

NAME

 paw_ihr

FUNCTION

  Calculate the PAW onsite contribution to the matrix elements of the i\nabla operator.
  in cartesian coordinates. Take also into account the contribution arising from the U 
  part of the Hamiltonian (if any)

INPUTS

  isppol=Spin index.
  nspinor=Number of spinori components.
  npw=Number of planewaves for this k-point.
  istwfk=Storage mode for the wavefunctions.
  kpoint(3)=k-point in reduced coordinates.
  Cryst<crystal_t>=Info on the crystal structure.
    %natom=Number of atoms in unit cell
    %typat(natom)
  Pawtab(ntypat)=Only for PAW, TABulated data initialized at start
    %lmn_size Number of (l,m,n) elements for the paw basis
    %nabla_ij(3,lmn_size,lmn_size)) Onsite contribution
      <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for each type
  ug1(nspinor*npwwfn)=Left wavefunction.
  ug2(nspinor*npwwfn)=Right wavefunction
  HUr(natom)=Commutator of the DFT+U part of the Hamiltonian with the position operator.
  Cprj_kb1(natom,nspinor),Cprj_kb2(natom,nspinor) <type(pawcprj_type)>=
   projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to 
   wavefunctions (k,b1,s) and (k,b2,s), respectively.

OUTPUT

  onsite(2,3)=Onsite contribution to  $i<ug1|\nabla|ug2>$

SOURCE

162 function paw_ihr(isppol,nspinor,npw,istwfk,kpoint,Cryst,Pawtab,ug1,ug2,gvec,Cprj_kb1,Cprj_kb2,HUr) result(ihr_comm)
163 
164  implicit none
165 
166 !Arguments ------------------------------------
167 !scalars
168  integer,intent(in) :: isppol,nspinor,npw,istwfk
169  complex(gwpc) :: ihr_comm(3,nspinor**2)
170  type(crystal_t),intent(in) :: Cryst
171 !arrays
172  integer,intent(in) :: gvec(3,npw)
173  real(dp),intent(in) :: kpoint(3)
174  complex(gwpc),intent(in) :: ug1(nspinor*npw),ug2(nspinor*npw)
175  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
176  type(pawcprj_type),intent(in) :: Cprj_kb1(Cryst%natom,nspinor),Cprj_kb2(Cryst%natom,nspinor)
177  type(pawhur_t),intent(in) :: Hur(Cryst%natom)
178 
179 !Local variables-------------------------------
180  integer :: iatom,itypat,lmn_size,ilmn,jlmn,isel
181  integer :: ig,iab,spad1,spad2
182  real(dp) :: re_p,im_p
183  complex(dpc) :: ctemp
184 !arrays
185  integer :: spinorwf_pad(2,4)
186  real(dp) :: hurc_ij(3),ons_cart(2,3) !,ons_comm_red(2,3)
187  real(dp) :: gspace_cart2red(3,3) !rs_cart2red(3,3),
188  real(dp), ABI_CONTIGUOUS pointer :: nabla_ij(:,:,:)
189  complex(gwpc) :: ihr_comm_cart(3,nspinor**2)
190 
191 ! *************************************************************************
192 
193  ! [H, r] = -\nabla + [V_{nl}, r] 
194  ! Note that V_nl is present only if the AE-Hamiltonian is non-local e.g. DFT+U or LEXX.
195  spinorwf_pad=RESHAPE((/0,0,npw,npw,0,npw,npw,0/),(/2,4/))
196  ihr_comm=zero
197 
198  ! -i <c,k|\nabla_r|v,k> = \sum_G u_{ck}^*(G) [k+G] u_{vk}(G) in reduced coordinates.
199  if (istwfk==1) then
200    do iab=1,nspinor**2
201      spad1 = spinorwf_pad(1,iab)
202      spad2 = spinorwf_pad(2,iab)
203      do ig=1,npw 
204        ctemp = CONJG(ug1(ig+spad1)) * ug2(ig+spad2)
205        ihr_comm(:,iab) = ihr_comm(:,iab) + ctemp* ( kpoint + gvec(:,ig))
206      end do
207    end do
208  else 
209    ! Symmetrized expression: \sum_G  (k+G) 2i Im [ u_a^*(G) u_b(G) ]. (k0,G0) term is null.
210    do ig=1,npw 
211      ctemp = CONJG(ug1(ig)) * ug2(ig)
212      ihr_comm(:,1) = ihr_comm(:,1) + two*j_dpc * AIMAG(ctemp) * (kpoint + gvec(:,ig))
213    end do
214  end if
215  !
216  ! Add on-site terms.
217  ons_cart=zero
218  ABI_CHECK(nspinor==1,"nspinor/=1 not coded")
219 
220  do iatom=1,Cryst%natom 
221    itypat=Cryst%typat(iatom)
222    lmn_size=Pawtab(itypat)%lmn_size
223    nabla_ij => Pawtab(itypat)%nabla_ij(:,:,:) 
224    !
225    !=== Unpacked loop over lmn channels ====
226    do jlmn=1,lmn_size
227      do ilmn=1,lmn_size
228        re_p =  Cprj_kb1(iatom,1)%cp(1,ilmn)*Cprj_kb2(iatom,1)%cp(1,jlmn) &
229 &             +Cprj_kb1(iatom,1)%cp(2,ilmn)*Cprj_kb2(iatom,1)%cp(2,jlmn) 
230  
231        im_p =  Cprj_kb1(iatom,1)%cp(1,ilmn)*Cprj_kb2(iatom,1)%cp(2,jlmn) &
232 &             -Cprj_kb1(iatom,1)%cp(2,ilmn)*Cprj_kb2(iatom,1)%cp(1,jlmn)
233 
234        ! Onsite contribution given by -i\nabla.
235        ons_cart(1,1)=ons_cart(1,1) + im_p*nabla_ij(1,ilmn,jlmn)
236        ons_cart(1,2)=ons_cart(1,2) + im_p*nabla_ij(2,ilmn,jlmn)
237        ons_cart(1,3)=ons_cart(1,3) + im_p*nabla_ij(3,ilmn,jlmn)
238 
239        ons_cart(2,1)=ons_cart(2,1) - re_p*nabla_ij(1,ilmn,jlmn)
240        ons_cart(2,2)=ons_cart(2,2) - re_p*nabla_ij(2,ilmn,jlmn)
241        ons_cart(2,3)=ons_cart(2,3) - re_p*nabla_ij(3,ilmn,jlmn)
242        !
243        if (Pawtab(itypat)%usepawu/=0) then ! Add i[V_u, r] 
244          isel=Hur(iatom)%ij_select(ilmn,jlmn,isppol)
245          if (isel>0) then
246            hurc_ij(:)=Hur(iatom)%commutator(:,isel,isppol)
247 
248            ons_cart(1,1)=ons_cart(1,1) - im_p*hurc_ij(1)
249            ons_cart(1,2)=ons_cart(1,2) - im_p*hurc_ij(2)
250            ons_cart(1,3)=ons_cart(1,3) - im_p*hurc_ij(3)
251 
252            ons_cart(2,1)=ons_cart(2,1) + re_p*hurc_ij(1)
253            ons_cart(2,2)=ons_cart(2,2) + re_p*hurc_ij(2)
254            ons_cart(2,3)=ons_cart(2,3) + re_p*hurc_ij(3)
255          end if
256        end if
257 
258      end do !ilmn
259    end do !jlmn
260  end do !iatom
261 
262  ! ons_cart is in Cartesian coordinates in real space 
263  ! while ihr_comm is in reduced coordinates in reciprocal space in terms of gprimd.
264  !rs_cart2red = TRANSPOSE(Cryst%gprimd) ! if <r> is in terms of real space vectors
265  gspace_cart2red = TRANSPOSE(Cryst%rprimd)
266 
267  !ons_comm_red(1,:)=MATMUL(rs_cart2red,ons_comm(1,:))
268  !ons_comm_red(2,:)=MATMUL(rs_cart2red,ons_comm(2,:))
269  !ihr_comm(:,1) = ihr_comm(:,1) + CMPLX(ons_comm_red(1,:),ons_comm_red(2,:),kind=gwpc)
270 
271  ihr_comm_cart(:,1) = two_pi*MATMUL(Cryst%gprimd,ihr_comm(:,1))
272  ihr_comm_cart(:,1) = ihr_comm_cart(:,1) + CMPLX(ons_cart(1,:),ons_cart(2,:),kind=gwpc)
273 
274  ! Final result is in reduced coordinates, in terms of gprimd.
275  ihr_comm(:,1) = MATMUL(gspace_cart2red, ihr_comm_cart(:,1))/two_pi
276 
277 end function paw_ihr

m_paw_hr/pawhur_free [ Functions ]

[ Top ] [ m_paw_hr ] [ Functions ]

NAME

 pawhur_free

FUNCTION

  Deallocate memory 

SOURCE

103 subroutine pawhur_free(Hur)
104 
105  implicit none
106 
107 !Arguments ------------------------------------
108  type(pawhur_t),intent(inout) :: Hur(:)
109 
110 !Local variables-------------------------------
111  integer :: iat
112 ! *************************************************************************
113 
114  do iat=1,SIZE(Hur)
115    if (allocated(Hur(iat)%ij_select)) then
116      ABI_FREE(Hur(iat)%ij_select)
117    end if
118    if (allocated(Hur(iat)%commutator)) then
119      ABI_FREE(Hur(iat)%commutator)
120    end if
121  end do
122 
123 end subroutine pawhur_free

m_paw_hr/pawhur_init [ Functions ]

[ Top ] [ m_paw_hr ] [ Functions ]

NAME

 pawhur_init

FUNCTION

  Creation method for the pawhur_t data type.

INPUTS

OUTPUT

SOURCE

407 subroutine pawhur_init(hur,nsppol,pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij)
408 
409  implicit none
410 
411 !Arguments ------------------------------------
412 !scalars
413  integer,intent(in) :: nsppol,pawprtvol
414  type(crystal_t),intent(in) :: Cryst
415  type(Pawang_type),intent(in) :: Pawang
416 !arrays
417  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
418  type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)
419  type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom)
420  type(pawhur_t),intent(inout) :: Hur(Cryst%natom)
421 
422 !Local variables-------------------------------
423 !scalars 
424  integer :: iatom,ij_idx,isel,itypat,isppol,lmn2_size_max,lmn2_size,lmn_size,lpawu
425  integer :: jlmn,jl,jm,jlm,jln,k0lmn,k0lm,k0ln,ilmn,il,im,ilm,iln
426  integer :: m2,m1,left_lmn,right_lmn,tot_lmn,nmax
427 !arrays
428  integer :: nsel(3,nsppol)
429  integer, ABI_CONTIGUOUS pointer :: indlmn(:,:)
430  real(dp) :: sumr_ij(3)
431  real(dp),allocatable :: rcart_onsite(:,:,:)
432  real(dp),allocatable :: rij_tmp(:,:,:),vpawu(:,:,:,:)
433 
434 ! *************************************************************************
435 
436  ! Get onsite matrix elements of the position operator.
437  lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size) 
438  ABI_MALLOC(rcart_onsite,(3,lmn2_size_max,Cryst%natom))
439 
440  call pawr(Pawtab,Pawrad,Pawang,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,lmn2_size_max,rcart_onsite)
441 
442  do iatom=1,Cryst%natom
443    itypat=Cryst%typat(iatom)
444    if (Pawtab(itypat)%usepawu==0) CYCLE
445    lmn2_size=Pawtab(itypat)%lmn2_size
446    lmn_size =Pawtab(itypat)%lmn_size
447    lpawu=Pawtab(itypat)%lpawu
448    Hur(iatom)%lmn2_size=lmn2_size
449    Hur(iatom)%lmn_size =lmn_size
450    Hur(iatom)%nsppol   =nsppol
451    indlmn => Pawtab(itypat)%indlmn
452 
453    ABI_MALLOC(rij_tmp,(3,lmn_size**2,nsppol))
454    rij_tmp=zero
455 
456    ! Get Vpawu^{\sigma}_{m1,m2}
457    ABI_MALLOC(vpawu,(Paw_ij(iatom)%cplex_dij,2*lpawu+1,2*lpawu+1,Paw_ij(iatom)%ndij))
458    call pawpupot(Paw_ij(iatom)%cplex_dij,Paw_ij(iatom)%ndij,&
459 &                Paw_ij(iatom)%noccmmp,Paw_ij(iatom)%nocctot,&
460 &                pawprtvol,Pawtab(itypat),vpawu)
461 
462    do isppol=1,nsppol ! spinor not implemented
463 
464      ! === Loop on (jl,jm,jn) channels ===
465      ij_idx=0
466      do jlmn=1,lmn_size
467        jl =indlmn(1,jlmn)
468        jm =indlmn(2,jlmn)
469        jlm=indlmn(4,jlmn)
470        jln=indlmn(5,jlmn)
471 
472        k0lmn=jlmn*(jlmn-1)/2 
473        k0lm =jlm *(jlm -1)/2
474        k0ln =jln *(jln -1)/2
475        !
476        ! === Loop on (il,im,in) channels === 
477        ! * Looping over all ij components. Elements are not symmetric.
478        do ilmn=1,lmn_size
479          il =indlmn(1,ilmn)
480          im =indlmn(2,ilmn)
481          ilm=indlmn(4,ilmn)
482          iln=indlmn(5,ilmn)
483 
484          ij_idx=ij_idx+1
485 
486          ! === Selection rules ===
487          if (il/=lpawu.and.jl/=lpawu) CYCLE 
488 
489          sumr_ij(:)=zero 
490          do m2=1,2*lpawu+1
491            do m1=1,2*lpawu+1
492              if (m1==(im-lpawu-1).and.il==lpawu) then 
493                left_lmn =ilmn-(il+im+1)+m2
494                right_lmn=jlmn
495                if (right_lmn>=left_lmn) then 
496                  tot_lmn=right_lmn*(right_lmn-1)/2 + left_lmn
497                else 
498                  tot_lmn=left_lmn*(left_lmn-1)/2 + right_lmn
499                end if
500                sumr_ij=sumr_ij+vpawu(1,m1,m2,isppol)*rcart_onsite(:,tot_lmn,iatom)
501              end if
502 
503              if (m2==(jm-lpawu-1).and.jl==lpawu) then 
504                left_lmn =ilmn
505                right_lmn=jlmn-(jl+jm+1)+m1
506                if (right_lmn>=left_lmn) then 
507                  tot_lmn=right_lmn*(right_lmn-1)/2 + left_lmn
508                else 
509                  tot_lmn=left_lmn*(left_lmn-1)/2 + right_lmn
510                end if
511                sumr_ij=sumr_ij+vpawu(1,m1,m2,isppol)*rcart_onsite(:,tot_lmn,iatom)
512              end if
513            end do !m1
514          end do !m2
515 
516          rij_tmp(:,ij_idx,isppol)=sumr_ij(:)
517 
518        end do !ilmn
519      end do !jlmn
520    end do !isppol
521 
522    ABI_FREE(vpawu)
523 
524    ! === Save values in packed form ===
525    ABI_MALLOC(Hur(iatom)%ij_select,(lmn_size,lmn_size,nsppol))
526    Hur(iatom)%ij_select=0
527    nsel(:,:)=COUNT(ABS(rij_tmp)>tol6,DIM=2)
528    nmax=MAXVAL(nsel)
529    ABI_MALLOC(Hur(iatom)%commutator,(3,nmax,nsppol))
530    do isppol=1,nsppol
531      ij_idx=0
532      isel  =0
533      do jlmn=1,lmn_size
534        do ilmn=1,lmn_size
535          ij_idx=ij_idx+1
536          if (ANY (ABS(rij_tmp(:,ij_idx,isppol))>tol6) ) then
537            isel=isel+1
538            Hur(iatom)%ij_select(ilmn,jlmn,isppol)=isel
539            Hur(iatom)%commutator(:,isel,isppol)=rij_tmp(:,ij_idx,isppol)
540          end if
541        end do
542      end do
543    end do
544 
545    ABI_FREE(rij_tmp)
546  end do !iatom
547 
548  ABI_FREE(rcart_onsite)
549 
550 end subroutine pawhur_init

m_paw_hr/pawhur_t [ Types ]

[ Top ] [ m_paw_hr ] [ Types ]

NAME

  pawhur_t

FUNCTION

  The pawhur_t data type stores basic dimensions and quantities 
  used in the GW part for the treatment of the non-analytic behavior of the 
  heads and wings of the irreducible polarizability in the long wave-length limit (i.e. q-->0).
  Note that, within the PAW formalism, a standard KS Hamiltonian has a semi-local contribution 
  arising from the kinetic operator (if we work in the AE representation). 
  When DFT+U is used, a fully non-local term is added to the Hamiltonian whose commutator with the position operator 
  has to be considered during the calculation of the heads and wings of the polarizability in the optical limit

SOURCE

66  type,public :: pawhur_t
67 
68   integer :: lmn_size
69   integer :: lmn2_size
70   integer :: nsppol
71   !integer :: nsel
72 
73   integer,allocatable :: ij_select(:,:,:)
74   ! ijselect(lmn_size,lmn_size,nsppol) 
75   ! Selection rules of ij matrix elements
76   ! Do not take into account selection on x-y-x for the time being.
77 
78   real(dp),allocatable :: commutator(:,:,:)
79   ! commutator(3,nsel,nsppol)
80  end type pawhur_t
81 
82  public ::  pawhur_init          ! Init object
83  public ::  pawhur_free          ! Deallocate memory
84  public ::  paw_ihr              
85  public ::  paw_cross_ihr_comm

m_paw_hr/pawr [ Functions ]

[ Top ] [ m_paw_hr ] [ Functions ]

NAME

 pawr

FUNCTION

 Evaluate matrix elements of the position operator between PAW AE partial waves.

INPUTS

  Pawtab(ntypat) <type(pawtab_type)>=paw tabulated data read at start:
     %lmn_size
     %lmn2_size
     %indklmn
     %phiphj
  Pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
     %mesh_size=Dimension of radial mesh
     %rad(mesh_size)=The coordinates of all the points of the radial mesh
  Pawang <type(pawang_type)>=paw angular mesh and related data
     %lmax=Maximum value of angular momentum l+1
     %gntselect((2*l_max-1)**2,l_max**2,l_max**2)= selection rules for Gaunt coefficients
     %realgnt
  natom=number of atoms in unit cell
  ntypat=number of types of atom
  typat(natom)=type of each atom
  xcart(3,natom)=cartesian coordinates

OUTPUT

  rcart_onsite(3,lmn2_size_max,natom)

SOURCE

585 subroutine pawr(Pawtab,Pawrad,Pawang,natom,ntypat,typat,xcart,lmn2_size_max,rcart_onsite)
586 
587  implicit none
588 
589 !Arguments ------------------------------------
590 !scalars
591  integer,intent(in) :: lmn2_size_max,natom,ntypat
592  type(Pawang_type),intent(in) :: Pawang
593 
594 !arrays
595  integer,intent(in) :: typat(natom)
596  real(dp),intent(in) :: xcart(3,natom)
597  real(dp),intent(inout) :: rcart_onsite(3,lmn2_size_max,natom)
598  type(Pawrad_type),intent(in) :: Pawrad(ntypat)
599  type(Pawtab_type),target,intent(in) :: Pawtab(ntypat)
600 
601 !Local variables-------------------------------
602 !scalars
603  integer,parameter :: ll1=1
604  integer :: iatom,idir,ignt,il,ilm,ilm_G,ilmn,iln,im,itypat,jl,jlm,jlmn,jln,jm,k0lm
605  integer :: k0lmn,k0ln,klm,klmn,kln,lmn_size,mesh_size,mm_G,lmn2_size
606  real(dp) :: fact,intff,rgnt
607 !arrays
608  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
609  real(dp),allocatable :: ff(:),rad(:),rc_tmp(:,:)
610 
611 ! *************************************************************************
612 
613  DBG_ENTER("COLL")
614 
615  fact=two*SQRT(pi/three)
616  rcart_onsite(:,:,:)=zero
617 
618  do itypat=1,ntypat
619    lmn_size  =Pawtab(itypat)%lmn_size
620    lmn2_size =Pawtab(itypat)%lmn2_size
621    mesh_size =Pawtab(itypat)%mesh_size
622    indlmn =>  Pawtab(itypat)%indlmn
623 
624    ABI_MALLOC(ff,(mesh_size))
625    ABI_MALLOC(rad,(mesh_size))
626    rad(1:mesh_size)=Pawrad(itypat)%rad(1:mesh_size)
627 
628    ABI_MALLOC(rc_tmp,(3,lmn2_size))
629    rc_tmp=zero
630    !
631    ! === Loop on (jl,jm,jn) channels
632    do jlmn=1,lmn_size
633      jl =indlmn(1,jlmn)
634      jm =indlmn(2,jlmn)
635      jlm=indlmn(4,jlmn)
636      jln=indlmn(5,jlmn)
637 
638      k0lmn=jlmn*(jlmn-1)/2
639      k0lm =jlm *(jlm -1)/2
640      k0ln =jln *(jln -1)/2
641      !
642      ! === Loop on (il,im,in) channels; klmn is the index for packed form ===
643      do ilmn=1,jlmn
644        il =indlmn(1,ilmn)
645        im =indlmn(2,ilmn)
646        ilm=indlmn(4,ilmn)
647        iln=indlmn(5,ilmn)
648 
649        klmn=k0lmn+ilmn
650        klm =k0lm +ilm
651        kln =k0ln +iln
652        !
653        ! === For each cartesian direction, use expansion in terms of RSH ===
654        ! TODO Add a check if l=1 is in the set
655        do idir=1,3
656          mm_G=0
657          if (idir==1) mm_G= 1
658          if (idir==2) mm_G=-1
659          if (idir==3) mm_G= 0
660          ilm_G=1+ll1**2+ll1+mm_G
661          ignt=Pawang%gntselect(ilm_G,klm)
662          if (ignt/=0) then
663            rgnt=Pawang%realgnt(ignt)
664            ff(1)=zero
665            !ff(2:mesh_size)=(Pawtab(itypat)%phiphj(2:mesh_size,kln)-Pawtab(itypat)%tphitphj(2:mesh_size,kln))*rad(2:mesh_size)
666            ff(2:mesh_size)=Pawtab(itypat)%phiphj(2:mesh_size,kln)*rad(2:mesh_size)
667            call simp_gen(intff,ff,Pawrad(itypat))
668            rc_tmp(idir,klmn)=fact*intff*rgnt
669          end if
670        end do !idir
671 
672      end do !ilmn
673    end do !jllmn
674 
675    ! === Make matrix elements for each atom of this type ===
676    do jlmn=1,lmn_size
677      jl =indlmn(1,jlmn)
678      jm =indlmn(2,jlmn)
679      jln=indlmn(5,jlmn)
680 
681      k0lmn=jlmn*(jlmn-1)/2
682      k0ln =jln *(jln -1)/2
683      do ilmn=1,jlmn
684        il =indlmn(1,ilmn)
685        im =indlmn(2,ilmn)
686        iln=indlmn(5,ilmn)
687 
688        klmn=k0lmn+ilmn
689        kln =k0ln +iln
690 
691        intff=zero
692        if (il==jl.and.jm==im) then
693          ff(1:mesh_size)=Pawtab(itypat)%phiphj(1:mesh_size,kln)
694          call simp_gen(intff,ff,Pawrad(itypat))
695        end if
696        do iatom=1,natom
697          if (typat(iatom)/=itypat) CYCLE
698          rcart_onsite(:,klmn,iatom)=rc_tmp(:,klmn) + xcart(:,iatom)*intff
699        end do
700 
701      end do ! ilmn
702    end do !jlmn
703 
704    ABI_FREE(ff)
705    ABI_FREE(rad)
706    ABI_FREE(rc_tmp)
707  end do !itypat
708 
709  DBG_EXIT("COLL")
710 
711 end subroutine pawr