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. LDA+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-2018 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

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

PARENTS

      cchi0q0

CHILDREN

      simp_gen

SOURCE

351 subroutine paw_cross_ihr_comm(ihr_comm,nspinor,nr,Cryst,Pawfgrtab,Paw_onsite,&
352 & ur_ae1,ur_ae2,ur_ae_onsite1,ur_ae_onsite2,Cprj_kb1,Cprj_kb2)
353 
354 
355 !This section has been created automatically by the script Abilint (TD).
356 !Do not modify the following lines by hand.
357 #undef ABI_FUNC
358 #define ABI_FUNC 'paw_cross_ihr_comm'
359 !End of the abilint section
360 
361  implicit none
362 
363 !Arguments ------------------------------------
364 !scalars
365  integer,intent(in) :: nspinor,nr
366  type(crystal_t),intent(in) :: Cryst
367 !arrays
368  complex(gwpc),intent(inout) :: ihr_comm(3,nspinor**2)
369  complex(gwpc),intent(in) :: ur_ae1(nr),ur_ae2(nr)
370  complex(gwpc),intent(in) :: ur_ae_onsite1(nr),ur_ae_onsite2(nr)
371  type(pawfgrtab_type),intent(in) :: Pawfgrtab(Cryst%natom)
372  type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
373  type(pawcprj_type),intent(in) :: Cprj_kb1(Cryst%natom,nspinor),Cprj_kb2(Cryst%natom,nspinor)
374 
375 !Local variables-------------------------------
376  integer :: iatom,lmn_size,ilmn,ifgd,ifftsph,nfgd
377  complex(dpc) :: cp1, cp2
378  complex(dpc) :: cross1,cross2
379 !arrays
380  real(dp) :: gspace_cart2red(3,3)
381  complex(gwpc) :: ihr_comm_cart(3,nspinor**2)
382  complex(dpc) :: dphigr(3), dphigr1(3),dphigr2(3)
383 
384 ! *************************************************************************
385 
386  ABI_CHECK(nspinor==1,"nspinor + pawcross not implemented")
387 
388  ! [H, r] = -\nabla + [V_{nl}, r]
389  ! The V_nl part, present in case of LDA+U, is omitted for the cross terms contribution
390  ! Recall that delta_rho_tw_ij = (psi_i - phi_i)* (phi_j - tphi_j) + (phi_i - tphi_i)* (psi_j - phi_j)
391  ihr_comm_cart(:,1) = czero
392 
393  do iatom=1,Cryst%natom
394    lmn_size = Paw_onsite(iatom)%lmn_size
395    nfgd = Pawfgrtab(iatom)%nfgd
396 
397    do ifgd=1,nfgd
398 
399      ifftsph = Pawfgrtab(iatom)%ifftsph(ifgd)
400 
401      cross1 = ur_ae1(ifftsph) - ur_ae_onsite1(ifftsph)
402      cross2 = ur_ae2(ifftsph) - ur_ae_onsite2(ifftsph)
403 
404      do ilmn=1,lmn_size
405 
406        dphigr(1:3) = Paw_onsite(iatom)%phi_gr(1:3,ifgd,ilmn) - Paw_onsite(iatom)%tphi_gr(1:3,ifgd,ilmn)
407 
408        cp1 = CMPLX(Cprj_kb1(iatom,1)%cp(1,ilmn),Cprj_kb1(iatom,1)%cp(2,ilmn)) * sqrt(Cryst%ucvol) ! that damn magic factor
409        cp2 = CMPLX(Cprj_kb2(iatom,1)%cp(1,ilmn),Cprj_kb2(iatom,1)%cp(2,ilmn)) * sqrt(Cryst%ucvol)
410 
411        dphigr1(1:3) = cp1 * dphigr(1:3)
412        dphigr2(1:3) = cp2 * dphigr(1:3)
413 
414        ihr_comm_cart(1,1) = ihr_comm_cart(1,1) - j_dpc * (CONJG(cross1) * dphigr2(1) - CONJG(dphigr1(1)) * cross2) / nr
415        ihr_comm_cart(2,1) = ihr_comm_cart(2,1) - j_dpc * (CONJG(cross1) * dphigr2(2) - CONJG(dphigr1(2)) * cross2) / nr
416        ihr_comm_cart(3,1) = ihr_comm_cart(3,1) - j_dpc * (CONJG(cross1) * dphigr2(3) - CONJG(dphigr1(3)) * cross2) / nr
417 
418      end do
419    end do
420  end do
421 
422  ! Go to reduced coordinate
423  gspace_cart2red = TRANSPOSE(Cryst%rprimd)
424  ihr_comm(:,1) = ihr_comm(:,1) +  MATMUL(gspace_cart2red, ihr_comm_cart(:,1)) / two_pi
425 
426 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 LDA+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>$

PARENTS

      cchi0q0,debug_tools,spectra

SOURCE

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

m_paw_hr/pawhur_free [ Functions ]

[ Top ] [ m_paw_hr ] [ Functions ]

NAME

 pawhur_free

FUNCTION

  Deallocate memory 

PARENTS

      bethe_salpeter,cchi0q0,cchi0q0_intraband

CHILDREN

      simp_gen

SOURCE

110 subroutine pawhur_free(Hur)
111 
112 
113 !This section has been created automatically by the script Abilint (TD).
114 !Do not modify the following lines by hand.
115 #undef ABI_FUNC
116 #define ABI_FUNC 'pawhur_free'
117 !End of the abilint section
118 
119  implicit none
120 
121 !Arguments ------------------------------------
122  type(pawhur_t),intent(inout) :: Hur(:)
123 
124 !Local variables-------------------------------
125  integer :: iat
126 ! *************************************************************************
127 
128  do iat=1,SIZE(Hur)
129    if (allocated(Hur(iat)%ij_select)) then
130      ABI_FREE(Hur(iat)%ij_select)
131    end if
132    if (allocated(Hur(iat)%commutator)) then
133      ABI_FREE(Hur(iat)%commutator)
134    end if
135  end do
136 
137 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

PARENTS

      bethe_salpeter,cchi0q0,cchi0q0_intraband

CHILDREN

      simp_gen

SOURCE

450 subroutine pawhur_init(hur,nsppol,pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij)
451 
452 
453 !This section has been created automatically by the script Abilint (TD).
454 !Do not modify the following lines by hand.
455 #undef ABI_FUNC
456 #define ABI_FUNC 'pawhur_init'
457 !End of the abilint section
458 
459  implicit none
460 
461 !Arguments ------------------------------------
462 !scalars
463  integer,intent(in) :: nsppol,pawprtvol
464  type(crystal_t),intent(in) :: Cryst
465  type(Pawang_type),intent(in) :: Pawang
466 !arrays
467  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
468  type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)
469  type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom)
470  type(pawhur_t),intent(inout) :: Hur(Cryst%natom)
471 
472 !Local variables-------------------------------
473 !scalars 
474  integer :: iatom,ij_idx,isel,itypat,isppol,lmn2_size_max,lmn2_size,lmn_size,lpawu
475  integer :: jlmn,jl,jm,jlm,jln,k0lmn,k0lm,k0ln,ilmn,il,im,ilm,iln
476  integer :: m2,m1,left_lmn,right_lmn,tot_lmn,nmax
477 !arrays
478  integer :: nsel(3,nsppol)
479  integer, ABI_CONTIGUOUS pointer :: indlmn(:,:)
480  real(dp) :: sumr_ij(3)
481  real(dp),allocatable :: rcart_onsite(:,:,:)
482  real(dp),allocatable :: rij_tmp(:,:,:),vpawu(:,:,:,:)
483 
484 ! *************************************************************************
485 
486  ! Get onsite matrix elements of the position operator.
487  lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size) 
488  ABI_MALLOC(rcart_onsite,(3,lmn2_size_max,Cryst%natom))
489 
490  call pawr(Pawtab,Pawrad,Pawang,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,lmn2_size_max,rcart_onsite)
491 
492  do iatom=1,Cryst%natom
493    itypat=Cryst%typat(iatom)
494    if (Pawtab(itypat)%usepawu==0) CYCLE
495    lmn2_size=Pawtab(itypat)%lmn2_size
496    lmn_size =Pawtab(itypat)%lmn_size
497    lpawu=Pawtab(itypat)%lpawu
498    Hur(iatom)%lmn2_size=lmn2_size
499    Hur(iatom)%lmn_size =lmn_size
500    Hur(iatom)%nsppol   =nsppol
501    indlmn => Pawtab(itypat)%indlmn
502 
503    ABI_MALLOC(rij_tmp,(3,lmn_size**2,nsppol))
504    rij_tmp=zero
505 
506    ! Get Vpawu^{\sigma}_{m1,m2}
507    ABI_MALLOC(vpawu,(Paw_ij(iatom)%cplex_dij,2*lpawu+1,2*lpawu+1,Paw_ij(iatom)%ndij))
508    call pawpupot(Paw_ij(iatom)%cplex_dij,Paw_ij(iatom)%ndij,&
509 &                Paw_ij(iatom)%noccmmp,Paw_ij(iatom)%nocctot,&
510 &                pawprtvol,Pawtab(itypat),vpawu)
511 
512    do isppol=1,nsppol ! spinor not implemented
513 
514      ! === Loop on (jl,jm,jn) channels ===
515      ij_idx=0
516      do jlmn=1,lmn_size
517        jl =indlmn(1,jlmn)
518        jm =indlmn(2,jlmn)
519        jlm=indlmn(4,jlmn)
520        jln=indlmn(5,jlmn)
521 
522        k0lmn=jlmn*(jlmn-1)/2 
523        k0lm =jlm *(jlm -1)/2
524        k0ln =jln *(jln -1)/2
525        !
526        ! === Loop on (il,im,in) channels === 
527        ! * Looping over all ij components. Elements are not symmetric.
528        do ilmn=1,lmn_size
529          il =indlmn(1,ilmn)
530          im =indlmn(2,ilmn)
531          ilm=indlmn(4,ilmn)
532          iln=indlmn(5,ilmn)
533 
534          ij_idx=ij_idx+1
535 
536          ! === Selection rules ===
537          if (il/=lpawu.and.jl/=lpawu) CYCLE 
538 
539          sumr_ij(:)=zero 
540          do m2=1,2*lpawu+1
541            do m1=1,2*lpawu+1
542              if (m1==(im-lpawu-1).and.il==lpawu) then 
543                left_lmn =ilmn-(il+im+1)+m2
544                right_lmn=jlmn
545                if (right_lmn>=left_lmn) then 
546                  tot_lmn=right_lmn*(right_lmn-1)/2 + left_lmn
547                else 
548                  tot_lmn=left_lmn*(left_lmn-1)/2 + right_lmn
549                end if
550                sumr_ij=sumr_ij+vpawu(1,m1,m2,isppol)*rcart_onsite(:,tot_lmn,iatom)
551              end if
552 
553              if (m2==(jm-lpawu-1).and.jl==lpawu) then 
554                left_lmn =ilmn
555                right_lmn=jlmn-(jl+jm+1)+m1
556                if (right_lmn>=left_lmn) then 
557                  tot_lmn=right_lmn*(right_lmn-1)/2 + left_lmn
558                else 
559                  tot_lmn=left_lmn*(left_lmn-1)/2 + right_lmn
560                end if
561                sumr_ij=sumr_ij+vpawu(1,m1,m2,isppol)*rcart_onsite(:,tot_lmn,iatom)
562              end if
563            end do !m1
564          end do !m2
565 
566          rij_tmp(:,ij_idx,isppol)=sumr_ij(:)
567 
568        end do !ilmn
569      end do !jlmn
570    end do !isppol
571 
572    ABI_FREE(vpawu)
573 
574    ! === Save values in packed form ===
575    ABI_MALLOC(Hur(iatom)%ij_select,(lmn_size,lmn_size,nsppol))
576    Hur(iatom)%ij_select=0
577    nsel(:,:)=COUNT(ABS(rij_tmp)>tol6,DIM=2)
578    nmax=MAXVAL(nsel)
579    ABI_MALLOC(Hur(iatom)%commutator,(3,nmax,nsppol))
580    do isppol=1,nsppol
581      ij_idx=0
582      isel  =0
583      do jlmn=1,lmn_size
584        do ilmn=1,lmn_size
585          ij_idx=ij_idx+1
586          if (ANY (ABS(rij_tmp(:,ij_idx,isppol))>tol6) ) then
587            isel=isel+1
588            Hur(iatom)%ij_select(ilmn,jlmn,isppol)=isel
589            Hur(iatom)%commutator(:,isel,isppol)=rij_tmp(:,ij_idx,isppol)
590          end if
591        end do
592      end do
593    end do
594 
595    ABI_FREE(rij_tmp)
596  end do !iatom
597 
598  ABI_FREE(rcart_onsite)
599 
600 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 LDA+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

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

PARENTS

      m_paw_hr

CHILDREN

      simp_gen

SOURCE

641 subroutine pawr(Pawtab,Pawrad,Pawang,natom,ntypat,typat,xcart,lmn2_size_max,rcart_onsite)
642 
643 
644 !This section has been created automatically by the script Abilint (TD).
645 !Do not modify the following lines by hand.
646 #undef ABI_FUNC
647 #define ABI_FUNC 'pawr'
648 !End of the abilint section
649 
650  implicit none
651 
652 !Arguments ------------------------------------
653 !scalars
654  integer,intent(in) :: lmn2_size_max,natom,ntypat
655  type(Pawang_type),intent(in) :: Pawang
656 
657 !arrays
658  integer,intent(in) :: typat(natom)
659  real(dp),intent(in) :: xcart(3,natom)
660  real(dp),intent(inout) :: rcart_onsite(3,lmn2_size_max,natom)
661  type(Pawrad_type),intent(in) :: Pawrad(ntypat)
662  type(Pawtab_type),target,intent(in) :: Pawtab(ntypat)
663 
664 !Local variables-------------------------------
665 !scalars
666  integer,parameter :: ll1=1
667  integer :: iatom,idir,ignt,il,ilm,ilm_G,ilmn,iln,im,itypat,jl,jlm,jlmn,jln,jm,k0lm
668  integer :: k0lmn,k0ln,klm,klmn,kln,lmn_size,mesh_size,mm_G,lmn2_size
669  real(dp) :: fact,intff,rgnt
670 !arrays
671  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
672  real(dp),allocatable :: ff(:),rad(:),rc_tmp(:,:)
673 
674 ! *************************************************************************
675 
676  DBG_ENTER("COLL")
677 
678  fact=two*SQRT(pi/three)
679  rcart_onsite(:,:,:)=zero
680 
681  do itypat=1,ntypat
682    lmn_size  =Pawtab(itypat)%lmn_size
683    lmn2_size =Pawtab(itypat)%lmn2_size
684    mesh_size =Pawtab(itypat)%mesh_size
685    indlmn =>  Pawtab(itypat)%indlmn
686 
687    ABI_MALLOC(ff,(mesh_size))
688    ABI_MALLOC(rad,(mesh_size))
689    rad(1:mesh_size)=Pawrad(itypat)%rad(1:mesh_size)
690 
691    ABI_MALLOC(rc_tmp,(3,lmn2_size))
692    rc_tmp=zero
693    !
694    ! === Loop on (jl,jm,jn) channels
695    do jlmn=1,lmn_size
696      jl =indlmn(1,jlmn)
697      jm =indlmn(2,jlmn)
698      jlm=indlmn(4,jlmn)
699      jln=indlmn(5,jlmn)
700 
701      k0lmn=jlmn*(jlmn-1)/2
702      k0lm =jlm *(jlm -1)/2
703      k0ln =jln *(jln -1)/2
704      !
705      ! === Loop on (il,im,in) channels; klmn is the index for packed form ===
706      do ilmn=1,jlmn
707        il =indlmn(1,ilmn)
708        im =indlmn(2,ilmn)
709        ilm=indlmn(4,ilmn)
710        iln=indlmn(5,ilmn)
711 
712        klmn=k0lmn+ilmn
713        klm =k0lm +ilm
714        kln =k0ln +iln
715        !
716        ! === For each cartesian direction, use expansion in terms of RSH ===
717        ! TODO Add a check if l=1 is in the set
718        do idir=1,3
719          mm_G=0
720          if (idir==1) mm_G= 1
721          if (idir==2) mm_G=-1
722          if (idir==3) mm_G= 0
723          ilm_G=1+ll1**2+ll1+mm_G
724          ignt=Pawang%gntselect(ilm_G,klm)
725          if (ignt/=0) then
726            rgnt=Pawang%realgnt(ignt)
727            ff(1)=zero
728            !ff(2:mesh_size)=(Pawtab(itypat)%phiphj(2:mesh_size,kln)-Pawtab(itypat)%tphitphj(2:mesh_size,kln))*rad(2:mesh_size)
729            ff(2:mesh_size)=Pawtab(itypat)%phiphj(2:mesh_size,kln)*rad(2:mesh_size)
730            call simp_gen(intff,ff,Pawrad(itypat))
731            rc_tmp(idir,klmn)=fact*intff*rgnt
732          end if
733        end do !idir
734 
735      end do !ilmn
736    end do !jllmn
737 
738    ! === Make matrix elements for each atom of this type ===
739    do jlmn=1,lmn_size
740      jl =indlmn(1,jlmn)
741      jm =indlmn(2,jlmn)
742      jln=indlmn(5,jlmn)
743 
744      k0lmn=jlmn*(jlmn-1)/2
745      k0ln =jln *(jln -1)/2
746      do ilmn=1,jlmn
747        il =indlmn(1,ilmn)
748        im =indlmn(2,ilmn)
749        iln=indlmn(5,ilmn)
750 
751        klmn=k0lmn+ilmn
752        kln =k0ln +iln
753 
754        intff=zero
755        if (il==jl.and.jm==im) then
756          ff(1:mesh_size)=Pawtab(itypat)%phiphj(1:mesh_size,kln)
757          call simp_gen(intff,ff,Pawrad(itypat))
758        end if
759        do iatom=1,natom
760          if (typat(iatom)/=itypat) CYCLE
761          rcart_onsite(:,klmn,iatom)=rc_tmp(:,klmn) + xcart(:,iatom)*intff
762        end do
763 
764      end do ! ilmn
765    end do !jlmn
766 
767    ABI_FREE(ff)
768    ABI_FREE(rad)
769    ABI_FREE(rc_tmp)
770  end do !itypat
771 
772  DBG_EXIT("COLL")
773 
774 end subroutine pawr