Line data Source code
1 : !!****m* ABINIT/m_paw_hr
2 : !! NAME
3 : !! m_paw_hr
4 : !!
5 : !! FUNCTION
6 : !! This module provides objects and methods to calculate the matrix elements
7 : !! of the commutator PAW [H,r] needed for the correct treatment of the optical limit q-->0
8 : !! in the matrix elements <k-q,b1|e^{-iqr}|k,b2>. As PAW is a full potential method
9 : !! the commutator reduces to the contribution given by the velocity operator.
10 : !! However, when the all-electron Hamiltonian is non-local (e.g. DFT+U or
11 : !! LEXX) additional on-site terms have to be considered in the calculation of the
12 : !! matrix elements of [H.r].
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 2008-2026 ABINIT group (MG)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !!
20 : !! SOURCE
21 :
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
47 : !!***
48 :
49 : !----------------------------------------------------------------------
50 :
51 : !!****t* m_paw_hr/pawhur_t
52 : !! NAME
53 : !! pawhur_t
54 : !!
55 : !! FUNCTION
56 : !! The pawhur_t data type stores basic dimensions and quantities
57 : !! used in the GW part for the treatment of the non-analytic behavior of the
58 : !! heads and wings of the irreducible polarizability in the long wave-length limit (i.e. q-->0).
59 : !! Note that, within the PAW formalism, a standard KS Hamiltonian has a semi-local contribution
60 : !! arising from the kinetic operator (if we work in the AE representation).
61 : !! When DFT+U is used, a fully non-local term is added to the Hamiltonian whose commutator with the position operator
62 : !! has to be considered during the calculation of the heads and wings of the polarizability in the optical limit
63 : !!
64 : !! SOURCE
65 :
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
86 :
87 : !!***
88 :
89 : CONTAINS !========================================================================================
90 : !!***
91 :
92 : !----------------------------------------------------------------------
93 :
94 : !!****f* m_paw_hr/pawhur_free
95 : !! NAME
96 : !! pawhur_free
97 : !!
98 : !! FUNCTION
99 : !! Deallocate memory
100 : !!
101 : !! SOURCE
102 :
103 34 : subroutine pawhur_free(Hur)
104 :
105 : !Arguments ------------------------------------
106 : type(pawhur_t),intent(inout) :: Hur(:)
107 :
108 : !Local variables-------------------------------
109 : integer :: iat
110 : ! *************************************************************************
111 :
112 69 : do iat=1,SIZE(Hur)
113 35 : ABI_SFREE(Hur(iat)%ij_select)
114 69 : ABI_SFREE(Hur(iat)%commutator)
115 : end do
116 :
117 34 : end subroutine pawhur_free
118 : !!***
119 :
120 : !----------------------------------------------------------------------
121 :
122 : !!****f* m_paw_hr/paw_ihr
123 : !! NAME
124 : !! paw_ihr
125 : !!
126 : !! FUNCTION
127 : !! Calculate the PAW onsite contribution to the matrix elements of the i\nabla operator.
128 : !! in cartesian coordinates. Take also into account the contribution arising from the U
129 : !! part of the Hamiltonian (if any)
130 : !!
131 : !! INPUTS
132 : !! isppol=Spin index.
133 : !! nspinor=Number of spinori components.
134 : !! npw=Number of planewaves for this k-point.
135 : !! istwfk=Storage mode for the wavefunctions.
136 : !! kpoint(3)=k-point in reduced coordinates.
137 : !! Cryst<crystal_t>=Info on the crystal structure.
138 : !! %natom=Number of atoms in unit cell
139 : !! %typat(natom)
140 : !! Pawtab(ntypat)=Only for PAW, TABulated data initialized at start
141 : !! %lmn_size Number of (l,m,n) elements for the paw basis
142 : !! %nabla_ij(3,lmn_size,lmn_size)) Onsite contribution
143 : !! <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for each type
144 : !! ug1(nspinor*npwwfn)=Left wavefunction.
145 : !! ug2(nspinor*npwwfn)=Right wavefunction
146 : !! HUr(natom)=Commutator of the DFT+U part of the Hamiltonian with the position operator.
147 : !! Cprj_kb1(natom,nspinor),Cprj_kb2(natom,nspinor) <type(pawcprj_type)>=
148 : !! projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to
149 : !! wavefunctions (k,b1,s) and (k,b2,s), respectively.
150 : !!
151 : !! OUTPUT
152 : !! onsite(2,3)=Onsite contribution to $i<ug1|\nabla|ug2>$
153 : !!
154 : !! SOURCE
155 :
156 3623 : function paw_ihr(isppol,nspinor,npw,istwfk,kpoint,Cryst,Pawtab,ug1,ug2,gvec,Cprj_kb1,Cprj_kb2,HUr) result(ihr_comm)
157 :
158 : !Arguments ------------------------------------
159 : !scalars
160 : integer,intent(in) :: isppol,nspinor,npw,istwfk
161 : complex(gwp) :: ihr_comm(3,nspinor**2)
162 : type(crystal_t),intent(in) :: Cryst
163 : !arrays
164 : integer,intent(in) :: gvec(3,npw)
165 : real(dp),intent(in) :: kpoint(3)
166 : complex(gwp),intent(in) :: ug1(nspinor*npw),ug2(nspinor*npw)
167 : type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
168 : type(pawcprj_type),intent(in) :: Cprj_kb1(Cryst%natom,nspinor),Cprj_kb2(Cryst%natom,nspinor)
169 : type(pawhur_t),intent(in) :: Hur(Cryst%natom)
170 :
171 : !Local variables-------------------------------
172 : integer :: iatom,itypat,lmn_size,ilmn,jlmn,isel
173 : integer :: ig,iab,spad1,spad2
174 : real(dp) :: re_p,im_p
175 : complex(dp) :: ctemp
176 : !arrays
177 : integer :: spinorwf_pad(2,4)
178 : real(dp) :: hurc_ij(3),ons_cart(2,3) !,ons_comm_red(2,3)
179 : real(dp) :: gspace_cart2red(3,3) !rs_cart2red(3,3),
180 3623 : real(dp), ABI_CONTIGUOUS pointer :: nabla_ij(:,:,:)
181 7246 : complex(gwp) :: ihr_comm_cart(3,nspinor**2)
182 : ! *************************************************************************
183 :
184 : ! [H, r] = -\nabla + [V_{nl}, r]
185 : ! Note that V_nl is present only if the AE-Hamiltonian is non-local e.g. DFT+U or LEXX.
186 32607 : spinorwf_pad=RESHAPE((/0,0,npw,npw,0,npw,npw,0/),(/2,4/))
187 18115 : ihr_comm=zero
188 :
189 : ! -i <c,k|\nabla_r|v,k> = \sum_G u_{ck}^*(G) [k+G] u_{vk}(G) in reduced coordinates.
190 3623 : if (istwfk==1) then
191 7246 : do iab=1,nspinor**2
192 3623 : spad1 = spinorwf_pad(1,iab)
193 3623 : spad2 = spinorwf_pad(2,iab)
194 1373867 : do ig=1,npw
195 1366621 : ctemp = CONJG(ug1(ig+spad1)) * ug2(ig+spad2)
196 5470107 : ihr_comm(:,iab) = ihr_comm(:,iab) + ctemp* ( kpoint + gvec(:,ig))
197 : end do
198 : end do
199 : else
200 : ! Symmetrized expression: \sum_G (k+G) 2i Im [ u_a^*(G) u_b(G) ]. (k0,G0) term is null.
201 0 : do ig=1,npw
202 0 : ctemp = CONJG(ug1(ig)) * ug2(ig)
203 0 : ihr_comm(:,1) = ihr_comm(:,1) + two*j_dpc * AIMAG(ctemp) * (kpoint + gvec(:,ig))
204 : end do
205 : end if
206 : !
207 : ! Add on-site terms.
208 3623 : ons_cart=zero
209 3623 : ABI_CHECK(nspinor==1,"nspinor/=1 not coded")
210 :
211 12885 : do iatom=1,Cryst%natom
212 9262 : itypat=Cryst%typat(iatom)
213 9262 : lmn_size=Pawtab(itypat)%lmn_size
214 9262 : nabla_ij => Pawtab(itypat)%nabla_ij(:,:,:)
215 : !
216 : !=== Unpacked loop over lmn channels ====
217 87332 : do jlmn=1,lmn_size
218 715600 : do ilmn=1,lmn_size
219 : re_p = Cprj_kb1(iatom,1)%cp(1,ilmn)*Cprj_kb2(iatom,1)%cp(1,jlmn) &
220 631891 : & +Cprj_kb1(iatom,1)%cp(2,ilmn)*Cprj_kb2(iatom,1)%cp(2,jlmn)
221 :
222 : im_p = Cprj_kb1(iatom,1)%cp(1,ilmn)*Cprj_kb2(iatom,1)%cp(2,jlmn) &
223 631891 : & -Cprj_kb1(iatom,1)%cp(2,ilmn)*Cprj_kb2(iatom,1)%cp(1,jlmn)
224 :
225 : ! Onsite contribution given by -i\nabla.
226 631891 : ons_cart(1,1)=ons_cart(1,1) + im_p*nabla_ij(1,ilmn,jlmn)
227 631891 : ons_cart(1,2)=ons_cart(1,2) + im_p*nabla_ij(2,ilmn,jlmn)
228 631891 : ons_cart(1,3)=ons_cart(1,3) + im_p*nabla_ij(3,ilmn,jlmn)
229 :
230 631891 : ons_cart(2,1)=ons_cart(2,1) - re_p*nabla_ij(1,ilmn,jlmn)
231 631891 : ons_cart(2,2)=ons_cart(2,2) - re_p*nabla_ij(2,ilmn,jlmn)
232 631891 : ons_cart(2,3)=ons_cart(2,3) - re_p*nabla_ij(3,ilmn,jlmn)
233 : !
234 706338 : if (Pawtab(itypat)%usepawu/=0) then ! Add i[V_u, r]
235 0 : isel=Hur(iatom)%ij_select(ilmn,jlmn,isppol)
236 0 : if (isel>0) then
237 0 : hurc_ij(:)=Hur(iatom)%commutator(:,isel,isppol)
238 :
239 0 : ons_cart(1,1)=ons_cart(1,1) - im_p*hurc_ij(1)
240 0 : ons_cart(1,2)=ons_cart(1,2) - im_p*hurc_ij(2)
241 0 : ons_cart(1,3)=ons_cart(1,3) - im_p*hurc_ij(3)
242 :
243 0 : ons_cart(2,1)=ons_cart(2,1) + re_p*hurc_ij(1)
244 0 : ons_cart(2,2)=ons_cart(2,2) + re_p*hurc_ij(2)
245 0 : ons_cart(2,3)=ons_cart(2,3) + re_p*hurc_ij(3)
246 : end if
247 : end if
248 :
249 : end do !ilmn
250 : end do !jlmn
251 : end do !iatom
252 :
253 : ! ons_cart is in Cartesian coordinates in real space
254 : ! while ihr_comm is in reduced coordinates in reciprocal space in terms of gprimd.
255 : !rs_cart2red = TRANSPOSE(Cryst%gprimd) ! if <r> is in terms of real space vectors
256 47099 : gspace_cart2red = TRANSPOSE(Cryst%rprimd)
257 :
258 : !ons_comm_red(1,:)=MATMUL(rs_cart2red,ons_comm(1,:))
259 : !ons_comm_red(2,:)=MATMUL(rs_cart2red,ons_comm(2,:))
260 : !ihr_comm(:,1) = ihr_comm(:,1) + CMPLX(ons_comm_red(1,:),ons_comm_red(2,:),kind=gwp)
261 :
262 123182 : ihr_comm_cart(:,1) = two_pi*MATMUL(Cryst%gprimd,ihr_comm(:,1))
263 14492 : ihr_comm_cart(:,1) = ihr_comm_cart(:,1) + CMPLX(ons_cart(1,:),ons_cart(2,:),kind=gwp)
264 :
265 : ! Final result is in reduced coordinates, in terms of gprimd.
266 123182 : ihr_comm(:,1) = MATMUL(gspace_cart2red, ihr_comm_cart(:,1))/two_pi
267 :
268 3623 : end function paw_ihr
269 : !!***
270 :
271 : !----------------------------------------------------------------------
272 :
273 : !!****f* m_paw_hr/paw_cross_ihr_comm
274 : !! NAME
275 : !! paw_cross_ihr_comm
276 : !!
277 : !! FUNCTION
278 : !! Adds the PAW cross term contribution to the matrix elements of the i\nabla operator.
279 : !! in cartesian coordinates. Should take also into account the contribution arising from the U
280 : !! part of the Hamiltonian (if any)
281 : !!
282 : !! INPUTS
283 : !! ihr_comm = the commutator [H,r] evaluated between states i and j, with only the plane-wave and
284 : !! the onsite parts included
285 : !! isppol=Spin index.
286 : !! nspinor=Number of spinori components.
287 : !! nr=Number real-space points on the fine fft grid for the ae wavefunctions
288 : !! kpoint(3)=k-point in reduced coordinates.
289 : !! Cryst<crystal_t>=Info on the crystal structure.
290 : !! %natom=Number of atoms in unit cell
291 : !! %typat(natom)
292 : !! Pawfgrtab(ntypat)= PAW tabulated data on the fine grid
293 : !! %lmn_size Number of (l,m,n) elements for the paw basis
294 : !! %nfgr Number of points on the fine grid
295 : !! %ifftsph Indexes of the fine-grid points on the fft mesh
296 : !! Paw_onsite(ntypat)= PAW tabulated data on the fine grid points inside the sphere
297 : !! %phi_gr(3,nfgr,lmn_size) gradient of phi in cartesian coordinates
298 : !! %tphi_gr(3,nfgr,lmn_size) gradient of tphi in cartesian coordinates
299 : !! ur_ae1(nr),ur_ae2(nr)=Left and right AE wavefunction.
300 : !! ur_ae_onsite1(nr),ur_ae_onsite2(nr)=Left and right AE onsite wavefunction.
301 : !! Cprj_kb1(natom,nspinor),Cprj_kb2(natom,nspinor) <type(pawcprj_type)>=
302 : !! projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to
303 : !! wavefunctions (k,b1,s) and (k,b2,s), respectively.
304 : !!
305 : !! OUTPUT
306 : !!
307 : !! SIDE EFFECTS
308 : !! The cross-term contribution is added to the commutator
309 : !!
310 : !! SOURCE
311 :
312 0 : subroutine paw_cross_ihr_comm(ihr_comm,nspinor,nr,Cryst,Pawfgrtab,Paw_onsite,&
313 0 : & ur_ae1,ur_ae2,ur_ae_onsite1,ur_ae_onsite2,Cprj_kb1,Cprj_kb2)
314 :
315 : !Arguments ------------------------------------
316 : !scalars
317 : integer,intent(in) :: nspinor,nr
318 : type(crystal_t),intent(in) :: Cryst
319 : !arrays
320 : complex(gwp),intent(inout) :: ihr_comm(3,nspinor**2)
321 : complex(gwp),intent(in) :: ur_ae1(nr),ur_ae2(nr)
322 : complex(gwp),intent(in) :: ur_ae_onsite1(nr),ur_ae_onsite2(nr)
323 : type(pawfgrtab_type),intent(in) :: Pawfgrtab(Cryst%natom)
324 : type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
325 : type(pawcprj_type),intent(in) :: Cprj_kb1(Cryst%natom,nspinor),Cprj_kb2(Cryst%natom,nspinor)
326 :
327 : !Local variables-------------------------------
328 : integer :: iatom,lmn_size,ilmn,ifgd,ifftsph,nfgd
329 : complex(dp) :: cp1, cp2
330 : complex(dp) :: cross1,cross2
331 : !arrays
332 : real(dp) :: gspace_cart2red(3,3)
333 0 : complex(gwp) :: ihr_comm_cart(3,nspinor**2)
334 : complex(dp) :: dphigr(3), dphigr1(3),dphigr2(3)
335 : ! *************************************************************************
336 :
337 0 : ABI_CHECK(nspinor==1,"nspinor + pawcross not implemented")
338 :
339 : ! [H, r] = -\nabla + [V_{nl}, r]
340 : ! The V_nl part, present in case of DFT+U, is omitted for the cross terms contribution
341 : ! Recall that delta_rho_tw_ij = (psi_i - phi_i)* (phi_j - tphi_j) + (phi_i - tphi_i)* (psi_j - phi_j)
342 0 : ihr_comm_cart(:,1) = czero
343 :
344 0 : do iatom=1,Cryst%natom
345 0 : lmn_size = Paw_onsite(iatom)%lmn_size
346 0 : nfgd = Pawfgrtab(iatom)%nfgd
347 :
348 0 : do ifgd=1,nfgd
349 :
350 0 : ifftsph = Pawfgrtab(iatom)%ifftsph(ifgd)
351 :
352 0 : cross1 = ur_ae1(ifftsph) - ur_ae_onsite1(ifftsph)
353 0 : cross2 = ur_ae2(ifftsph) - ur_ae_onsite2(ifftsph)
354 :
355 0 : do ilmn=1,lmn_size
356 :
357 0 : dphigr(1:3) = Paw_onsite(iatom)%phi_gr(1:3,ifgd,ilmn) - Paw_onsite(iatom)%tphi_gr(1:3,ifgd,ilmn)
358 :
359 0 : cp1 = CMPLX(Cprj_kb1(iatom,1)%cp(1,ilmn),Cprj_kb1(iatom,1)%cp(2,ilmn)) * sqrt(Cryst%ucvol) ! that damn magic factor
360 0 : cp2 = CMPLX(Cprj_kb2(iatom,1)%cp(1,ilmn),Cprj_kb2(iatom,1)%cp(2,ilmn)) * sqrt(Cryst%ucvol)
361 :
362 0 : dphigr1(1:3) = cp1 * dphigr(1:3)
363 0 : dphigr2(1:3) = cp2 * dphigr(1:3)
364 :
365 0 : ihr_comm_cart(1,1) = ihr_comm_cart(1,1) - j_dpc * (CONJG(cross1) * dphigr2(1) - CONJG(dphigr1(1)) * cross2) / nr
366 0 : ihr_comm_cart(2,1) = ihr_comm_cart(2,1) - j_dpc * (CONJG(cross1) * dphigr2(2) - CONJG(dphigr1(2)) * cross2) / nr
367 0 : ihr_comm_cart(3,1) = ihr_comm_cart(3,1) - j_dpc * (CONJG(cross1) * dphigr2(3) - CONJG(dphigr1(3)) * cross2) / nr
368 :
369 : end do
370 : end do
371 : end do
372 :
373 : ! Go to reduced coordinate
374 0 : gspace_cart2red = TRANSPOSE(Cryst%rprimd)
375 0 : ihr_comm(:,1) = ihr_comm(:,1) + MATMUL(gspace_cart2red, ihr_comm_cart(:,1)) / two_pi
376 :
377 0 : end subroutine paw_cross_ihr_comm
378 : !!***
379 :
380 : !----------------------------------------------------------------------
381 :
382 : !!****f* m_paw_hr/pawhur_init
383 : !! NAME
384 : !! pawhur_init
385 : !!
386 : !! FUNCTION
387 : !! Creation method for the pawhur_t data type.
388 : !!
389 : !! INPUTS
390 : !!
391 : !! OUTPUT
392 : !!
393 : !! SOURCE
394 :
395 0 : subroutine pawhur_init(hur,nsppol,pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij)
396 :
397 : !Arguments ------------------------------------
398 : !scalars
399 : integer,intent(in) :: nsppol,pawprtvol
400 : type(crystal_t),intent(in) :: Cryst
401 : type(Pawang_type),intent(in) :: Pawang
402 : !arrays
403 : type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
404 : type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)
405 : type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom)
406 : type(pawhur_t),intent(inout) :: Hur(Cryst%natom)
407 :
408 : !Local variables-------------------------------
409 : !scalars
410 : integer :: iatom,ij_idx,isel,itypat,isppol,lmn2_size_max,lmn2_size,lmn_size,lpawu
411 : integer :: jlmn,jl,jm,jlm,jln,k0lmn,k0lm,k0ln,ilmn,il,im,ilm,iln
412 : integer :: m2,m1,left_lmn,right_lmn,tot_lmn,nmax
413 : !arrays
414 0 : integer :: nsel(3,nsppol)
415 0 : integer, ABI_CONTIGUOUS pointer :: indlmn(:,:)
416 : real(dp) :: sumr_ij(3)
417 0 : real(dp),allocatable :: rcart_onsite(:,:,:)
418 0 : real(dp),allocatable :: rij_tmp(:,:,:),vpawu(:,:,:,:)
419 :
420 : ! *************************************************************************
421 :
422 : ! Get onsite matrix elements of the position operator.
423 0 : lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size)
424 0 : ABI_MALLOC(rcart_onsite,(3,lmn2_size_max,Cryst%natom))
425 :
426 0 : call pawr(Pawtab,Pawrad,Pawang,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,lmn2_size_max,rcart_onsite)
427 :
428 0 : do iatom=1,Cryst%natom
429 0 : itypat=Cryst%typat(iatom)
430 0 : if (Pawtab(itypat)%usepawu==0) CYCLE
431 0 : lmn2_size=Pawtab(itypat)%lmn2_size
432 0 : lmn_size =Pawtab(itypat)%lmn_size
433 0 : lpawu=Pawtab(itypat)%lpawu
434 0 : Hur(iatom)%lmn2_size=lmn2_size
435 0 : Hur(iatom)%lmn_size =lmn_size
436 0 : Hur(iatom)%nsppol =nsppol
437 0 : indlmn => Pawtab(itypat)%indlmn
438 :
439 0 : ABI_MALLOC(rij_tmp,(3,lmn_size**2,nsppol))
440 0 : rij_tmp=zero
441 :
442 : ! Get Vpawu^{\sigma}_{m1,m2}
443 0 : ABI_MALLOC(vpawu,(Paw_ij(iatom)%cplex_dij,2*lpawu+1,2*lpawu+1,Paw_ij(iatom)%ndij))
444 : call pawpupot(Paw_ij(iatom)%cplex_dij,Paw_ij(iatom)%ndij,&
445 : & Paw_ij(iatom)%noccmmp,Paw_ij(iatom)%nocctot,&
446 0 : & pawprtvol,Pawtab(itypat),vpawu)
447 :
448 0 : do isppol=1,nsppol ! spinor not implemented
449 :
450 : ! === Loop on (jl,jm,jn) channels ===
451 : ij_idx=0
452 0 : do jlmn=1,lmn_size
453 0 : jl =indlmn(1,jlmn)
454 0 : jm =indlmn(2,jlmn)
455 0 : jlm=indlmn(4,jlmn)
456 0 : jln=indlmn(5,jlmn)
457 :
458 0 : k0lmn=jlmn*(jlmn-1)/2
459 0 : k0lm =jlm *(jlm -1)/2
460 0 : k0ln =jln *(jln -1)/2
461 : !
462 : ! === Loop on (il,im,in) channels ===
463 : ! * Looping over all ij components. Elements are not symmetric.
464 0 : do ilmn=1,lmn_size
465 0 : il =indlmn(1,ilmn)
466 0 : im =indlmn(2,ilmn)
467 0 : ilm=indlmn(4,ilmn)
468 0 : iln=indlmn(5,ilmn)
469 :
470 0 : ij_idx=ij_idx+1
471 :
472 : ! === Selection rules ===
473 0 : if (il/=lpawu.and.jl/=lpawu) CYCLE
474 :
475 0 : sumr_ij(:)=zero
476 0 : do m2=1,2*lpawu+1
477 0 : do m1=1,2*lpawu+1
478 0 : if (m1==(im-lpawu-1).and.il==lpawu) then
479 0 : left_lmn =ilmn-(il+im+1)+m2
480 0 : right_lmn=jlmn
481 0 : if (right_lmn>=left_lmn) then
482 0 : tot_lmn=right_lmn*(right_lmn-1)/2 + left_lmn
483 : else
484 0 : tot_lmn=left_lmn*(left_lmn-1)/2 + right_lmn
485 : end if
486 0 : sumr_ij=sumr_ij+vpawu(1,m1,m2,isppol)*rcart_onsite(:,tot_lmn,iatom)
487 : end if
488 :
489 0 : if (m2==(jm-lpawu-1).and.jl==lpawu) then
490 0 : left_lmn =ilmn
491 0 : right_lmn=jlmn-(jl+jm+1)+m1
492 0 : if (right_lmn>=left_lmn) then
493 0 : tot_lmn=right_lmn*(right_lmn-1)/2 + left_lmn
494 : else
495 0 : tot_lmn=left_lmn*(left_lmn-1)/2 + right_lmn
496 : end if
497 0 : sumr_ij=sumr_ij+vpawu(1,m1,m2,isppol)*rcart_onsite(:,tot_lmn,iatom)
498 : end if
499 : end do !m1
500 : end do !m2
501 :
502 0 : rij_tmp(:,ij_idx,isppol)=sumr_ij(:)
503 :
504 : end do !ilmn
505 : end do !jlmn
506 : end do !isppol
507 :
508 0 : ABI_FREE(vpawu)
509 :
510 : ! === Save values in packed form ===
511 0 : ABI_MALLOC(Hur(iatom)%ij_select,(lmn_size,lmn_size,nsppol))
512 0 : Hur(iatom)%ij_select=0
513 0 : nsel(:,:)=COUNT(ABS(rij_tmp)>tol6,DIM=2)
514 0 : nmax=MAXVAL(nsel)
515 0 : ABI_MALLOC(Hur(iatom)%commutator,(3,nmax,nsppol))
516 0 : do isppol=1,nsppol
517 : ij_idx=0
518 : isel =0
519 0 : do jlmn=1,lmn_size
520 0 : do ilmn=1,lmn_size
521 0 : ij_idx=ij_idx+1
522 0 : if (ANY (ABS(rij_tmp(:,ij_idx,isppol))>tol6) ) then
523 0 : isel=isel+1
524 0 : Hur(iatom)%ij_select(ilmn,jlmn,isppol)=isel
525 0 : Hur(iatom)%commutator(:,isel,isppol)=rij_tmp(:,ij_idx,isppol)
526 : end if
527 : end do
528 : end do
529 : end do
530 :
531 0 : ABI_FREE(rij_tmp)
532 : end do !iatom
533 :
534 0 : ABI_FREE(rcart_onsite)
535 :
536 0 : end subroutine pawhur_init
537 : !!***
538 :
539 : !----------------------------------------------------------------------
540 :
541 : !!****f* m_paw_hr/pawr
542 : !! NAME
543 : !! pawr
544 : !!
545 : !! FUNCTION
546 : !! Evaluate matrix elements of the position operator between PAW AE partial waves.
547 : !!
548 : !! INPUTS
549 : !! Pawtab(ntypat) <type(pawtab_type)>=paw tabulated data read at start:
550 : !! %lmn_size
551 : !! %lmn2_size
552 : !! %indklmn
553 : !! %phiphj
554 : !! Pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
555 : !! %mesh_size=Dimension of radial mesh
556 : !! %rad(mesh_size)=The coordinates of all the points of the radial mesh
557 : !! Pawang <type(pawang_type)>=paw angular mesh and related data
558 : !! %lmax=Maximum value of angular momentum l+1
559 : !! %gntselect((2*l_max-1)**2,l_max**2,l_max**2)= selection rules for Gaunt coefficients
560 : !! %realgnt
561 : !! natom=number of atoms in unit cell
562 : !! ntypat=number of types of atom
563 : !! typat(natom)=type of each atom
564 : !! xcart(3,natom)=cartesian coordinates
565 : !!
566 : !! OUTPUT
567 : !! rcart_onsite(3,lmn2_size_max,natom)
568 : !!
569 : !! SOURCE
570 :
571 0 : subroutine pawr(Pawtab,Pawrad,Pawang,natom,ntypat,typat,xcart,lmn2_size_max,rcart_onsite)
572 :
573 : !Arguments ------------------------------------
574 : !scalars
575 : integer,intent(in) :: lmn2_size_max,natom,ntypat
576 : type(Pawang_type),intent(in) :: Pawang
577 :
578 : !arrays
579 : integer,intent(in) :: typat(natom)
580 : real(dp),intent(in) :: xcart(3,natom)
581 : real(dp),intent(inout) :: rcart_onsite(3,lmn2_size_max,natom)
582 : type(Pawrad_type),intent(in) :: Pawrad(ntypat)
583 : type(Pawtab_type),target,intent(in) :: Pawtab(ntypat)
584 :
585 : !Local variables-------------------------------
586 : !scalars
587 : integer,parameter :: ll1=1
588 : integer :: iatom,idir,ignt,il,ilm,ilm_G,ilmn,iln,im,itypat,jl,jlm,jlmn,jln,jm,k0lm
589 : integer :: k0lmn,k0ln,klm,klmn,kln,lmn_size,mesh_size,mm_G,lmn2_size
590 : real(dp) :: fact,intff,rgnt
591 : !arrays
592 0 : integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
593 0 : real(dp),allocatable :: ff(:),rad(:),rc_tmp(:,:)
594 :
595 : ! *************************************************************************
596 :
597 : DBG_ENTER("COLL")
598 :
599 0 : fact=two*SQRT(pi/three)
600 0 : rcart_onsite(:,:,:)=zero
601 :
602 0 : do itypat=1,ntypat
603 0 : lmn_size =Pawtab(itypat)%lmn_size
604 0 : lmn2_size =Pawtab(itypat)%lmn2_size
605 0 : mesh_size =Pawtab(itypat)%mesh_size
606 0 : indlmn => Pawtab(itypat)%indlmn
607 :
608 0 : ABI_MALLOC(ff,(mesh_size))
609 0 : ABI_MALLOC(rad,(mesh_size))
610 0 : rad(1:mesh_size)=Pawrad(itypat)%rad(1:mesh_size)
611 :
612 0 : ABI_MALLOC(rc_tmp,(3,lmn2_size))
613 0 : rc_tmp=zero
614 : !
615 : ! === Loop on (jl,jm,jn) channels
616 0 : do jlmn=1,lmn_size
617 0 : jl =indlmn(1,jlmn)
618 0 : jm =indlmn(2,jlmn)
619 0 : jlm=indlmn(4,jlmn)
620 0 : jln=indlmn(5,jlmn)
621 :
622 0 : k0lmn=jlmn*(jlmn-1)/2
623 0 : k0lm =jlm *(jlm -1)/2
624 0 : k0ln =jln *(jln -1)/2
625 : !
626 : ! === Loop on (il,im,in) channels; klmn is the index for packed form ===
627 0 : do ilmn=1,jlmn
628 0 : il =indlmn(1,ilmn)
629 0 : im =indlmn(2,ilmn)
630 0 : ilm=indlmn(4,ilmn)
631 0 : iln=indlmn(5,ilmn)
632 :
633 0 : klmn=k0lmn+ilmn
634 0 : klm =k0lm +ilm
635 0 : kln =k0ln +iln
636 : !
637 : ! === For each cartesian direction, use expansion in terms of RSH ===
638 : ! TODO Add a check if l=1 is in the set
639 0 : do idir=1,3
640 0 : mm_G=0
641 0 : if (idir==1) mm_G= 1
642 0 : if (idir==2) mm_G=-1
643 0 : if (idir==3) mm_G= 0
644 0 : ilm_G=1+ll1**2+ll1+mm_G
645 0 : ignt=Pawang%gntselect(ilm_G,klm)
646 0 : if (ignt/=0) then
647 0 : rgnt=Pawang%realgnt(ignt)
648 0 : ff(1)=zero
649 : !ff(2:mesh_size)=(Pawtab(itypat)%phiphj(2:mesh_size,kln)-Pawtab(itypat)%tphitphj(2:mesh_size,kln))*rad(2:mesh_size)
650 0 : ff(2:mesh_size)=Pawtab(itypat)%phiphj(2:mesh_size,kln)*rad(2:mesh_size)
651 0 : call simp_gen(intff,ff,Pawrad(itypat))
652 0 : rc_tmp(idir,klmn)=fact*intff*rgnt
653 : end if
654 : end do !idir
655 :
656 : end do !ilmn
657 : end do !jllmn
658 :
659 : ! === Make matrix elements for each atom of this type ===
660 0 : do jlmn=1,lmn_size
661 0 : jl =indlmn(1,jlmn)
662 0 : jm =indlmn(2,jlmn)
663 0 : jln=indlmn(5,jlmn)
664 :
665 0 : k0lmn=jlmn*(jlmn-1)/2
666 0 : k0ln =jln *(jln -1)/2
667 0 : do ilmn=1,jlmn
668 0 : il =indlmn(1,ilmn)
669 0 : im =indlmn(2,ilmn)
670 0 : iln=indlmn(5,ilmn)
671 :
672 0 : klmn=k0lmn+ilmn
673 0 : kln =k0ln +iln
674 :
675 0 : intff=zero
676 0 : if (il==jl.and.jm==im) then
677 0 : ff(1:mesh_size)=Pawtab(itypat)%phiphj(1:mesh_size,kln)
678 0 : call simp_gen(intff,ff,Pawrad(itypat))
679 : end if
680 0 : do iatom=1,natom
681 0 : if (typat(iatom)/=itypat) CYCLE
682 0 : rcart_onsite(:,klmn,iatom)=rc_tmp(:,klmn) + xcart(:,iatom)*intff
683 : end do
684 :
685 : end do ! ilmn
686 : end do !jlmn
687 :
688 0 : ABI_FREE(ff)
689 0 : ABI_FREE(rad)
690 0 : ABI_FREE(rc_tmp)
691 : end do !itypat
692 :
693 : DBG_EXIT("COLL")
694 :
695 0 : end subroutine pawr
696 : !!***
697 :
698 : !----------------------------------------------------------------------
699 :
700 0 : END MODULE m_paw_hr
701 : !!***
|
