TABLE OF CONTENTS
- ABINIT/m_paw_hr
- m_paw_hr/paw_cross_ihr_comm
- m_paw_hr/paw_ihr
- m_paw_hr/pawhur_free
- m_paw_hr/pawhur_init
- m_paw_hr/pawhur_t
- m_paw_hr/pawr
ABINIT/m_paw_hr [ 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