Line data Source code
1 : !!****m* ABINIT/m_melemts
2 : !! NAME
3 : !! m_melemts
4 : !!
5 : !! FUNCTION
6 : !! This module defines an object used as database to store matrix
7 : !! elements of several potentials and operators between two Bloch states.
8 : !! These values are used in the GW part of abinit to evaluate QP energies
9 : !! using the perturbative approach.
10 : !!
11 : !! COPYRIGHT
12 : !! Copyright (C) 2008-2026 ABINIT group (MG)
13 : !! This file is distributed under the terms of the
14 : !! GNU General Public License, see ~abinit/COPYING
15 : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : !!
17 : !! TODO
18 :
19 : !! NOTES
20 : !! * This module is supposed to be used only in the GW part to facilitate
21 : !! we might decide to use ragged arrays
22 : !!
23 : !! Mels(nkcalc, nsppol*nspinor**2)%data
24 : !!
25 : !! or replaced nkibz with nkcalc with reduce memory
26 : !!
27 : !! SOURCE
28 :
29 : #if defined HAVE_CONFIG_H
30 : #include "config.h"
31 : #endif
32 :
33 : #include "abi_common.h"
34 :
35 : module m_melemts
36 :
37 : use defs_basis
38 : use m_errors
39 : use m_xmpi
40 : use m_abicore
41 :
42 : use m_fstrings, only : tolower, sjoin
43 : use m_numeric_tools, only : print_arr
44 :
45 : implicit none
46 :
47 : private
48 : !!***
49 :
50 : !----------------------------------------------------------------------
51 :
52 : !!****t* m_melemts/melflags_t
53 : !! NAME
54 : !!
55 : !! FUNCTION
56 : !! Container for the flags defining the status of the corresponding
57 : !! pointer defined in the type melements_t. Possible values are:
58 : !! * 0 if the correspondending array is not allocated.
59 : !! * 1 if allocated but not yet calculated.
60 : !! * 2 if allocated and calculated
61 : !!
62 : !! SOURCE
63 :
64 : type,public :: melflags_t
65 :
66 : integer :: has_kinetic = 0
67 : integer :: has_hbare = 0
68 : integer :: has_lexexch = 0
69 : integer :: has_sxcore = 0
70 : integer :: has_vhartree = 0
71 : integer :: has_vu = 0
72 : integer :: has_vxc = 0
73 : integer :: has_vxcval = 0
74 : integer :: has_vxcval_hybrid = 0
75 : integer :: only_diago = 0
76 : ! 1 if only diagonal elements are calculated
77 :
78 : contains
79 : procedure :: reset => melflags_reset ! Reset the value of the flags.
80 : procedure :: copy => melflags_copy ! Copy the object
81 : end type melflags_t
82 : !!***
83 :
84 : !----------------------------------------------------------------------
85 :
86 : !!****t* m_melemts/melements_t
87 : !! NAME
88 : !!
89 : !! FUNCTION
90 : !! Structure defining a database to store the matrix elements of operators
91 : !! needed for GW calculations.
92 : !!
93 : !! SOURCE
94 :
95 : type,public :: melements_t
96 :
97 : integer :: nkibz
98 : ! Number of k-points in the IBZ.
99 :
100 : integer :: nsppol
101 : ! Number of independent spin-polarizations.
102 :
103 : integer :: nspinor
104 : ! 1 for collinear, 2 for noncollinear.
105 :
106 : integer :: nspden
107 : ! Number of independent spin-density components.
108 :
109 : integer :: bmin, bmax
110 : ! min and Max band index over k-points and spin.
111 : ! Used to dimension the arrays below.
112 :
113 : integer, allocatable :: bands_idx(:,:,:)
114 : ! (2, nkibz, nsppol)
115 : ! min and Max band index for each k-point and spin.
116 :
117 : integer, allocatable :: iscalc(:,:)
118 : ! (nkibz, nsppol)
119 : ! 1 if this k-point and spin has been calculated, 0 otherwise.
120 :
121 : real(dp), allocatable :: kibz(:,:)
122 : ! (3, nkibz)
123 : ! The list of k-points in reduced coordinates.
124 :
125 : complex(dp), allocatable :: kinetic(:,:,:,:)
126 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
127 : ! Matrix elements of the kinetic energy.
128 :
129 : complex(dp), allocatable :: hbare(:,:,:,:)
130 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
131 : ! Matrix elements of the bare Hamiltonian.
132 :
133 : complex(dp), allocatable :: sxcore(:,:,:,:)
134 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
135 : ! Matrix elements of the Fock operator generated by core electrons.
136 :
137 : complex(dp), allocatable :: vhartree(:,:,:,:)
138 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
139 : ! Matrix elements of the Hartree potential.
140 :
141 : complex(dp), allocatable :: vlexx(:,:,:,:)
142 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
143 : ! Matrix elements of the local exact exchange potential.
144 :
145 : complex(dp), allocatable :: vu(:,:,:,:)
146 : ! vu(b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
147 : ! Matrix elements of the U Hamiltonian.
148 :
149 : complex(dp), allocatable :: vxc(:,:,:,:)
150 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
151 : ! Matrix elements of XC potential, including model core if present.
152 :
153 : complex(dp), allocatable :: vxcval(:,:,:,:)
154 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
155 : ! Matrix elements of the XC potential, valence-only contribution.
156 :
157 : complex(dp), allocatable :: vxcval_hybrid(:,:,:,:)
158 : ! (b1:b2, b1:b2, nkibz, nsppol*nspinor**2)
159 : ! Matrix elements of the XC potential for hybrid calculations, valence-only contribution.
160 :
161 : type(melflags_t) :: flags
162 :
163 : contains
164 :
165 : procedure :: init => melements_init ! Initialize the object
166 : procedure :: free => melements_free ! Free memory
167 : procedure :: herm => melements_herm ! Construct the lower triangle from the upper triangle
168 : procedure :: mpisum => melements_mpisum ! Perform a collective SUM within the MPI communicator comm
169 : procedure :: print => melements_print ! Print matrix elements
170 : procedure :: zero => melements_zero ! Set matrix elements connecting states with different irrep to zero.
171 : !procedure :: mels_get_exene_core
172 :
173 : end type melements_t
174 : !!***
175 :
176 : integer,parameter,private :: NNAMES = 9
177 : integer,parameter,private :: NAMELEN = 13
178 :
179 : ! List of matrix element names.
180 : ! IMPORTANT: Remember to add a new entry here when you add a new matrix element.
181 : ! then change the routine my_select_melements accordingly.
182 : character(len=NAMELEN),parameter,private :: ANAMES(NNAMES) = [ &
183 : "vxc ", &
184 : "vxcval ", &
185 : "vxcval_hybrid", &
186 : "sxcore ", &
187 : "vu ", &
188 : "vlexx ", &
189 : "vhartree ", &
190 : "kinetic ", &
191 : "hbare " &
192 : ]
193 :
194 : CONTAINS !========================================================================================
195 :
196 : !!****f* m_melemts/melflags_reset
197 : !! NAME
198 : !! melflags_reset
199 : !!
200 : !! FUNCTION
201 : !! Set all flags in melflags_t to 0.
202 : !!
203 839 : subroutine melflags_reset(Mflags)
204 :
205 : !Arguments ------------------------------------
206 : class(melflags_t),intent(inout) :: Mflags
207 : ! *************************************************************************
208 :
209 839 : Mflags%has_kinetic = 0
210 839 : Mflags%has_hbare = 0
211 839 : Mflags%has_sxcore = 0
212 839 : Mflags%has_vhartree = 0
213 839 : Mflags%has_vu = 0
214 839 : Mflags%has_vxc = 0
215 839 : Mflags%has_vxcval = 0
216 839 : Mflags%has_vxcval_hybrid = 0
217 839 : Mflags%has_lexexch = 0
218 839 : Mflags%only_diago = 0
219 :
220 839 : end subroutine melflags_reset
221 : !!***
222 :
223 : !----------------------------------------------------------------------
224 :
225 : !!****f* m_melemts/melflags_copy
226 : !! NAME
227 : !! melflags_copy
228 : !!
229 : !! FUNCTION
230 : !! Copy an object storing the flags.
231 : !!
232 : !! INPUTS
233 : !! Mflags_in=The flags to be copied.
234 : !!
235 : !! OUTPUT
236 : !! Mflags_out=The new set of flags.
237 : !!
238 : !! SOURCE
239 :
240 281 : subroutine melflags_copy(Mflags_in, Mflags_out)
241 :
242 : !Arguments ------------------------------------
243 : class(melflags_t),intent(in) :: Mflags_in
244 : class(melflags_t),intent(inout) :: Mflags_out
245 : ! *************************************************************************
246 :
247 281 : call Mflags_out%reset()
248 :
249 : ! @melflags_t
250 281 : Mflags_out%has_kinetic = Mflags_in%has_kinetic
251 281 : Mflags_out%has_hbare = Mflags_in%has_hbare
252 281 : Mflags_out%has_sxcore = Mflags_in%has_sxcore
253 281 : Mflags_out%has_vhartree = Mflags_in%has_vhartree
254 281 : Mflags_out%has_vu = Mflags_in%has_vu
255 281 : Mflags_out%has_vxc = Mflags_in%has_vxc
256 281 : Mflags_out%has_vxcval = Mflags_in%has_vxcval
257 281 : Mflags_out%has_vxcval_hybrid = Mflags_in%has_vxcval_hybrid
258 281 : Mflags_out%has_lexexch = Mflags_in%has_lexexch
259 281 : Mflags_out%only_diago = Mflags_in%only_diago
260 :
261 281 : end subroutine melflags_copy
262 : !!***
263 :
264 : !----------------------------------------------------------------------
265 :
266 : !!****f* m_melemts/melements_free
267 : !! NAME
268 : !! melements_free
269 : !!
270 : !! FUNCTION
271 : !! Free all dynamic memory of the database
272 : !!
273 : !! SOURCE
274 :
275 281 : subroutine melements_free(Mels)
276 :
277 : !Arguments ------------------------------------
278 : class(melements_t),intent(inout) :: Mels
279 : ! *************************************************************************
280 :
281 : ! integer arrays
282 281 : ABI_SFREE(Mels%bands_idx)
283 281 : ABI_SFREE(Mels%iscalc)
284 :
285 : ! real arrays
286 281 : ABI_SFREE(Mels%kibz)
287 :
288 : ! complex arrays
289 281 : ABI_SFREE(Mels%kinetic)
290 281 : ABI_SFREE(Mels%hbare)
291 281 : ABI_SFREE(Mels%sxcore)
292 281 : ABI_SFREE(Mels%vhartree)
293 281 : ABI_SFREE(Mels%vlexx)
294 281 : ABI_SFREE(Mels%vu)
295 281 : ABI_SFREE(Mels%vxc)
296 281 : ABI_SFREE(Mels%vxcval)
297 281 : ABI_SFREE(Mels%vxcval_hybrid)
298 :
299 : ! Reset all has_* flags.
300 281 : call Mels%flags%reset()
301 :
302 281 : end subroutine melements_free
303 : !!***
304 :
305 : !----------------------------------------------------------------------
306 :
307 : !!****f* m_melemts/my_select_melements
308 : !! NAME
309 : !! my_select_melements
310 : !!
311 : !! FUNCTION
312 : !! Helper function returning a pointer to the array "aname" as well as the status of the array
313 : !!
314 : !! INPUTS
315 : !! Mels<melements_t>=The database.
316 : !! aname=String with the name of the array.
317 : !!
318 : !! OUTPUT
319 : !! flag_p=Pointer to the integer defining the status of the array, see melflags_t.
320 : !! arr_p=The pointer to the array.
321 : !!
322 : !! SOURCE
323 :
324 7551 : subroutine my_select_melements(Mels, aname, flag_p, arr_p)
325 :
326 : !Arguments ------------------------------------
327 : !scalars
328 : integer,pointer :: flag_p
329 : character(len=*),intent(in) :: aname
330 : type(melements_t),target,intent(in) :: Mels
331 : !arrays
332 : complex(dp),contiguous, pointer :: arr_p(:,:,:,:)
333 : ! *************************************************************************
334 :
335 15941 : SELECT CASE (tolower(aname))
336 : CASE ("kinetic")
337 839 : flag_p => Mels%flags%has_kinetic
338 839 : arr_p => Mels%kinetic
339 : CASE ("hbare")
340 839 : flag_p => Mels%flags%has_hbare
341 839 : arr_p => Mels%hbare
342 : CASE ("sxcore")
343 839 : flag_p => Mels%flags%has_sxcore
344 839 : arr_p => Mels%sxcore
345 : CASE ("vhartree")
346 839 : flag_p => Mels%flags%has_vhartree
347 839 : arr_p => Mels%vhartree
348 : CASE ("vlexx")
349 839 : flag_p => Mels%flags%has_lexexch
350 839 : arr_p => Mels%vlexx
351 : CASE ("vu")
352 839 : flag_p => Mels%flags%has_vu
353 839 : arr_p => Mels%vu
354 : CASE ("vxc")
355 839 : flag_p => Mels%flags%has_vxc
356 839 : arr_p => Mels%vxc
357 : CASE ("vxcval")
358 839 : flag_p => Mels%flags%has_vxcval
359 839 : arr_p => Mels%vxcval
360 : CASE ("vxcval_hybrid")
361 839 : flag_p => Mels%flags%has_vxcval_hybrid
362 839 : arr_p => Mels%vxcval_hybrid
363 : CASE DEFAULT
364 7551 : ABI_ERROR(sjoin("Wrong aname: ", aname))
365 : END SELECT
366 :
367 7551 : end subroutine my_select_melements
368 : !!***
369 :
370 : !----------------------------------------------------------------------
371 :
372 : !!****f* m_melemts/melements_init
373 : !! NAME
374 : !! melements_init
375 : !!
376 : !! FUNCTION
377 : !! Initialize the database, allocate arrays according to
378 : !! Mflags_in, zeroing the content of the allocated arrays.
379 : !!
380 : !! INPUTS
381 : !! nsppol=Number of independent spin polarizations.
382 : !! nspden=Number of spin-density components
383 : !! nspinor=Number of spinor components
384 : !! band_idx=min and Max band index for each ik_ibz and spin
385 : !!
386 : !! OUTPUT
387 : !! Mels=The initialized database with dimensions and allocated memory.
388 : !!
389 : !! SOURCE
390 :
391 281 : subroutine melements_init(Mels, Mflags_in, nsppol, nspden, nspinor, nkibz, kibz, bands_idx)
392 :
393 : !Arguments ------------------------------------
394 : !scalars
395 : class(melements_t),intent(out) :: Mels
396 : integer,intent(in) :: nspinor,nspden,nsppol,nkibz
397 : type(melflags_t),intent(in) :: Mflags_in
398 : !arrays
399 : integer,intent(in) :: bands_idx(2, nkibz, nsppol)
400 : real(dp),intent(in) :: kibz(3,nkibz)
401 :
402 : !Local variables-------------------------------
403 : integer :: ikibz,isppol,bmin,bmax,b1,b2
404 : ! *************************************************************************
405 :
406 : ! Copy flags.
407 281 : call Mflags_in%copy(Mels%flags)
408 :
409 : ! Copy dimensions.
410 281 : Mels%nkibz = nkibz
411 281 : Mels%nsppol = nsppol
412 281 : Mels%nspinor = nspinor
413 281 : Mels%nspden = nspden
414 :
415 1124 : ABI_MALLOC(Mels%bands_idx, (2, nkibz, nsppol))
416 5976 : Mels%bands_idx = bands_idx
417 :
418 1124 : ABI_MALLOC(Mels%iscalc,(nkibz, nsppol))
419 2277 : Mels%iscalc = 0
420 :
421 281 : bmin = HUGE(1); bmax =-HUGE(1)
422 568 : do isppol=1,Mels%nsppol
423 2277 : do ikibz=1,Mels%nkibz
424 3274 : if (ANY(Mels%bands_idx(:,ikibz,isppol)/=0)) then
425 1070 : b1 = Mels%bands_idx(1,ikibz,isppol)
426 1070 : b2 = Mels%bands_idx(2,ikibz,isppol)
427 1070 : Mels%iscalc(ikibz, isppol)=1
428 1070 : bmin = MIN(bmin,b1)
429 1070 : bmax = MAX(bmax,b2)
430 1070 : ABI_CHECK(b2 >= b1 .and. b1 > 0, "Wrong b1, b2")
431 : end if
432 : end do
433 : end do
434 :
435 281 : if (bmin==HUGE(1).or.bmax==-HUGE(1)) then
436 0 : ABI_BUG("Wrong bands_idx")
437 : end if
438 :
439 281 : Mels%bmin = bmin
440 281 : Mels%bmax = bmax
441 :
442 281 : b1 = Mels%bmin; b2 = Mels%bmax
443 :
444 : ! real arrays
445 843 : ABI_MALLOC(Mels%kibz, (3,nkibz))
446 7350 : Mels%kibz = kibz
447 :
448 : ! complex arrays
449 281 : if (Mels%flags%has_kinetic == 1) then
450 3909 : ABI_CALLOC(Mels%kinetic, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
451 : end if
452 281 : if (Mels%flags%has_hbare == 1) then
453 3909 : ABI_CALLOC(Mels%hbare, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
454 : end if
455 281 : if (Mels%flags%has_sxcore == 1) then
456 1090 : ABI_CALLOC(Mels%sxcore, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
457 : end if
458 281 : if (Mels%flags%has_vhartree == 1) then
459 160810 : ABI_CALLOC(Mels%vhartree, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
460 : end if
461 281 : if (Mels%flags%has_lexexch == 1) then
462 0 : ABI_CALLOC(Mels%vlexx, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
463 : end if
464 281 : if (Mels%flags%has_vu == 1) then
465 0 : ABI_CALLOC(Mels%vu, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
466 : end if
467 281 : if (Mels%flags%has_vxc == 1) then
468 107903 : ABI_CALLOC(Mels%vxc, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
469 : end if
470 281 : if (Mels%flags%has_vxcval == 1) then
471 107903 : ABI_CALLOC(Mels%vxcval, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
472 : end if
473 281 : if (Mels%flags%has_vxcval_hybrid == 1) then
474 28278 : ABI_CALLOC(Mels%vxcval_hybrid, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
475 : end if
476 :
477 281 : end subroutine melements_init
478 : !!***
479 :
480 : !----------------------------------------------------------------------
481 :
482 : !!****f* m_melemts/melements_herm
483 : !! NAME
484 : !! melements_herm
485 : !!
486 : !! FUNCTION
487 : !! Reconstruc the lower triangle of all calculated arrays.
488 : !! Assuming Hermitian operator. Works both for collinear and non-collinear case.
489 : !!
490 : !! INPUTS
491 : !! Mels=The database
492 : !! [aname]=The name of the array to be symmetrized, by default all calculated arrays are filled.
493 : !!
494 : !! SIDE EFFECTS
495 : !! All arrays whose flag is 2, are filled assuming an Hermitian operator.
496 : !!
497 : !! SOURCE
498 :
499 281 : subroutine melements_herm(Mels, aname)
500 :
501 : !Arguments ------------------------------------
502 : !scalars
503 : class(melements_t),intent(inout) :: Mels
504 : character(len=*),optional,intent(in) :: aname
505 :
506 : !Local variables-------------------------------
507 : integer :: is,ik,ib,jb,iab,iab_tr,iname
508 : integer,pointer :: flag_p
509 : character(len=NAMELEN) :: key
510 : !arrays
511 : integer,parameter :: trsp_idx(2:4) = [2,4,3]
512 281 : complex(dp),contiguous, pointer :: arr_p(:,:,:,:)
513 : ! *************************************************************************
514 :
515 : ! === Symmetrize matrix elements ===
516 : ! * In the collinear case, generate the lower triangle by just doing a complex conjugate.
517 : ! * In the noncollinear case do also a transposition since A_{12}^{ab} = A_{21}^{ba}^*
518 : ! 2-->2, 3-->4, 4-->3
519 : !
520 2810 : do iname=1,NNAMES
521 2529 : key = ANAMES(iname)
522 2529 : if (PRESENT(aname)) then
523 0 : if (key /= aname) CYCLE
524 : end if
525 :
526 2529 : call my_select_melements(Mels, key, flag_p, arr_p)
527 :
528 2810 : if (flag_p>0) then
529 5230 : do ik=1,Mels%nkibz
530 9759 : do is=1,Mels%nsppol
531 :
532 43495 : do jb=Mels%bmin,Mels%bmax
533 215484 : do ib=Mels%bmin,jb ! Upper triangle
534 :
535 210955 : if (ib/=jb) then
536 142017 : arr_p(jb,ib,ik,is)=CONJG(arr_p(ib,jb,ik,is))
537 142017 : if (Mels%nspinor==2) then
538 29052 : do iab=2,4
539 21789 : iab_tr=trsp_idx(iab)
540 29052 : arr_p(jb,ib,ik,iab)=CONJG(arr_p(ib,jb,ik,iab_tr))
541 : end do
542 : end if
543 : else ! For ib==jb force real-valued
544 34469 : arr_p(jb,ib,ik,is)=half*(arr_p(jb,ib,ik,is)+CONJG(arr_p(jb,ib,ik,is)))
545 34469 : if (Mels%nspinor==2) arr_p(jb,ib,ik,2)=half*(arr_p(ib,jb,ik,2)+CONJG(arr_p(ib,jb,ik,2)))
546 : end if
547 :
548 : end do !ib
549 : end do !jb
550 :
551 : end do !is
552 : end do !ik
553 : end if
554 :
555 : end do !inames
556 :
557 281 : end subroutine melements_herm
558 : !!***
559 :
560 : !----------------------------------------------------------------------
561 :
562 : !!****f* m_melemts/melements_mpisum
563 : !! NAME
564 : !! melements_mpisum
565 : !!
566 : !! FUNCTION
567 : !! Perform a collective SUM within the MPI communicator comm
568 : !! of the matrix elements stored in the database.
569 : !!
570 : !! INPUTS
571 : !! Mels=The database
572 : !! [aname]=The name of a particular array to be summed, by default
573 : !! all allocated arrays are considered.
574 : !!
575 : !! SIDE EFFECTS
576 : !! All arrays whose flag==1 are summed within the MPI communicator comm.
577 : !! In output the corresponding flas is set to 2.
578 : !!
579 : !! SOURCE
580 :
581 281 : subroutine melements_mpisum(Mels, comm, aname)
582 :
583 : !Arguments ------------------------------------
584 : !scalars
585 : class(melements_t),intent(inout) :: Mels
586 : integer,intent(in) :: comm
587 : character(len=*),optional,intent(in) :: aname
588 :
589 : !Local variables-------------------------------
590 : integer :: iname,ierr
591 : integer,pointer :: flag_p
592 : character(len=NAMELEN) :: key
593 : !character(len=500) :: msg
594 : !arrays
595 281 : complex(dp),contiguous, pointer :: arr_p(:,:,:,:)
596 : ! *************************************************************************
597 :
598 2810 : do iname=1,NNAMES
599 2529 : key = ANAMES(iname)
600 2529 : if (PRESENT(aname)) then
601 0 : if (key /= aname) CYCLE
602 : end if
603 :
604 2529 : call my_select_melements(Mels, key, flag_p, arr_p)
605 :
606 2810 : if (flag_p == 1) then
607 733 : call xmpi_sum(arr_p, comm, ierr)
608 733 : flag_p = 2 ! Tag this array as calculated
609 : end if
610 : end do
611 :
612 281 : end subroutine melements_mpisum
613 : !!***
614 :
615 : !----------------------------------------------------------------------
616 :
617 : !!****f* m_melemts/melements_print
618 : !! NAME
619 : !! melements_print
620 : !!
621 : !! FUNCTION
622 : !! Printout of the content of all calculated array.
623 : !! Optionally, it is possible to print the content of a single entry of the database.
624 : !!
625 : !! INPUTS
626 : !! Mels=The database
627 : !! [unit]=the unit number for output, defaults to std_out
628 : !! [prtvol]=verbosity level, defaults to 0
629 : !! [mode_paral]=either "COLL" or "PERS", default to "COLL"
630 : !! [header]=title for info
631 : !!
632 : !! OUTPUT
633 : !! Only writing
634 : !!
635 : !! SOURCE
636 :
637 0 : subroutine melements_print(Mels, names_list, header, unit, prtvol, mode_paral)
638 :
639 : !Arguments ------------------------------------
640 : !scalars
641 : class(melements_t),intent(in) :: Mels
642 : integer,optional,intent(in) :: prtvol,unit
643 : character(len=*),optional,intent(in) :: names_list(:)
644 : character(len=*),optional,intent(in) :: header
645 : character(len=4),optional,intent(in) :: mode_paral
646 :
647 : !Local variables-------------------------------
648 : integer :: my_unt,my_prtvol,max_r,max_c,ii, isppol,ikibz,iab,ib,b1,b2,my_nkeys,ikey
649 : integer,pointer :: flag_p
650 : character(len=4) :: my_mode
651 : character(len=NAMELEN) :: key
652 : character(len=500) :: msg,str,fmt
653 : !arrays
654 277 : integer,allocatable :: tab(:)
655 277 : character(len=NAMELEN),allocatable :: my_keys(:)
656 277 : complex(dp),allocatable :: mat(:,:)
657 :
658 : type rarr_dpc4
659 : complex(dp),contiguous, pointer :: arr_p(:,:,:,:)
660 : end type rarr_dpc4
661 277 : type(rarr_dpc4),allocatable :: data_p(:)
662 : ! *************************************************************************
663 :
664 : !@melements_t
665 277 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
666 277 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
667 277 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
668 :
669 277 : if (Mels%nspinor == 2) ABI_WARNING("nspinor=2 not coded")
670 :
671 277 : if (present(names_list)) then
672 0 : my_nkeys=SIZE(names_list)
673 0 : ABI_MALLOC(my_keys, (my_nkeys))
674 0 : my_keys = names_list
675 : else
676 277 : my_nkeys = NNAMES
677 277 : ABI_MALLOC(my_keys, (NNAMES))
678 3047 : my_keys = ANAMES
679 : end if
680 :
681 831 : ABI_MALLOC(data_p, (my_nkeys))
682 831 : ABI_MALLOC(tab, (my_nkeys))
683 2770 : tab = 0
684 :
685 277 : my_nkeys=0; str = " ib"; ii=4
686 2770 : do ikey=1,size(my_keys)
687 2493 : key = my_keys(ikey)
688 2493 : call my_select_melements(Mels,key,flag_p,data_p(ikey)%arr_p)
689 2770 : if (flag_p == 2) then
690 713 : my_nkeys = my_nkeys+1
691 713 : tab(my_nkeys) = ikey
692 713 : str(ii+1:)=" "//TRIM(tolower(key))
693 713 : ii = ii+MAX(1+LEN_TRIM(key),10)
694 713 : ABI_CHECK(ii <490, "I'm gonna SIGFAULT!")
695 : end if
696 : end do
697 :
698 277 : write(msg,'(2a)')ch10,' === Matrix Elements stored in Mels% [eV] === '
699 277 : if (PRESENT(header)) write(msg,'(4a)')ch10,' === '//TRIM(ADJUSTL(header))//' [eV] === '
700 277 : call wrtout(my_unt,msg,my_mode)
701 277 : if (Mels%nspinor == 2) then
702 7 : call wrtout(my_unt, "Sum_ab M_ab, M_11, M_22, Re(M_12), IM(Re_12)" ,my_mode)
703 : end if
704 :
705 277 : if (my_nkeys==0) GOTO 10
706 277 : write(fmt,'(a,i4,a)')'(1x,i3,',my_nkeys,'(1x,f9.5))' ! width of 10 chars
707 :
708 560 : do isppol=1,Mels%nsppol
709 2245 : do ikibz=1,Mels%nkibz
710 1685 : if (Mels%iscalc(ikibz,isppol)/=1) CYCLE
711 :
712 1046 : write(msg,'(a,3es16.8,a,i2,a)')" kpt= (",Mels%kibz(:,ikibz),") spin=",isppol,":"
713 1046 : call wrtout(my_unt,msg,my_mode)
714 :
715 1046 : b1 = Mels%bands_idx(1,ikibz,isppol)
716 1046 : b2 = Mels%bands_idx(2,ikibz,isppol)
717 :
718 1329 : if (Mels%flags%only_diago==1 .or. my_prtvol==0) then
719 : ! Print only the diagonal.
720 814 : write(msg,'(a)')str
721 814 : call wrtout(my_unt,msg,my_mode)
722 7195 : do ib=b1,b2
723 7195 : if (Mels%nspinor == 1) then
724 20888 : write(msg,fmt)ib,(REAL(data_p(tab(ikey))%arr_p(ib,ib,ikibz,1))*Ha_eV, ikey=1,my_nkeys)
725 6237 : call wrtout(my_unt,msg,my_mode)
726 : else
727 : ! Write sum_ab, then diagonal elements, finally Re_12, Im_12
728 2304 : write(msg,fmt)ib,(real(sum(data_p(tab(ikey))%arr_p(ib,ib,ikibz,:)))*Ha_eV, ikey=1,my_nkeys)
729 144 : call wrtout(my_unt,msg,my_mode)
730 144 : if (my_prtvol > 0) then
731 0 : write(msg,fmt)ib,(real(data_p(tab(ikey))%arr_p(ib,ib,ikibz,1))*Ha_eV, ikey=1,my_nkeys)
732 0 : call wrtout(my_unt,msg,my_mode)
733 0 : write(msg,fmt)ib,(real(data_p(tab(ikey))%arr_p(ib,ib,ikibz,2))*Ha_eV, ikey=1,my_nkeys)
734 0 : call wrtout(my_unt,msg,my_mode)
735 0 : write(msg,fmt)ib,(real(data_p(tab(ikey))%arr_p(ib,ib,ikibz,3))*Ha_eV, ikey=1,my_nkeys)
736 0 : call wrtout(my_unt,msg,my_mode)
737 0 : write(msg,fmt)ib,(aimag(data_p(tab(ikey))%arr_p(ib,ib,ikibz,3))*Ha_eV, ikey=1,my_nkeys)
738 0 : call wrtout(my_unt,msg,my_mode)
739 : end if
740 : end if
741 : end do
742 :
743 : else
744 : ! Print full matrix.
745 232 : max_r = b2-b1+1
746 232 : max_c = MIN(b2-b1+1, 9)
747 928 : ABI_MALLOC(mat,(b1:b2,b1:b2))
748 786 : do ikey=1,my_nkeys
749 554 : write(msg,'(3a)')" **** Off-diagonal elements of ",TRIM(my_keys(tab(ikey)))," **** "
750 554 : call wrtout(my_unt,msg,my_mode)
751 1108 : do iab=1,Mels%nspinor**2
752 111072 : mat = data_p(tab(ikey))%arr_p(b1:b2,b1:b2,ikibz,iab) * Ha_eV
753 1662 : call print_arr([my_unt], mat, max_r, max_c)
754 : end do
755 554 : write(msg,'(a)')ch10
756 786 : call wrtout(my_unt,msg,my_mode)
757 : end do
758 232 : ABI_FREE(mat)
759 : end if
760 :
761 : end do !ikibz
762 : end do ! isppol
763 :
764 : 10 continue
765 :
766 277 : ABI_FREE(my_keys)
767 277 : ABI_FREE(data_p)
768 277 : ABI_FREE(tab)
769 :
770 554 : end subroutine melements_print
771 : !!***
772 :
773 : !----------------------------------------------------------------------
774 :
775 : !!****f* m_melemts/melements_zero
776 : !! NAME
777 : !! melements_zero
778 : !!
779 : !! FUNCTION
780 : !! Set matrix elements connecting states with different irreducible representation to zero.
781 : !!
782 : !! INPUTS
783 : !! irrep_tab=Array used to select the entries that have to be set to zero.
784 : !! irrep_tab(ib,ik,is)=gives the index of the irreducible representation associated to state (ib,ik,is).
785 : !! [aname]=The name of the array to be symmetrized, by default
786 : !! all calculated arrays are filled.
787 : !!
788 : !! SIDE EFFECTS
789 : !! Mels= All arrays elements connecting states belonging to different irreps are set to zero.
790 : !!
791 : !! SOURCE
792 :
793 0 : subroutine melements_zero(Mels, irrep_tab, aname)
794 :
795 : !Arguments ------------------------------------
796 : !scalars
797 : class(melements_t),intent(inout) :: Mels
798 : character(len=*),optional,intent(in) :: aname
799 : !arrays
800 : integer,intent(in) :: irrep_tab(Mels%bmin:Mels%bmax,Mels%nkibz,Mels%nsppol)
801 :
802 : !Local variables-------------------------------
803 : integer :: is,ik,ib,jb,iname,irrep_j,irrep_i
804 : integer,pointer :: flag_p
805 : character(len=NAMELEN) :: key
806 : !arrays
807 : complex(dp),contiguous, pointer :: arr_p(:,:,:,:)
808 : ! *************************************************************************
809 :
810 0 : do iname=1,NNAMES
811 0 : key = ANAMES(iname)
812 0 : if (PRESENT(aname)) then
813 0 : if (key /= aname) CYCLE
814 : end if
815 :
816 0 : call my_select_melements(Mels,key,flag_p,arr_p)
817 :
818 0 : if (flag_p>0) then
819 0 : do is=1,Mels%nsppol
820 0 : do ik=1,Mels%nkibz
821 :
822 0 : do jb=Mels%bmin,Mels%bmax
823 0 : irrep_j = irrep_tab(jb,ik,is)
824 0 : do ib=Mels%bmin,Mels%bmax
825 0 : irrep_i = irrep_tab(ib,ik,is)
826 : !
827 : ! Set this matrix element to zero if the irreps are known and they differ.
828 0 : if (irrep_i/=irrep_j .and. ALL((/irrep_i,irrep_j/) /=0) ) then
829 : !write(std_out,*)"setting to zero ",ib,jb,ik,is
830 0 : if (Mels%nspinor==2) then
831 0 : arr_p(ib,jb,ik,is)=czero
832 : else
833 0 : arr_p(ib,jb,ik,:)=czero
834 : end if
835 : end if
836 :
837 : end do !ib
838 : end do !jb
839 :
840 : end do !is
841 : end do !ik
842 : end if
843 :
844 : end do !inames
845 :
846 0 : end subroutine melements_zero
847 : !!***
848 :
849 : !----------------------------------------------------------------------
850 :
851 : !!****f* m_sigma/mels_get_exene_core
852 : !! NAME
853 : !! mels_get_exene_core
854 : !!
855 : !! FUNCTION
856 : !! Compute exchange energy.
857 : !!
858 : !! INPUTS
859 : !! mels<melements_t>=Matrix elements.
860 : !! kmesh<kmesh_t>=BZ sampling.
861 : !! bands<band_t>=Bands with occupation factors
862 : !!
863 : !! SOURCE
864 :
865 : !pure function mels_get_exene_core(mels,kmesh,bands) result(ex_energy)
866 : !
867 : !!Arguments ------------------------------------
868 : !!scalars
869 : ! real(dp) :: ex_energy
870 : ! type(melements_t),intent(in) :: mels
871 : ! type(kmesh_t),intent(in) :: kmesh
872 : ! type(ebands_t),intent(in) :: bands
873 : !
874 : !!Local variables-------------------------------
875 : !!scalars
876 : ! integer :: ik,ib,spin
877 : ! real(dp) :: wtk,occ_bks
878 : !! *************************************************************************
879 : !
880 : ! ex_energy = zero
881 : !
882 : ! do spin=1,mels%nsppol
883 : ! do ik=1,mels%nkibz
884 : ! wtk = kmesh%wt(ik)
885 : ! do ib=mels%bmin,mels%bmax
886 : ! occ_bks = bands%occ(ib,ik,spin)
887 : ! if (mels%nspinor==1) then
888 : ! ex_energy = ex_energy + half * occ_bks * wtk * mels%sxcore(ib,ib,ik,spin)
889 : ! else
890 : ! ex_energy = ex_energy + half * occ_bks wtk *SUM(mels%sxcore(ib,ib,ik,:))
891 : ! end if
892 : ! end do
893 : ! end do
894 : ! end do
895 : !
896 : !end function mels_get_exene_core
897 : !!!***
898 :
899 843 : end module m_melemts
900 : !!***
|
