Branch data 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 : 7210 : 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 : 7210 : real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:)
148 : : !************************************************************************
149 : :
150 : : !@vkbr_t
151 : 7210 : vkbr%istwfk = istwfk
152 : 7210 : vkbr%ntypat = cryst%ntypat
153 : 7210 : vkbr%natom = cryst%natom
154 : 7210 : vkbr%mpsang = psps%mpsang
155 : 7210 : vkbr%npw = npw
156 : 7210 : vkbr%inclvkb = inclvkb
157 [ + + ]: 28840 : 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 [ - + + - : 28840 : ABI_MALLOC(vkbsign, (psps%lnmax, cryst%ntypat))
- + + - +
- + - - +
- + ]
167 [ - + + - : 36050 : ABI_MALLOC(vkb, (npw, psps%lnmax, cryst%ntypat))
- + + - -
+ + - + -
+ - - + -
+ ]
168 [ - + + - : 28840 : ABI_MALLOC(vkbd, (npw, psps%lnmax, cryst%ntypat))
- + + - -
+ + - - +
- + ]
169 : 7210 : call calc_vkb(cryst,psps,kpoint,npw,npw,gvec,vkbsign,vkb,vkbd)
170 : :
171 : 7210 : select case (inclvkb)
172 : : case (2)
173 : : ! Complex spherical harmonics (CPU and mem \propto npw).
174 : 7210 : write(msg,'(a,f12.1)')'out-of-memory in fnl; Mb= ',one*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwp*b2Mb
175 [ - + + - : 43260 : ABI_STAT_MALLOC(vkbr%fnl,(npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr)
- + + - +
+ + - - +
+ - + - +
+ + - + -
- + + - ]
176 [ - + ]: 7210 : ABI_CHECK(ierr==0, msg)
177 : :
178 : 7210 : write(msg,'(a,f12.1)')'out-of-memory in fnld; Mb= ',three*npw*psps%mpsang**2*psps%mproj*cryst%natom*2*gwp*b2Mb
179 [ - + + - : 43260 : ABI_STAT_MALLOC(vkbr%fnld,(3,npw,psps%mpsang**2,psps%mproj,cryst%natom), ierr)
- + + - +
+ + - - +
+ - + - +
+ + - + -
- + + - ]
180 [ - + ]: 7210 : ABI_CHECK(ierr==0, msg)
181 : :
182 : 7210 : call ccgradvnl_ylm(cryst,psps,npw,gvec,kpoint,vkbsign,vkb,vkbd,vkbr%fnl,vkbr%fnld)
183 : :
184 : : case default
185 [ + - ]: 7210 : ABI_ERROR(sjoin("Wrong inclvkb= ",itoa(inclvkb)))
186 : : end select
187 : :
188 [ - + ]: 7210 : ABI_FREE(vkbsign)
189 [ - + ]: 7210 : ABI_FREE(vkb)
190 [ - + ]: 7210 : ABI_FREE(vkbd)
191 : :
192 : 7210 : 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 : 8164 : subroutine vkbr_free_0D(vkbr)
207 : :
208 : : !Arguments ------------------------------------
209 : : type(vkbr_t),intent(inout) :: vkbr
210 : : !************************************************************************
211 : :
212 : : !complex
213 [ + + ]: 8164 : ABI_SFREE(vkbr%fnl)
214 [ + + ]: 8164 : ABI_SFREE(vkbr%fnld)
215 : :
216 : 8164 : 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 [ - + ]: 148 : 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 [ + + ]: 1002 : do ii=1,SIZE(vkbr)
242 : 1002 : call vkbr_free_0D(vkbr(ii))
243 : : end do
244 : :
245 : 148 : 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 : 262147 : 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 [ - + ]: 262147 : ABI_CHECK(nspinor == 1, "inclvkb > 0 with nspinor == 2 is not coded")
306 : :
307 : : ! Adding term i <c,k|[Vnl,r]|v,k> ===
308 : 524294 : select case (vkbr%inclvkb)
309 : : case (2)
310 : : ! Complex spherical harmonics (much faster!).
311 : 262147 : dum=czero_gw; gamma_term=czero
312 : :
313 [ + + ]: 796947 : do iat=1,vkbr%natom
314 : 534800 : itypat = cryst%typat(iat)
315 [ + + + + ]: 2402434 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
316 : : iln0 = 0
317 [ + + ]: 2639687 : do ilmn=1,nlmn
318 : 1842740 : il = 1 + psps%indlmn(1,ilmn,itypat)
319 : 1842740 : in = psps%indlmn(3,ilmn,itypat)
320 : 1842740 : iln = psps%indlmn(5,ilmn,itypat)
321 [ + - ]: 1842740 : if (iln <= iln0) cycle
322 : 1842740 : iln0 = iln
323 : : !if (indlmn(6,ilmn,itypat) /= 1 .or. vkbsign(iln,itypat) == zero) cycle
324 : : !in = 1
325 [ + + ]: 7183924 : do im=1,2*(il-1)+1
326 : : ! Index of im and il
327 : 4806384 : ilm = im + (il-1)*(il-1)
328 : 4806384 : cta1 = czero_gw; cta2(:) = czero_gw
329 : 4806384 : cta4 = czero_gw; cta3(:) = czero_gw
330 [ + + ]: 2081861720 : do ig=1,npw
331 : : ! Here we take advantage of the property Y_(l-m)= (-i)^m Y_lm^*.
332 : 2077055336 : cta1 = cta1 + ug1(ig) * vkbr%fnl (ig,ilm,in,iat)
333 [ + + ]: 8308221344 : cta2(:)= cta2(:) + ug2(ig) * vkbr%fnld(:,ig,ilm,in,iat)
334 [ + + ]: 8308221344 : cta3(:)= cta3(:) + ug1(ig) * vkbr%fnld(:,ig,ilm,in,iat)
335 : 2077055336 : cta4 = cta4 + ug2(ig) * vkbr%fnl (ig,ilm,in,iat)
336 [ + + + + ]: 2096280872 : if (ig==1) gamma_term = gamma_term + CONJG(cta1)*cta2(:) +CONJG(cta3(:))*cta4
337 : : end do
338 [ + + ]: 21068276 : dum(:)= dum(:) + CONJG(cta1)*cta2(:) + CONJG(cta3(:))*cta4
339 : : end do
340 : :
341 : : end do
342 : : end do
343 : :
344 [ + + ]: 262147 : 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 [ + + ]: 1048588 : rhotwx(:,1) = rhotwx(:,1) + dum(:)
348 : :
349 : : case default
350 [ + - ]: 262147 : ABI_ERROR(sjoin("Wrong inclvkb:", itoa(vkbr%inclvkb)))
351 : : end select
352 : :
353 : 262147 : 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 : 7476 : 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 : 7476 : real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:),ylm_gr(:,:,:),ylm_k(:,:)
404 : : ! *************************************************************************
405 : :
406 : : DBG_ENTER("COLL")
407 [ - + ]: 7476 : ABI_CHECK(psps%usepaw==0, "You should not be here!")
408 [ - + ]: 7476 : ABI_CHECK(psps%useylm==0, "useylm/=0 not considered!")
409 : :
410 : : ! Compute KB dyadic sign.
411 [ + + + + ]: 43581 : vkbsign=zero
412 [ + + ]: 15032 : do itypat=1,psps%ntypat
413 : 7556 : iln0 = 0
414 [ + + + + ]: 36105 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
415 [ + + ]: 43486 : do ilmn=1,nlmn
416 : 28454 : iln = psps%indlmn(5,ilmn,itypat)
417 [ + - ]: 28454 : if (iln <= iln0) cycle
418 : 28454 : iln0 = iln
419 [ - + ]: 36010 : 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 : 7476 : idir=0; nkpg=0; ider=1; dimffnl=2 ! To retrieve the first derivative.
426 : :
427 : : ! Quantities used only if useylm==1
428 [ - + + - : 29904 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2*Psps%useylm))
+ - - - +
- - + - +
- + ]
429 [ - + + - : 37380 : ABI_MALLOC(ylm_gr, (npw_k, 3+6*(ider/2),psps%mpsang**2*Psps%useylm))
+ - + - -
- + - - +
- + - + ]
430 [ - + + - : 14952 : ABI_MALLOC(kpg_dum, (npw_k, nkpg))
- + - + ]
431 [ - + + - : 44856 : 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 : 7476 : psps%ntypat,Psps%pspso,Psps%qgrid_ff,cryst%rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr)
436 : :
437 [ - + ]: 7476 : ABI_FREE(ylm_k)
438 [ - + ]: 7476 : ABI_FREE(ylm_gr)
439 [ - + ]: 7476 : ABI_FREE(kpg_dum)
440 : :
441 [ - + + - : 22428 : ABI_MALLOC(modkplusg, (npw_k))
+ - + - -
+ - + ]
442 : 7476 : effmass_free = one; ecutsm = zero; ecut = huge(one)
443 : 7476 : call mkkin(ecut,ecutsm,effmass_free,cryst%gmet,kg_k,modkplusg,kpoint,npw_k,0,0)
444 [ + + ]: 2390156 : modkplusg(:) = SQRT(half/pi**2*modkplusg(:))
445 [ + + ]: 2390156 : modkplusg(:) = MAX(modkplusg(:),tol10)
446 : :
447 : : ! Calculate matrix elements.
448 [ + + + + : 21327942 : vkb=zero; vkbd=zero
+ + + + +
+ + + ]
449 : :
450 [ + + ]: 15032 : do itypat=1,psps%ntypat
451 : 7556 : iln0 = 0
452 [ + + + + ]: 36105 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
453 [ + + ]: 43486 : do ilmn=1,nlmn
454 : 28454 : il = 1 + psps%indlmn(1,ilmn,itypat)
455 : 28454 : in = psps%indlmn(3,ilmn,itypat)
456 : 28454 : iln = psps%indlmn(5,ilmn,itypat)
457 : : !write(*,*)ilmn, iln, il, in
458 [ + - ]: 28454 : if (iln <= iln0) cycle
459 : 28454 : iln0 = iln
460 : : !if (vkbsign(iln,itypat) == zero) cycle
461 [ + - ]: 36010 : if (ABS(psps%ekb(iln,itypat)) > 1.0d-10) then
462 [ - + ]: 28454 : ABI_CHECK(iln == ilmn, "iln != ilmn")
463 : : !ABI_CHECK(il == iln, "il != iln")
464 [ + + ]: 28454 : if (il==1) then
465 [ + + ]: 3568883 : vkb (1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)
466 [ + + ]: 3568883 : vkbd(1:npw_k,iln,itypat) = ffnl(:,2,iln,itypat)*modkplusg(:)/two_pi
467 [ + + ]: 18591 : else if (il==2) then
468 [ + + ]: 3946017 : vkb(1:npw_k,iln,itypat) = ffnl(:,1,iln,itypat)*modkplusg(:)
469 [ + + ]: 3946017 : do ig=1,npw_k
470 : : vkbd(ig,iln,itypat) = ((ffnl(ig,2,iln,itypat)*modkplusg(ig)*modkplusg(ig))+&
471 : 3946017 : 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 [ + + ]: 10620065 : vkb (:,iln,itypat) = SQRT(4*pi/cryst%ucvol*(2*il-1)*ABS(psps%ekb(iln,itypat)))*vkb (:,iln,itypat)
483 [ + + ]: 10620065 : 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 [ - + ]: 7476 : ABI_FREE(ffnl)
489 [ - + ]: 7476 : ABI_FREE(modkplusg)
490 : :
491 : : DBG_EXIT("COLL")
492 : :
493 : 7476 : 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 [ - + ]: 306253 : 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 [ + + ]: 2756277 : spinorwf_pad = reshape([0, 0, npw, npw, 0, npw, npw, 0], [2, 4])
565 [ + + + + ]: 1552001 : 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 [ + + ]: 306253 : if (istwfk == 1) then
570 [ + + ]: 598746 : do iab=1,nspinor
571 : 300237 : spad1 = spinorwf_pad(1,iab); spad2 = spinorwf_pad(2,iab)
572 [ + + ]: 95547266 : do ig=1,npw
573 : 94948520 : c_tmp = GWPC_CONJG(ug1(ig+spad1)) * ug2(ig+spad2)
574 [ + + ]: 380094317 : 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 [ + + ]: 306253 : if (inclvkb /= 0) then
588 [ - + ]: 262147 : ABI_CHECK(istwfk == vkbr%istwfk, "input istwfk /= vkbr%istwfk")
589 : 262147 : call add_vnlr_commutator(vkbr,cryst,psps,npw,nspinor,ug1,ug2,ihr_comm)
590 : : end if
591 : :
592 : 306253 : 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 : 7210 : 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 [ - + ]: 7210 : 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 [ + + + + ]: 50470 : a1=cryst%rprimd(:,1); b1=two_pi*Cryst%gprimd(:,1)
666 [ + + + + ]: 50470 : a2=cryst%rprimd(:,2); b2=two_pi*Cryst%gprimd(:,2)
667 [ + + + + ]: 50470 : a3=cryst%rprimd(:,3); b3=two_pi*Cryst%gprimd(:,3)
668 : :
669 : : ! Calculate Kleiman-Bylander factor and first derivative.
670 [ + + + + : 302572550 : fnl=czero_gw; fnld=czero_gw
+ + + + +
+ + + + +
+ + + + ]
671 : :
672 [ + + ]: 2303804 : do ig=1,npw
673 : : ! Get kcart = k+G in Cartesian coordinates.
674 [ + + ]: 9186376 : kg(:)= kpoint(:) + REAL(gvec(:,ig))
675 [ + + ]: 9186376 : kcart(:) = kg(1)*b1(:) + kg(2)*b2(:) + kg(3)*b3(:)
676 : : ! Solve the problem with sinth=0. or sinphi=0
677 [ + + ]: 2296594 : if (ABS(kcart(2))<ppad) kcart(2) = kcart(2) + ppad
678 : :
679 : 2296594 : mkg2 = kcart(1)**2+kcart(2)**2+kcart(3)**2
680 : 2296594 : mkg = SQRT(mkg2)
681 : : ! The next to solve the problem with k=Gamma.
682 : : !if (mkg < 0.0001) cycle
683 : :
684 : 2296594 : sq=SQRT(kcart(1)**2+kcart(2)**2)
685 : :
686 : : gcart(:)= REAL(gvec(1,ig))*b1(:)&
687 : : & +REAL(gvec(2,ig))*b2(:)&
688 [ + + ]: 9186376 : & +REAL(gvec(3,ig))*b3(:)
689 : :
690 : : ! Calculate spherical coordinates (th, phi).
691 : 2296594 : costh = kcart(3)/mkg
692 : 2296594 : sinth = sq/mkg
693 : 2296594 : cosphi= kcart(1)/sq
694 : 2296594 : sinphi= kcart(2)/sq
695 : :
696 : 2296594 : gradth(1) = kcart(1)*kcart(3)/mkg**3/sinth
697 : 2296594 : gradth(2) = kcart(2)*kcart(3)/mkg**3/sinth
698 : 2296594 : gradth(3) = -(one/mkg-kcart(3)**2/mkg**3)/sinth
699 : 2296594 : gradphi(1) = -(one/sq - kcart(1)**2/sq**3)/sinphi
700 : 2296594 : gradphi(2) = kcart(2)*kcart(1)/sq**3/sinphi
701 : 2296594 : gradphi(3) = czero
702 : :
703 [ + + ]: 6916395 : do iat=1,cryst%natom
704 : 4612591 : itypat = cryst%typat(iat)
705 : 4612591 : 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 : 4612591 : sfac=CMPLX(COS(xdotg), SIN(xdotg), kind=dp)
709 : :
710 : 4612591 : iln0 = 0
711 [ + + + + ]: 22305581 : nlmn = count(psps%indlmn(3,:,itypat) > 0)
712 [ + + ]: 24543415 : do ilmn=1,nlmn
713 : 17634230 : il = 1 + psps%indlmn(1,ilmn,itypat)
714 : 17634230 : in = psps%indlmn(3,ilmn,itypat)
715 : 17634230 : 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 [ - + ]: 17634230 : if (iln <= iln0) cycle
719 : 17634230 : iln0 = iln
720 [ + - ]: 17634230 : if (vkbsign(iln,itypat) == zero) cycle
721 : : !if (spin /= 1 .or. vkbsign(iln,itypat) == zero) cycle
722 : 17634230 : factor = SQRT(four_pi/REAL(2*(il-1)+1))
723 [ + + ]: 70069627 : do im=1,2*(il-1)+1
724 : : ! Index of im and il
725 : 47822806 : ilm = im + (il-1)*(il-1)
726 : :
727 : : ! Calculate the first KB factor, note that fnl is simple precision complex
728 : 47822806 : 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 : 47822806 : call ylmcd(il-1,im-il,kcart,dth,dphi)
733 [ + + ]: 191291224 : 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 : 47822806 : dylmcrys(1) = (a1(1)*dylmcart(1)+a1(2)*dylmcart(2)+a1(3)*dylmcart(3))/(two_pi)
743 : 47822806 : dylmcrys(2) = (a2(1)*dylmcart(1)+a2(2)*dylmcart(2)+a2(3)*dylmcart(3))/(two_pi)
744 : 47822806 : 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 [ + + ]: 208925454 : do ii=1,3
748 : : fnld(ii,ig,ilm,in,iat) = factor*sfac* &
749 : 191291224 : ( 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 : 7210 : end subroutine ccgradvnl_ylm
760 : : !!***
761 : :
762 : : !----------------------------------------------------------------------
763 : :
764 [ # # # # : 0 : END MODULE m_vkbr
# # # # #
# # # #
# ]
765 : : !!***
|
