Line data Source code
1 : !!****m* ABINIT/m_vkbr
2 : !! NAME
3 : !! m_vkbr
4 : !!
5 : !! FUNCTION
6 : !! This module provides objects and methods used to calculate the matrix elements
7 : !! of the commutator [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> when non-local pseudopotentials are used.
9 : !!
10 : !! NOTES
11 : !! This module is deprecated. Use ddkop_t in m_ddk.F90
12 : !!
13 : !! COPYRIGHT
14 : !! Copyright (C) 2008-2026 ABINIT group (MG, FB)
15 : !! This file is distributed under the terms of the
16 : !! GNU General Public License, see ~abinit/COPYING
17 : !! or http://www.gnu.org/copyleft/gpl.txt .
18 : !!
19 : !! SOURCE
20 :
21 : #if defined HAVE_CONFIG_H
22 : #include "config.h"
23 : #endif
24 :
25 : #include "abi_common.h"
26 :
27 : MODULE m_vkbr
28 :
29 : use defs_basis
30 : use m_hide_blas
31 : use m_errors
32 : use m_abicore
33 :
34 : use defs_datatypes, only : pseudopotential_type
35 : use m_gwdefs, only : czero_gw
36 : use m_fstrings, only : sjoin, itoa
37 : use m_paw_sphharm, only : ylmc, ylmcd
38 : use m_geometry, only : normv
39 : use m_crystal, only : crystal_t
40 : use m_kg, only : mkkin
41 : use m_mkffnl, only : mkffnl
42 :
43 : implicit none
44 :
45 : private
46 : !!***
47 :
48 : !----------------------------------------------------------------------
49 :
50 : !!****t* m_vkbr/vkbr_t
51 : !! NAME
52 : !!
53 : !! FUNCTION
54 : !! Matrix elements in |k+G> space needed for the
55 : !! evaluation of the matrix elements of the commutator [Vnl,r] for the
56 : !! optical limit in <kb1|e^{-iqr}|kb2>.
57 : !!
58 : !! SOURCE
59 :
60 : type,public :: vkbr_t
61 :
62 : integer :: istwfk
63 : ! Storage mode of the G vectors for this k-point.
64 :
65 : integer :: ntypat
66 : ! Number of type of atoms
67 :
68 : integer :: natom
69 : ! Number of atoms
70 :
71 : integer :: mpsang
72 : ! Max l+1 over atoms
73 :
74 : integer :: npw
75 : ! Number of G-vectors.
76 :
77 : integer :: inclvkb
78 : ! Option for calculating the matrix elements of [Vnl,r].
79 : ! 0 to exclude commutator, 2 to include it
80 :
81 : real(dp) :: kpoint(3)
82 : ! The k-point in reduced coordinates.
83 :
84 : complex(gwp),allocatable :: fnl(:,:,:,:)
85 : ! fnl(npw,mpsang**2,mproj,natom)
86 :
87 : complex(gwp),allocatable :: fnld(:,:,:,:,:)
88 : ! fnld(3,npw,mpsang**2,mproj,natom)
89 :
90 : end type vkbr_t
91 :
92 : public :: vkbr_init ! vkbr_t Constructor
93 : public :: vkbr_free ! Free memory
94 : public :: nc_ihr_comm ! Compute matrix elements of the commutator i[H,r] for NC pseudos
95 : public :: calc_vkb ! Kleynman-Bylander form factors and derivatives.
96 : !!***
97 :
98 : interface vkbr_free
99 : module procedure vkbr_free_0D
100 : module procedure vkbr_free_1D
101 : end interface vkbr_free
102 :
103 : CONTAINS !========================================================================================
104 :
105 : !----------------------------------------------------------------------
106 :
107 : !!****f* m_vkbr/vkbr_init
108 : !! NAME
109 : !! vkbr_init
110 : !!
111 : !! FUNCTION
112 : !! Creation method the the vkbr_t structures datatype.
113 : !!
114 : !! INPUTS
115 : !! cryst<crystal_t>=Datatype gathering info on the crystal structure.
116 : !! psps<pseudopotential_type>Structure gathering info on the pseudopotentials.
117 : !! inclvkb=Option defining the algorithm used for the application of [Vnl,r].
118 : !! 2 for Spherical harmonics
119 : !! istwfk=Storage mode for the wavefunctions at this k-point.
120 : !! npw=Number of planewaves in <k+G1|[Vnl,r]|k+G2>
121 : !! kpoint(3)=K-point of interest in reduced coordinates.
122 : !! gvec(3,npw)=Reduced coordinates of the G-vectors.
123 : !!
124 : !! OUTPUT
125 : !! vkbr<vkbr_t>=Structure containing arrays needed for calculating <\psi_1|[Vnl,r]\psi_2>.
126 : !! Completely initialized in output.
127 : !!
128 : !! SOURCE
129 :
130 7208 : subroutine vkbr_init(vkbr,cryst,psps,inclvkb,istwfk,npw,kpoint,gvec)
131 :
132 : !Arguments ------------------------------------
133 : !scalars
134 : integer,intent(in) :: npw,inclvkb,istwfk
135 : type(crystal_t),intent(in) :: cryst
136 : type(vkbr_t),intent(inout) :: vkbr
137 : type(pseudopotential_type),intent(in) :: psps
138 : !arrays
139 : integer,intent(in) :: gvec(3,npw)
140 : real(dp),intent(in) :: kpoint(3)
141 :
142 : !Local variables-------------------------------
143 : !scalars
144 : integer :: ierr
145 : character(len=500) :: msg
146 : !arrays
147 7208 : real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:)
148 : !************************************************************************
149 :
150 : !@vkbr_t
151 7208 : vkbr%istwfk = istwfk
152 7208 : vkbr%ntypat = cryst%ntypat
153 7208 : vkbr%natom = cryst%natom
154 7208 : vkbr%mpsang = psps%mpsang
155 7208 : vkbr%npw = npw
156 7208 : vkbr%inclvkb = inclvkb
157 28832 : vkbr%kpoint = kpoint
158 :
159 : ! Calculate KB form factors and derivatives.
160 : ! The arrays are allocated with lnmax to support pseudos with more than projector.
161 : ! Note that lnmax takes into account lloc hence arrays are in packed form and one should be
162 : ! accessed with the indices provided by indlmn.
163 : ! TODO: they should be calculated on-the-fly using calc_vkb
164 : ! For the moment, we opt for a quick an dirty implementation.
165 :
166 28832 : ABI_MALLOC(vkbsign, (psps%lnmax, cryst%ntypat))
167 36040 : ABI_MALLOC(vkb, (npw, psps%lnmax, cryst%ntypat))
168 28832 : ABI_MALLOC(vkbd, (npw, psps%lnmax, cryst%ntypat))
169 7208 : call calc_vkb(cryst,psps,kpoint,npw,npw,gvec,vkbsign,vkb,vkbd)
170 :
171 7208 : select case (inclvkb)
172 : case (2)
173 : ! Complex spherical harmonics (CPU and mem \propto npw).
174 7208 : write(msg,'(a,f12.1)')'out-of-memory in fnl; Mb= ',one*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwp*b2Mb
175 43248 : ABI_STAT_MALLOC(vkbr%fnl,(npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr)
176 7208 : ABI_CHECK(ierr==0, msg)
177 :
178 7208 : write(msg,'(a,f12.1)')'out-of-memory in fnld; Mb= ',three*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwp*b2Mb
179 43248 : ABI_STAT_MALLOC(vkbr%fnld,(3,npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr)
180 7208 : ABI_CHECK(ierr==0, msg)
181 :
182 7208 : call ccgradvnl_ylm(cryst,psps,npw,gvec,kpoint,vkbsign,vkb,vkbd,vkbr%fnl,vkbr%fnld)
183 :
184 : case default
185 7208 : ABI_ERROR(sjoin("Wrong inclvkb= ",itoa(inclvkb)))
186 : end select
187 :
188 7208 : ABI_FREE(vkbsign)
189 7208 : ABI_FREE(vkb)
190 7208 : ABI_FREE(vkbd)
191 :
192 7208 : end subroutine vkbr_init
193 : !!***
194 :
195 : !----------------------------------------------------------------------
196 :
197 : !!****f* m_vkbr/vkbr_free_0D
198 : !! NAME
199 : !! vkbr_free_0D
200 : !!
201 : !! FUNCTION
202 : !! Free all memory allocated in a structure of type vkbr_t
203 : !!
204 : !! SOURCE
205 :
206 8085 : subroutine vkbr_free_0D(vkbr)
207 :
208 : !Arguments ------------------------------------
209 : type(vkbr_t),intent(inout) :: vkbr
210 : !************************************************************************
211 :
212 : !complex
213 8085 : ABI_SFREE(vkbr%fnl)
214 8085 : ABI_SFREE(vkbr%fnld)
215 :
216 8085 : end subroutine vkbr_free_0D
217 : !!***
218 :
219 : !----------------------------------------------------------------------
220 :
221 : !!****f* m_vkbr/vkbr_free_1D
222 : !! NAME
223 : !! vkbr_free_1D
224 : !!
225 : !! FUNCTION
226 : !! Free all memory allocated in a structure of type vkbr_t
227 : !!
228 : !! SOURCE
229 :
230 146 : subroutine vkbr_free_1D(vkbr)
231 :
232 : !Arguments ------------------------------------
233 : !arrays
234 : type(vkbr_t),intent(inout) :: vkbr(:)
235 :
236 : !Local variables ------------------------------
237 : !scalars
238 : integer :: ii
239 : !************************************************************************
240 :
241 995 : do ii=1,SIZE(vkbr)
242 995 : call vkbr_free_0D(vkbr(ii))
243 : end do
244 :
245 146 : end subroutine vkbr_free_1D
246 : !!***
247 :
248 : !----------------------------------------------------------------------
249 :
250 : !!****f* m_vkbr/add_vnlr_commutator
251 : !! NAME
252 : !! add_vnlr_commutator
253 : !!
254 : !! FUNCTION
255 : !! Calculate the matrix elements of the dipole operator <phi1|r|phi2>.
256 : !! For norm conserving potentials the commutator [Vnl,r] is included according to inclvkb.
257 : !!
258 : !! INPUTS
259 : !! vkbr<vkbr_t>
260 : !! cryst<crystal_t>=Datatype gathering info on the crystal structure.
261 : !! psps<pseudopotential_type>Structure gathering info on the pseudopotentials.
262 : !! npw=Number of G for wavefunctions.
263 : !! nspinor=Number of spinorial components.
264 : !! ug1(npw*nspinor)=Left wavefunction.
265 : !! ug2(npw*nspinor)=Right wavefunction
266 : !!
267 : !! SIDE EFFECTS
268 : !! rhotwx(3,nspinor**2)= Updated. Matrix elements in reduced coordinates, see NOTES below.
269 : !!
270 : !! NOTES
271 : !! 1) <k b1|e^{-iq.r}|k b2> = \delta_{b1 b2} -iq <k b1|r|k b2> = \delta_{b1 b2} -iq ( <k b1| [H,r] |k b2> / (e1-e2) ).
272 : !!
273 : !! This routine calculates the matrix elements of ir*(e1-e2)
274 : !! Remember that [H,r] = -\nabla + [V_nl,r]
275 : !!
276 : !! 2) The Fourier transform of a two-point real function f(r1,r2) satisfies:
277 : !! a) f_{\Gamma}(G1,G2) = f_{\Gamma}(-G1,-G2)^*
278 : !! b) f_{G0/2} (G1,G2) = f_{G0/2}(-G1-G0,-G2-G0)^*
279 : !!
280 : !! TODO
281 : !! *) Spinorial case is not implemented.
282 : !!
283 : !! SOURCE
284 :
285 262099 : subroutine add_vnlr_commutator(vkbr,cryst,psps,npw,nspinor,ug1,ug2,rhotwx)
286 :
287 : !Arguments ------------------------------------
288 : !scalars
289 : integer,intent(in) :: npw,nspinor
290 : type(vkbr_t),intent(in) :: vkbr
291 : type(crystal_t),intent(in) :: cryst
292 : type(pseudopotential_type),intent(in) :: psps
293 : !arrays
294 : complex(gwp),target,intent(in) :: ug1(npw*nspinor),ug2(npw*nspinor)
295 : complex(gwp),intent(inout) :: rhotwx(3,nspinor**2)
296 :
297 : !Local variables ------------------------------
298 : !scalars
299 : integer :: iat,ig,ilm,itypat,nlmn,ilmn,iln0,iln,il,in,im
300 : complex(gwp) :: cta1,cta4
301 : !arrays
302 : complex(gwp) :: dum(3),cta2(3),cta3(3),gamma_term(3)
303 : !************************************************************************
304 :
305 262099 : ABI_CHECK(nspinor == 1, "inclvkb > 0 with nspinor == 2 is not coded")
306 :
307 : ! Adding term i <c,k|[Vnl,r]|v,k> ===
308 524198 : select case (vkbr%inclvkb)
309 : case (2)
310 : ! Complex spherical harmonics (much faster!).
311 262099 : dum=czero_gw; gamma_term=czero
312 :
313 796803 : do iat=1,vkbr%natom
314 534704 : itypat = cryst%typat(iat)
315 2402146 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
316 : iln0 = 0
317 2639399 : do ilmn=1,nlmn
318 1842596 : il = 1 + psps%indlmn(1,ilmn,itypat)
319 1842596 : in = psps%indlmn(3,ilmn,itypat)
320 1842596 : iln = psps%indlmn(5,ilmn,itypat)
321 1842596 : if (iln <= iln0) cycle
322 1842596 : iln0 = iln
323 : !if (indlmn(6,ilmn,itypat) /= 1 .or. vkbsign(iln,itypat) == zero) cycle
324 : !in = 1
325 7183444 : do im=1,2*(il-1)+1
326 : ! Index of im and il
327 4806144 : ilm = im + (il-1)*(il-1)
328 4806144 : cta1 = czero_gw; cta2(:) = czero_gw
329 4806144 : cta4 = czero_gw; cta3(:) = czero_gw
330 2081840840 : do ig=1,npw
331 : ! Here we take advantage of the property Y_(l-m)= (-i)^m Y_lm^*.
332 2077034696 : cta1 = cta1 + ug1(ig) * vkbr%fnl (ig,ilm,in,iat)
333 8308138784 : cta2(:)= cta2(:) + ug2(ig) * vkbr%fnld(:,ig,ilm,in,iat)
334 8308138784 : cta3(:)= cta3(:) + ug1(ig) * vkbr%fnld(:,ig,ilm,in,iat)
335 2077034696 : cta4 = cta4 + ug2(ig) * vkbr%fnl (ig,ilm,in,iat)
336 2096259272 : if (ig==1) gamma_term = gamma_term + CONJG(cta1)*cta2(:) +CONJG(cta3(:))*cta4
337 : end do
338 21067172 : dum(:)= dum(:) + CONJG(cta1)*cta2(:) + CONJG(cta3(:))*cta4
339 : end do
340 :
341 : end do
342 : end do
343 :
344 262099 : if (vkbr%istwfk>1) then
345 16896 : dum = two * j_dpc * AIMAG(dum); if (vkbr%istwfk==2) dum = dum - j_dpc * AIMAG(gamma_term)
346 : end if
347 1048396 : rhotwx(:,1) = rhotwx(:,1) + dum(:)
348 :
349 : case default
350 262099 : ABI_ERROR(sjoin("Wrong inclvkb:", itoa(vkbr%inclvkb)))
351 : end select
352 :
353 262099 : end subroutine add_vnlr_commutator
354 : !!***
355 :
356 : !----------------------------------------------------------------------
357 :
358 : !!****f* m_vkbr/calc_vkb
359 : !! NAME
360 : !! calc_vkb
361 : !!
362 : !! FUNCTION
363 : !! This routine calculates the Kleynman-Bylander form factors and its derivatives
364 : !! needed for the evaluation of the matrix elements of the dipole operator <phi1|r|phi2>.
365 : !!
366 : !! INPUTS
367 : !! cryst<crystal_t>=Crystalline structure
368 : !! psps<pseudopotential_type>=Structured datatype gathering information on the pseudopotentials.
369 : !! kpoint(3)=The k-point in reduced coordinates.
370 : !! npw_k=Number of plane waves for this k-point.
371 : !! kg_k(3,npw_k)=Reduced coordinates of the G-vectors.
372 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
373 : !!
374 : !! OUTPUT
375 : !! vkb (npw_k, %lnmax, %ntypat)=KB form factors.
376 : !! vkbd(npw_k, %lnmax, %ntypat)=KB form factor derivatives.
377 : !! vkbsign(%lnmax, %ntypat) =KS dyadic sign.
378 : !!
379 : !! TODO
380 : !! SOC not implemented.
381 : !!
382 : !! SOURCE
383 :
384 7474 : subroutine calc_vkb(cryst,psps,kpoint,npw_k,mpw,kg_k,vkbsign,vkb,vkbd)
385 :
386 : !Arguments ------------------------------------
387 : !scalars
388 : integer,intent(in) :: npw_k, mpw
389 : type(crystal_t),intent(in) :: cryst
390 : type(pseudopotential_type),intent(in) :: psps
391 : !arrays
392 : integer,intent(in) :: kg_k(3,npw_k)
393 : real(dp),intent(in) :: kpoint(3)
394 : real(dp),intent(out) :: vkb (mpw,psps%lnmax,psps%ntypat)
395 : real(dp),intent(out) :: vkbd(mpw,psps%lnmax,psps%ntypat)
396 : real(dp),intent(out) :: vkbsign(psps%lnmax,psps%ntypat)
397 :
398 : !Local variables ------------------------------
399 : !scalars
400 : integer :: dimffnl,ider,idir,itypat,nkpg,in,il,ilmn,ig,iln,iln0,nlmn
401 : real(dp) :: effmass_free,ecutsm,ecut
402 : !arrays
403 7474 : real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:),ylm_gr(:,:,:),ylm_k(:,:)
404 : ! *************************************************************************
405 :
406 : DBG_ENTER("COLL")
407 7474 : ABI_CHECK(psps%usepaw==0, "You should not be here!")
408 7474 : ABI_CHECK(psps%useylm==0, "useylm/=0 not considered!")
409 :
410 : ! Compute KB dyadic sign.
411 43567 : vkbsign=zero
412 15026 : do itypat=1,psps%ntypat
413 7552 : iln0 = 0
414 36093 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
415 43474 : do ilmn=1,nlmn
416 28448 : iln = psps%indlmn(5,ilmn,itypat)
417 28448 : if (iln <= iln0) cycle
418 28448 : iln0 = iln
419 36000 : if (abs(psps%ekb(iln,itypat)) > 1.0d-10) vkbsign(iln,itypat) = dsign(one, psps%ekb(iln,itypat))
420 : end do
421 : end do
422 :
423 : ! Allocate KB form factor and derivative wrt k+G
424 : ! Here we do not use correct ordering for dimensions
425 7474 : idir=0; nkpg=0; ider=1; dimffnl=2 ! To retrieve the first derivative.
426 :
427 : ! Quantities used only if useylm==1
428 29896 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2*Psps%useylm))
429 37370 : ABI_MALLOC(ylm_gr, (npw_k, 3+6*(ider/2),psps%mpsang**2*Psps%useylm))
430 14948 : ABI_MALLOC(kpg_dum, (npw_k, nkpg))
431 44844 : ABI_MALLOC(ffnl, (npw_k,dimffnl, psps%lmnmax, psps%ntypat))
432 :
433 : call mkffnl(psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,cryst%gmet,cryst%gprimd,ider,idir,Psps%indlmn,&
434 : kg_k,kpg_dum,kpoint,psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,&
435 7474 : psps%ntypat,Psps%pspso,Psps%qgrid_ff,cryst%rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr)
436 :
437 7474 : ABI_FREE(ylm_k)
438 7474 : ABI_FREE(ylm_gr)
439 7474 : ABI_FREE(kpg_dum)
440 :
441 22422 : ABI_MALLOC(modkplusg, (npw_k))
442 7474 : effmass_free = one; ecutsm = zero; ecut = huge(one)
443 7474 : call mkkin(ecut,ecutsm,effmass_free,cryst%gmet,kg_k,modkplusg,kpoint,npw_k,0,0)
444 2389982 : modkplusg(:) = SQRT(half/pi**2*modkplusg(:))
445 2389982 : modkplusg(:) = MAX(modkplusg(:),tol10)
446 :
447 : ! Calculate matrix elements.
448 21326540 : vkb=zero; vkbd=zero
449 :
450 15026 : do itypat=1,psps%ntypat
451 7552 : iln0 = 0
452 36093 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
453 43474 : do ilmn=1,nlmn
454 28448 : il = 1 + psps%indlmn(1,ilmn,itypat)
455 28448 : in = psps%indlmn(3,ilmn,itypat)
456 28448 : iln = psps%indlmn(5,ilmn,itypat)
457 : !write(*,*)ilmn, iln, il, in
458 28448 : if (iln <= iln0) cycle
459 28448 : iln0 = iln
460 : !if (vkbsign(iln,itypat) == zero) cycle
461 36000 : if (ABS(psps%ekb(iln,itypat)) > 1.0d-10) then
462 28448 : ABI_CHECK(iln == ilmn, "iln != ilmn")
463 : !ABI_CHECK(il == iln, "il != iln")
464 28448 : if (il==1) then
465 3568535 : vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)
466 3568535 : vkbd(1:npw_k,iln,itypat) = ffnl(:,2,iln,itypat)*modkplusg(:)/two_pi
467 18589 : else if (il==2) then
468 3945843 : vkb(1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)*modkplusg(:)
469 3945843 : do ig=1,npw_k
470 : vkbd(ig,iln,itypat) = ((ffnl(ig,2,iln,itypat)*modkplusg(ig)*modkplusg(ig))+&
471 3945843 : ffnl(ig,1,iln,itypat) )/two_pi
472 : end do
473 7328 : else if (il==3) then
474 2606737 : vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)*modkplusg(:)**2
475 : vkbd(1:npw_k,iln,itypat) = (ffnl(:,2,iln,itypat)*modkplusg(:)**3+&
476 2606737 : 2*ffnl(:,1,iln,itypat)*modkplusg(:) )/two_pi
477 1536 : else if (il==4) then
478 435120 : vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)*modkplusg(:)**3
479 : vkbd(1:npw_k,iln,itypat) = (ffnl(:,2,iln,itypat)*modkplusg(:)**4+&
480 435120 : 3*ffnl(:,1,iln,itypat)*modkplusg(:)**2 )/two_pi
481 : end if
482 10619543 : vkb (:,iln,itypat) = SQRT(4*pi/cryst%ucvol*(2*il-1)*ABS(psps%ekb(iln,itypat)))*vkb (:,iln,itypat)
483 10619543 : vkbd(:,iln,itypat) = SQRT(4*pi/cryst%ucvol*(2*il-1)*ABS(psps%ekb(iln,itypat)))*vkbd(:,iln,itypat)
484 : end if
485 : end do
486 : end do
487 :
488 7474 : ABI_FREE(ffnl)
489 7474 : ABI_FREE(modkplusg)
490 :
491 : DBG_EXIT("COLL")
492 :
493 7474 : end subroutine calc_vkb
494 : !!***
495 :
496 : !----------------------------------------------------------------------
497 :
498 : !!****f* m_vkbr/nc_ihr_comm
499 : !! NAME
500 : !! nc_ihr_comm
501 : !!
502 : !! FUNCTION
503 : !! Calculate the matrix elements of the commutator i[H,r]
504 : !! For NC pseudppotentials, the commutator i[Vnl,r] is included depending on inclvkb.
505 : !!
506 : !! INPUTS
507 : !! vkbr<vkbr_t>
508 : !! cryst<crystal_t>=Unit cell and symmetries
509 : !! psps<pseudopotential_type>Structure gathering info on the pseudopotentials.
510 : !! nspinor=Number of spinorial components.
511 : !! npw=Number of G for wavefunctions.
512 : !! istwfk=Storage mode for wavefunctions.
513 : !! inclvkb=Option defining whether [Vnl,r] is added or not.
514 : !! kpoint(3)=k-point in reduced coordinates.
515 : !! ug1(npw*nspinor)=Left wavefunction.
516 : !! ug2(npw*nspinor)=Right wavefunction
517 : !! gvec(3,npw)=Planes waves for wavefunctions.
518 : !!
519 : !! OUTPUT
520 : !! ihr_comm(3,nspinor**2)= Matrix elements of the commutator i[H,r] between the input states.
521 : !! Result is in reduced coordinates. ug1 and ug2 are supposed to be orthogonal.
522 : !!
523 : !! NOTES
524 : !! <k b1|e^{-iq.r}|k b2> = \delta_{b1 b2} -iq <k b1|r|k b2> = \delta_{b1 b2} -iq ( <k b1| [H,r] |k b2> / (e1-e2) ).
525 : !! Remember that [H,r] = -\nabla + [V_nl,r]
526 : !!
527 : !! TODO
528 : !! *) Spinorial case is not implemented.
529 : !!
530 : !! SOURCE
531 :
532 306205 : function nc_ihr_comm(vkbr, cryst, psps, npw, nspinor, istwfk, inclvkb, kpoint, ug1, ug2, gvec) result(ihr_comm)
533 :
534 : !Arguments ------------------------------------
535 : !scalars
536 : integer,intent(in) :: npw,nspinor,inclvkb,istwfk
537 : type(vkbr_t),intent(in) :: vkbr
538 : type(crystal_t),intent(in) :: cryst
539 : type(Pseudopotential_type),intent(in) :: psps
540 : !arrays
541 : integer,intent(in) :: gvec(3,npw)
542 : real(dp),intent(in) :: kpoint(3)
543 : complex(gwp),intent(in) :: ug1(npw*nspinor),ug2(npw*nspinor)
544 : complex(gwp) :: ihr_comm(3,nspinor**2)
545 :
546 : !Local variables ------------------------------
547 : !scalars
548 : integer :: ig,iab,spad1,spad2
549 : complex(dp) :: c_tmp
550 : !arrays
551 : integer :: spinorwf_pad(2,4)
552 : !************************************************************************
553 :
554 : ! [H, r] = -\nabla + [V_{nl}, r]
555 : ! V_nl is present only in the case of NC pseudos but
556 : ! not in PAW unless even the AE Hamiltonian in non-local e.g. DFT+U or LEXX.
557 :
558 : ! -i <c,k|\nabla_r|v,k> in reduced coordinates is always included.
559 : ! -i <c,k|\nabla_r|v,k> = \sum_G u_{ck}^*(G) [k+G] u_{vk}(G)
560 : ! Note that here we assume c/=v, moreover the ug are supposed to be orthonormal and
561 : ! hence k+G can be replaced by G.
562 : ! HM 03/08/2018: we need band velocities so we don't assume c/=v anymore and we use k+G.
563 :
564 2755845 : spinorwf_pad = reshape([0, 0, npw, npw, 0, npw, npw, 0], [2, 4])
565 1551761 : ihr_comm = czero
566 :
567 : ! -i <c,k|\nabla_r|v,k> in reduced coordinates.
568 : ! This term is spin diagonal if nspinor == 2
569 306205 : if (istwfk == 1) then
570 598650 : do iab=1,nspinor
571 300189 : spad1 = spinorwf_pad(1,iab); spad2 = spinorwf_pad(2,iab)
572 95543042 : do ig=1,npw
573 94944392 : c_tmp = GWPC_CONJG(ug1(ig+spad1)) * ug2(ig+spad2)
574 380077757 : ihr_comm(:,iab) = ihr_comm(:,iab) + c_tmp * (kpoint + gvec(:,ig))
575 : end do
576 : end do
577 : else
578 : ! Symmetrized expression: \sum_G (k+G) 2i Im [ u_a^*(G) u_b(G) ]. (k0, G0) term is null.
579 7744 : ABI_CHECK(nspinor == 1, "nspinor != 1")
580 4078928 : do ig=1,npw
581 4071184 : c_tmp = GWPC_CONJG(ug1(ig)) * ug2(ig)
582 16292480 : ihr_comm(:,1) = ihr_comm(:,1) + two*j_dpc * AIMAG(c_tmp) * (kpoint + gvec(:,ig))
583 : end do
584 : end if
585 :
586 : ! Add second term $i <c,k|[Vnl,r]|v,k> $ in reduced cordinates.
587 306205 : if (inclvkb /= 0) then
588 262099 : ABI_CHECK(istwfk == vkbr%istwfk, "input istwfk /= vkbr%istwfk")
589 262099 : call add_vnlr_commutator(vkbr,cryst,psps,npw,nspinor,ug1,ug2,ihr_comm)
590 : end if
591 :
592 306205 : end function nc_ihr_comm
593 : !!***
594 :
595 : !----------------------------------------------------------------------
596 :
597 : !!****f* m_vkbr/ccgradvnl_ylm
598 : !! NAME
599 : !! ccgradvnl_ylm
600 : !!
601 : !! FUNCTION
602 : !! Compute Vnl(K) and grad_K Vnl(K) three reciprocal lattice units components
603 : !! using spherical harmonics instead of Legendre polynomials
604 : !! Needed for chi0(q=0)
605 : !!
606 : !! INPUTS
607 : !! cryst<crystal_t>=Unit cell and symmetries
608 : !! psps<pseudopotential_type>Structure gathering info on the pseudopotentials.
609 : !! npw=number of planewaves for wavefunctions
610 : !! gvec(3,npw)=integer coordinates of each plane wave in reciprocal space
611 : !! kpoint(3)=K-point in reduced coordinates.
612 : !! vkbsign(lnmax,ntypat)=sign of each KB dyadic product
613 : !! vkb(npw,lnmax,ntypat)=KB projector function
614 : !! vkbd(npw,lnmax,ntypat)=derivative of the KB projector function in reciprocal space
615 : !!
616 : !! OUTPUT
617 : !! fnl(npw,mpsang*2,natom),
618 : !! fnld(3,npw,mpsang*2,natom)
619 : !!
620 : !! NOTES
621 : !! Subroutine taken from the EXC code
622 : !! All the calculations are done in double precision, but the output arrays fnl and fnld
623 : !! are in single precision, should use double precision after modification of the other subroutines
624 : !!
625 : !! SOURCE
626 :
627 7208 : subroutine ccgradvnl_ylm(cryst,psps,npw,gvec,kpoint,vkbsign,vkb,vkbd,fnl,fnld)
628 :
629 : !Arguments ------------------------------------
630 : !scalars
631 : integer,intent(in) :: npw
632 : type(crystal_t),intent(in) :: cryst
633 : type(pseudopotential_type),intent(in) :: psps
634 : !arrays
635 : integer,intent(in) :: gvec(3,npw)
636 : real(dp),intent(in) :: kpoint(3)
637 : real(dp),intent(in) :: vkb(npw,psps%lnmax,cryst%ntypat)
638 : real(dp),intent(in) :: vkbd(npw,psps%lnmax,cryst%ntypat)
639 : real(dp),intent(in) :: vkbsign(psps%lnmax,cryst%ntypat)
640 : complex(gwp),intent(out) :: fnl(npw,psps%mpsang**2,psps%mproj,cryst%natom)
641 : complex(gwp),intent(out) :: fnld(3,npw,psps%mpsang**2,psps%mproj,cryst%natom)
642 :
643 : !Local variables-------------------------------
644 : !scalars
645 : integer :: ii,iat,ig,il,im,ilm,itypat,nlmn,iln0,iln,ilmn,in
646 : real(dp),parameter :: ppad=tol6
647 : real(dp) :: cosphi,costh,factor,mkg,mkg2,sinphi,sinth,sq,xdotg
648 : complex(dp) :: dphi,dth,sfac
649 : character(len=500) :: msg
650 : !arrays
651 : real(dp) :: gcart(3),kcart(3),kg(3)
652 : real(dp) :: b1(3),b2(3),b3(3),a1(3),a2(3),a3(3)
653 : complex(dp) :: dylmcart(3),dylmcrys(3),gradphi(3),gradth(3)
654 : !************************************************************************
655 :
656 : DBG_ENTER("COLL")
657 :
658 7208 : if (psps%mpsang > 4) then
659 : write(msg,'(3a)')&
660 0 : 'Number of angular momentum components bigger than programmed.',ch10,&
661 0 : 'Taking into account only s p d f '
662 0 : ABI_ERROR(msg)
663 : end if
664 :
665 50456 : a1=cryst%rprimd(:,1); b1=two_pi*Cryst%gprimd(:,1)
666 50456 : a2=cryst%rprimd(:,2); b2=two_pi*Cryst%gprimd(:,2)
667 50456 : a3=cryst%rprimd(:,3); b3=two_pi*Cryst%gprimd(:,3)
668 :
669 : ! Calculate Kleiman-Bylander factor and first derivative.
670 302556980 : fnl=czero_gw; fnld=czero_gw
671 :
672 2303630 : do ig=1,npw
673 : ! Get kcart = k+G in Cartesian coordinates.
674 9185688 : kg(:)= kpoint(:) + REAL(gvec(:,ig))
675 9185688 : kcart(:) = kg(1)*b1(:) + kg(2)*b2(:) + kg(3)*b3(:)
676 : ! Solve the problem with sinth=0. or sinphi=0
677 2296422 : if (ABS(kcart(2))<ppad) kcart(2) = kcart(2) + ppad
678 :
679 2296422 : mkg2 = kcart(1)**2+kcart(2)**2+kcart(3)**2
680 2296422 : mkg = SQRT(mkg2)
681 : ! The next to solve the problem with k=Gamma.
682 : !if (mkg < 0.0001) cycle
683 :
684 2296422 : sq=SQRT(kcart(1)**2+kcart(2)**2)
685 :
686 : gcart(:)= REAL(gvec(1,ig))*b1(:)&
687 : & +REAL(gvec(2,ig))*b2(:)&
688 9185688 : & +REAL(gvec(3,ig))*b3(:)
689 :
690 : ! Calculate spherical coordinates (th, phi).
691 2296422 : costh = kcart(3)/mkg
692 2296422 : sinth = sq/mkg
693 2296422 : cosphi= kcart(1)/sq
694 2296422 : sinphi= kcart(2)/sq
695 :
696 2296422 : gradth(1) = kcart(1)*kcart(3)/mkg**3/sinth
697 2296422 : gradth(2) = kcart(2)*kcart(3)/mkg**3/sinth
698 2296422 : gradth(3) = -(one/mkg-kcart(3)**2/mkg**3)/sinth
699 2296422 : gradphi(1) = -(one/sq - kcart(1)**2/sq**3)/sinphi
700 2296422 : gradphi(2) = kcart(2)*kcart(1)/sq**3/sinphi
701 2296422 : gradphi(3) = czero
702 :
703 6915877 : do iat=1,cryst%natom
704 4612247 : itypat = cryst%typat(iat)
705 4612247 : xdotg = gcart(1)*cryst%xcart(1,iat)+gcart(2)*Cryst%xcart(2,iat)+gcart(3)*Cryst%xcart(3,iat)
706 : ! Remember that in the GW code the reciprocal vectors
707 : ! are defined such as a_i*b_j = 2pi delta_ij, no need to introduce 2pi
708 4612247 : sfac=CMPLX(COS(xdotg), SIN(xdotg), kind=dp)
709 :
710 4612247 : iln0 = 0
711 22304549 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
712 24542383 : do ilmn=1,nlmn
713 17633714 : il = 1 + psps%indlmn(1,ilmn,itypat)
714 17633714 : in = psps%indlmn(3,ilmn,itypat)
715 17633714 : iln = psps%indlmn(5,ilmn,itypat)
716 : ! spin = 1 if scalar term (spin diagonal), 2 if SOC term.
717 : !spin = psps%indlmn(6, ilmn, itypat)
718 17633714 : if (iln <= iln0) cycle
719 17633714 : iln0 = iln
720 17633714 : if (vkbsign(iln,itypat) == zero) cycle
721 : !if (spin /= 1 .or. vkbsign(iln,itypat) == zero) cycle
722 17633714 : factor = SQRT(four_pi/REAL(2*(il-1)+1))
723 70067907 : do im=1,2*(il-1)+1
724 : ! Index of im and il
725 47821946 : ilm = im + (il-1)*(il-1)
726 :
727 : ! Calculate the first KB factor, note that fnl is simple precision complex
728 47821946 : fnl(ig,ilm,in,iat) = factor*sfac*ylmc(il-1,im-il,kcart) * vkb(ig,iln,itypat) * vkbsign(iln,itypat)
729 :
730 : ! Calculate the second KB factor (involving first derivatives)
731 : ! dYlm/dK = dYlm/dth * grad_K th + dYlm/dphi + grad_K phi
732 47821946 : call ylmcd(il-1,im-il,kcart,dth,dphi)
733 191287784 : dylmcart(:) = dth*gradth(:) + dphi*gradphi(:)
734 :
735 : ! Cartesian to crystallographic axis
736 : ! Notice: a bug was discovered by Marco Cazzaniga, december 2009
737 : ! the transformation matrix A=(a1,a2,a3) must act *on its left* on the
738 : ! covariant vector dylmcart (a *row* vector). The previous implementation assumed A
739 : ! acting on its right on a column vector, yielding wrong results for the (small)
740 : ! non local contributions to the spectra, such as a spurious anisotropy in isotropic systems.
741 : ! This is the correct version:
742 47821946 : dylmcrys(1) = (a1(1)*dylmcart(1)+a1(2)*dylmcart(2)+a1(3)*dylmcart(3))/(two_pi)
743 47821946 : dylmcrys(2) = (a2(1)*dylmcart(1)+a2(2)*dylmcart(2)+a2(3)*dylmcart(3))/(two_pi)
744 47821946 : dylmcrys(3) = (a3(1)*dylmcart(1)+a3(2)*dylmcart(2)+a3(3)*dylmcart(3))/(two_pi)
745 :
746 : ! Note that fnld is simple precision complex, it could be possible to use double precision
747 208921498 : do ii=1,3
748 : fnld(ii,ig,ilm,in,iat) = factor*sfac* &
749 191287784 : ( kg(ii)/mkg*ylmc(il-1,im-il,kcart)*vkbd(ig,iln,itypat) + dylmcrys(ii)*vkb(ig,iln,itypat) )
750 : end do
751 :
752 : end do !im
753 : end do !il
754 : end do !iat
755 : end do !ig
756 :
757 : DBG_EXIT("COLL")
758 :
759 7208 : end subroutine ccgradvnl_ylm
760 : !!***
761 :
762 : !----------------------------------------------------------------------
763 :
764 0 : END MODULE m_vkbr
765 : !!***
|
