Branch data 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 : : !!***
|
