Branch data 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 : 854 : subroutine melflags_reset(Mflags)
204 : :
205 : : !Arguments ------------------------------------
206 : : class(melflags_t),intent(inout) :: Mflags
207 : : ! *************************************************************************
208 : :
209 : 854 : Mflags%has_kinetic = 0
210 : 854 : Mflags%has_hbare = 0
211 : 854 : Mflags%has_sxcore = 0
212 : 854 : Mflags%has_vhartree = 0
213 : 854 : Mflags%has_vu = 0
214 : 854 : Mflags%has_vxc = 0
215 : 854 : Mflags%has_vxcval = 0
216 : 854 : Mflags%has_vxcval_hybrid = 0
217 : 854 : Mflags%has_lexexch = 0
218 : 854 : Mflags%only_diago = 0
219 : :
220 : 854 : 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 : 286 : 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 : 286 : call Mflags_out%reset()
248 : :
249 : : ! @melflags_t
250 : 286 : Mflags_out%has_kinetic = Mflags_in%has_kinetic
251 : 286 : Mflags_out%has_hbare = Mflags_in%has_hbare
252 : 286 : Mflags_out%has_sxcore = Mflags_in%has_sxcore
253 : 286 : Mflags_out%has_vhartree = Mflags_in%has_vhartree
254 : 286 : Mflags_out%has_vu = Mflags_in%has_vu
255 : 286 : Mflags_out%has_vxc = Mflags_in%has_vxc
256 : 286 : Mflags_out%has_vxcval = Mflags_in%has_vxcval
257 : 286 : Mflags_out%has_vxcval_hybrid = Mflags_in%has_vxcval_hybrid
258 : 286 : Mflags_out%has_lexexch = Mflags_in%has_lexexch
259 : 286 : Mflags_out%only_diago = Mflags_in%only_diago
260 : :
261 : 286 : 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 : 286 : subroutine melements_free(Mels)
276 : :
277 : : !Arguments ------------------------------------
278 : : class(melements_t),intent(inout) :: Mels
279 : : ! *************************************************************************
280 : :
281 : : ! integer arrays
282 [ + - ]: 286 : ABI_SFREE(Mels%bands_idx)
283 [ + - ]: 286 : ABI_SFREE(Mels%iscalc)
284 : :
285 : : ! real arrays
286 [ + - ]: 286 : ABI_SFREE(Mels%kibz)
287 : :
288 : : ! complex arrays
289 [ + + ]: 286 : ABI_SFREE(Mels%kinetic)
290 [ + + ]: 286 : ABI_SFREE(Mels%hbare)
291 [ + + ]: 286 : ABI_SFREE(Mels%sxcore)
292 [ + - ]: 286 : ABI_SFREE(Mels%vhartree)
293 [ - + ]: 286 : ABI_SFREE(Mels%vlexx)
294 [ - + ]: 286 : ABI_SFREE(Mels%vu)
295 [ + + ]: 286 : ABI_SFREE(Mels%vxc)
296 [ + + ]: 286 : ABI_SFREE(Mels%vxcval)
297 [ + + ]: 286 : ABI_SFREE(Mels%vxcval_hybrid)
298 : :
299 : : ! Reset all has_* flags.
300 : 286 : call Mels%flags%reset()
301 : :
302 : 286 : 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 : 7686 : 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 : 16226 : SELECT CASE (tolower(aname))
336 : : CASE ("kinetic")
337 : 854 : flag_p => Mels%flags%has_kinetic
338 : 854 : arr_p => Mels%kinetic
339 : : CASE ("hbare")
340 : 854 : flag_p => Mels%flags%has_hbare
341 : 854 : arr_p => Mels%hbare
342 : : CASE ("sxcore")
343 : 854 : flag_p => Mels%flags%has_sxcore
344 : 854 : arr_p => Mels%sxcore
345 : : CASE ("vhartree")
346 : 854 : flag_p => Mels%flags%has_vhartree
347 : 854 : arr_p => Mels%vhartree
348 : : CASE ("vlexx")
349 : 854 : flag_p => Mels%flags%has_lexexch
350 : 854 : arr_p => Mels%vlexx
351 : : CASE ("vu")
352 : 854 : flag_p => Mels%flags%has_vu
353 : 854 : arr_p => Mels%vu
354 : : CASE ("vxc")
355 : 854 : flag_p => Mels%flags%has_vxc
356 : 854 : arr_p => Mels%vxc
357 : : CASE ("vxcval")
358 : 854 : flag_p => Mels%flags%has_vxcval
359 : 854 : arr_p => Mels%vxcval
360 : : CASE ("vxcval_hybrid")
361 : 854 : flag_p => Mels%flags%has_vxcval_hybrid
362 : 854 : arr_p => Mels%vxcval_hybrid
363 : : CASE DEFAULT
364 [ + + + + : 7686 : ABI_ERROR(sjoin("Wrong aname: ", aname))
+ + + + +
- ]
365 : : END SELECT
366 : :
367 : 7686 : 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 [ + - ]: 286 : 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 : 286 : call Mflags_in%copy(Mels%flags)
408 : :
409 : : ! Copy dimensions.
410 : 286 : Mels%nkibz = nkibz
411 : 286 : Mels%nsppol = nsppol
412 : 286 : Mels%nspinor = nspinor
413 : 286 : Mels%nspden = nspden
414 : :
415 [ - + + - : 1144 : ABI_MALLOC(Mels%bands_idx, (2, nkibz, nsppol))
- + + - +
- + - - +
- + - + ]
416 [ - + + + : 6036 : Mels%bands_idx = bands_idx
+ + + + ]
417 : :
418 [ - + + - : 1144 : ABI_MALLOC(Mels%iscalc,(nkibz, nsppol))
- + + - +
- + - - +
- + - + ]
419 [ + + + + ]: 2302 : Mels%iscalc = 0
420 : :
421 : 286 : bmin = HUGE(1); bmax =-HUGE(1)
422 [ + + ]: 578 : do isppol=1,Mels%nsppol
423 [ + + ]: 2302 : do ikibz=1,Mels%nkibz
424 [ + + + + : 3314 : if (ANY(Mels%bands_idx(:,ikibz,isppol)/=0)) then
+ + ]
425 : 1075 : b1 = Mels%bands_idx(1,ikibz,isppol)
426 : 1075 : b2 = Mels%bands_idx(2,ikibz,isppol)
427 : 1075 : Mels%iscalc(ikibz, isppol)=1
428 : 1075 : bmin = MIN(bmin,b1)
429 : 1075 : bmax = MAX(bmax,b2)
430 [ - + ]: 1075 : ABI_CHECK(b2 >= b1 .and. b1 > 0, "Wrong b1, b2")
431 : : end if
432 : : end do
433 : : end do
434 : :
435 [ - + ]: 286 : if (bmin==HUGE(1).or.bmax==-HUGE(1)) then
436 : 0 : ABI_BUG("Wrong bands_idx")
437 : : end if
438 : :
439 : 286 : Mels%bmin = bmin
440 : 286 : Mels%bmax = bmax
441 : :
442 : 286 : b1 = Mels%bmin; b2 = Mels%bmax
443 : :
444 : : ! real arrays
445 [ - + + - : 858 : ABI_MALLOC(Mels%kibz, (3,nkibz))
+ - + - -
+ - + -
+ ]
446 [ - + + + : 7420 : Mels%kibz = kibz
+ + ]
447 : :
448 : : ! complex arrays
449 [ + + ]: 286 : 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 [ + + ]: 286 : 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 [ + + ]: 286 : if (Mels%flags%has_sxcore == 1) then
456 [ - + + - : 1160 : ABI_CALLOC(Mels%sxcore, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
457 : : end if
458 [ + - ]: 286 : if (Mels%flags%has_vhartree == 1) then
459 [ - + + - : 161034 : ABI_CALLOC(Mels%vhartree, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
460 : : end if
461 [ - + ]: 286 : 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 [ - + ]: 286 : 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 [ + + ]: 286 : if (Mels%flags%has_vxc == 1) then
468 [ - + + - : 108127 : ABI_CALLOC(Mels%vxc, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
469 : : end if
470 [ + + ]: 286 : if (Mels%flags%has_vxcval == 1) then
471 [ - + + - : 108127 : ABI_CALLOC(Mels%vxcval, (b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + + + +
+ + + +
+ ]
472 : : end if
473 [ + + ]: 286 : 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 : 286 : 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 : 286 : 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 : 286 : 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 [ + + ]: 2860 : do iname=1,NNAMES
521 : 2574 : key = ANAMES(iname)
522 [ - + ]: 2574 : if (PRESENT(aname)) then
523 [ # # ]: 0 : if (key /= aname) CYCLE
524 : : end if
525 : :
526 : 2574 : call my_select_melements(Mels, key, flag_p, arr_p)
527 : :
528 [ + + ]: 2860 : if (flag_p>0) then
529 [ + + ]: 5294 : do ik=1,Mels%nkibz
530 [ + + ]: 9871 : do is=1,Mels%nsppol
531 : :
532 [ + + ]: 43729 : do jb=Mels%bmin,Mels%bmax
533 [ + + ]: 215961 : do ib=Mels%bmin,jb ! Upper triangle
534 : :
535 [ + + ]: 211384 : if (ib/=jb) then
536 : 142170 : arr_p(jb,ib,ik,is)=CONJG(arr_p(ib,jb,ik,is))
537 [ + + ]: 142170 : 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 : 34607 : arr_p(jb,ib,ik,is)=half*(arr_p(jb,ib,ik,is)+CONJG(arr_p(jb,ib,ik,is)))
545 [ + + ]: 34607 : 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 : 286 : 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 : 286 : 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 : 286 : complex(dp),contiguous, pointer :: arr_p(:,:,:,:)
596 : : ! *************************************************************************
597 : :
598 [ + + ]: 2860 : do iname=1,NNAMES
599 : 2574 : key = ANAMES(iname)
600 [ - + ]: 2574 : if (PRESENT(aname)) then
601 [ # # ]: 0 : if (key /= aname) CYCLE
602 : : end if
603 : :
604 : 2574 : call my_select_melements(Mels, key, flag_p, arr_p)
605 : :
606 [ + + ]: 2860 : if (flag_p == 1) then
607 : 749 : call xmpi_sum(arr_p, comm, ierr)
608 : 749 : flag_p = 2 ! Tag this array as calculated
609 : : end if
610 : : end do
611 : :
612 : 286 : 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 : 282 : integer,allocatable :: tab(:)
655 : 282 : character(len=NAMELEN),allocatable :: my_keys(:)
656 : 282 : complex(dp),allocatable :: mat(:,:)
657 : :
658 : : type rarr_dpc4
659 : : complex(dp),contiguous, pointer :: arr_p(:,:,:,:)
660 : : end type rarr_dpc4
661 : 282 : type(rarr_dpc4),allocatable :: data_p(:)
662 : : ! *************************************************************************
663 : :
664 : : !@melements_t
665 [ - + ]: 282 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
666 [ + - ]: 282 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
667 [ - + ]: 282 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
668 : :
669 [ + + ]: 282 : if (Mels%nspinor == 2) ABI_WARNING("nspinor=2 not coded")
670 : :
671 [ - + ]: 282 : 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 : 282 : my_nkeys = NNAMES
677 [ - + ]: 282 : ABI_MALLOC(my_keys, (NNAMES))
678 [ - + + + ]: 3102 : my_keys = ANAMES
679 : : end if
680 : :
681 [ - + + - : 846 : ABI_MALLOC(data_p, (my_nkeys))
+ - + - -
+ - + ]
682 [ - + + - : 846 : ABI_MALLOC(tab, (my_nkeys))
+ - + - -
+ - + ]
683 [ + + ]: 2820 : tab = 0
684 : :
685 : 282 : my_nkeys=0; str = " ib"; ii=4
686 [ + + ]: 2820 : do ikey=1,size(my_keys)
687 : 2538 : key = my_keys(ikey)
688 : 2538 : call my_select_melements(Mels,key,flag_p,data_p(ikey)%arr_p)
689 [ + + ]: 2820 : if (flag_p == 2) then
690 : 729 : my_nkeys = my_nkeys+1
691 : 729 : tab(my_nkeys) = ikey
692 [ + - + - ]: 729 : str(ii+1:)=" "//TRIM(tolower(key))
693 : 729 : ii = ii+MAX(1+LEN_TRIM(key),10)
694 [ - + ]: 729 : ABI_CHECK(ii <490, "I'm gonna SIGFAULT!")
695 : : end if
696 : : end do
697 : :
698 : 282 : write(msg,'(2a)')ch10,' === Matrix Elements stored in Mels% [eV] === '
699 [ + - + - ]: 282 : if (PRESENT(header)) write(msg,'(4a)')ch10,' === '//TRIM(ADJUSTL(header))//' [eV] === '
700 : 282 : call wrtout(my_unt,msg,my_mode)
701 [ + + ]: 282 : 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 [ + - ]: 282 : if (my_nkeys==0) GOTO 10
706 : 282 : write(fmt,'(a,i4,a)')'(1x,i3,',my_nkeys,'(1x,f9.5))' ! width of 10 chars
707 : :
708 [ + + ]: 570 : do isppol=1,Mels%nsppol
709 [ + + ]: 2270 : do ikibz=1,Mels%nkibz
710 [ + + ]: 1700 : if (Mels%iscalc(ikibz,isppol)/=1) CYCLE
711 : :
712 : 1051 : write(msg,'(a,3es16.8,a,i2,a)')" kpt= (",Mels%kibz(:,ikibz),") spin=",isppol,":"
713 : 1051 : call wrtout(my_unt,msg,my_mode)
714 : :
715 : 1051 : b1 = Mels%bands_idx(1,ikibz,isppol)
716 : 1051 : b2 = Mels%bands_idx(2,ikibz,isppol)
717 : :
718 [ + + + + ]: 1339 : if (Mels%flags%only_diago==1 .or. my_prtvol==0) then
719 : : ! Print only the diagonal.
720 : 819 : write(msg,'(a)')str
721 : 819 : call wrtout(my_unt,msg,my_mode)
722 [ + + ]: 7214 : do ib=b1,b2
723 [ + + ]: 7214 : if (Mels%nspinor == 1) then
724 [ + - + + ]: 20948 : write(msg,fmt)ib,(REAL(data_p(tab(ikey))%arr_p(ib,ib,ikibz,1))*Ha_eV, ikey=1,my_nkeys)
725 : 6251 : 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 [ - + ]: 282 : ABI_FREE(my_keys)
767 [ - + ]: 282 : ABI_FREE(data_p)
768 [ - + ]: 282 : ABI_FREE(tab)
769 : :
770 [ - + - + ]: 564 : 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 [ - + - - : 858 : end module m_melemts
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + - + -
+ - + - +
- + ]
900 : : !!***
|
