Branch data Line data Source code
1 : : !!****m* ABINIT/m_matlu
2 : : !! NAME
3 : : !! m_matlu
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! INPUTS
14 : : !!
15 : : !! OUTPUT
16 : : !!
17 : : !! NOTES
18 : : !! subroutines in this module must never call
19 : : !! a subroutine of m_oper, m_green, m_self
20 : : !! in order to avoid circular dependancies
21 : : !!
22 : : !! SOURCE
23 : :
24 : : #if defined HAVE_CONFIG_H
25 : : #include "config.h"
26 : : #endif
27 : :
28 : :
29 : : #include "abi_common.h"
30 : :
31 : : MODULE m_matlu
32 : :
33 : : use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
34 : : use defs_basis
35 : : use m_errors
36 : : use m_abicore
37 : : use m_gputk
38 : : use m_abi_linalg
39 : :
40 : : #ifdef HAVE_GPU
41 : : use m_gpu_toolbox
42 : : #endif
43 : :
44 : : use m_abi_linalg, only : abi_xgemm
45 : : use m_fstrings, only : int2char4
46 : : use m_hide_lapack, only : xginv
47 : : use m_io_tools, only : flush_unit
48 : : use m_matrix, only : blockdiago_fordsyev,blockdiago_forzheev
49 : : use m_paw_dmft, only : paw_dmft_type
50 : : use m_xmpi, only : xmpi_bcast,xmpi_sum
51 : :
52 : : implicit none
53 : :
54 : : private
55 : :
56 : : public :: init_matlu
57 : : public :: inverse_matlu
58 : : public :: destroy_matlu
59 : : public :: diff_matlu
60 : : public :: add_matlu
61 : : public :: print_matlu
62 : : public :: sym_matlu
63 : : public :: copy_matlu
64 : : public :: copy_matlu_from_ndat
65 : : public :: copy_matlu_to_ndat
66 : : public :: gather_matlu
67 : : public :: zero_matlu
68 : : public :: trace_matlu
69 : : public :: diag_matlu
70 : : public :: rotate_matlu
71 : : public :: shift_matlu
72 : : public :: checkdiag_matlu
73 : : public :: checkreal_matlu
74 : : public :: prod_matlu
75 : : public :: conjg_matlu
76 : : public :: ln_matlu
77 : : public :: slm2ylm_matlu
78 : : public :: fac_matlu
79 : : public :: printplot_matlu
80 : : public :: identity_matlu
81 : : public :: magmomforb_matlu
82 : : public :: magmomfspin_matlu
83 : : public :: magmomfzeeman_matlu
84 : : public :: chi_matlu
85 : : public :: trace_prod_matlu
86 : : public :: xmpi_matlu
87 : : public :: symmetrize_matlu
88 : : public :: ylm2jmj_matlu
89 : : public :: magnfield_matlu
90 : : public :: magmomjmj_matlu
91 : : !!***
92 : :
93 : : !!****t* m_matlu/matlu_type
94 : : !! NAME
95 : : !! matlu_type
96 : : !!
97 : : !! FUNCTION
98 : : !! This structured datatype contains a matrix for the correlated subspace
99 : : !!
100 : : !! SOURCE
101 : :
102 : : type, public :: matlu_type ! for each atom
103 : :
104 : : integer :: lpawu
105 : : ! Value of the angular momentum for each correlated electrons
106 : :
107 : : ! integer :: natom
108 : : ! number of atoms (given for each atom, not useful..could be changed)
109 : : !
110 : : ! integer :: mband
111 : : ! ! Number of bands
112 : : !
113 : : ! integer :: mbandc
114 : : ! ! Total number of bands in the Kohn-Sham Basis for PAW+DMFT
115 : : !
116 : : ! integer :: natpawu ! Number of correlated atoms
117 : : !
118 : : ! integer :: nkpt
119 : : ! ! Number of k-point in the IBZ.
120 : : !character(len=12) :: whichmatlu
121 : : ! describe the type of local matrix computed (greenDFT, etc..)
122 : : !
123 : : integer :: gpu_option
124 : : ! Wether ks and matlu are stored on GPU
125 : : !
126 : : integer :: ndat
127 : : ! Number of elements computed in batch
128 : : !
129 : : integer :: nspinor
130 : : ! Number of spinorial components
131 : : !
132 : : integer :: nsppol
133 : : ! Number of polarizations
134 : :
135 : : complex(dp), allocatable :: mat(:,:,:)
136 : : ! Local quantity
137 : :
138 : : end type matlu_type
139 : :
140 : : !----------------------------------------------------------------------
141 : :
142 : :
143 : : CONTAINS !========================================================================================
144 : : !!***
145 : :
146 : : !!****f* m_matlu/init_matlu
147 : : !! NAME
148 : : !! init_matlu
149 : : !!
150 : : !! FUNCTION
151 : : !! Allocate variables used in type matlu_type.
152 : : !!
153 : : !! INPUTS
154 : : !! natom = number of atoms
155 : : !! nspinor = number of spinorial components
156 : : !! nsppol = number of polarisation components
157 : : !! lpawu_natom(natom) = value of lpawu for every atom
158 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
159 : : !!
160 : : !! OUTPUTS
161 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
162 : : !!
163 : : !! SOURCE
164 : :
165 : 361373 : subroutine init_matlu(natom,nspinor,nsppol,lpawu_natom,matlu,gpu_option,ndat)
166 : :
167 : : !Arguments ------------------------------------
168 : : integer, intent(in) :: natom,nspinor,nsppol
169 : : integer, intent(in) :: lpawu_natom(natom)
170 : : integer, intent(in), optional :: gpu_option,ndat
171 : : type(matlu_type), target, intent(inout) :: matlu(natom)
172 : : !Local variables ------------------------------------
173 : : integer :: iatom,lpawu,ndim,l_gpu_option,l_ndat
174 : : !************************************************************************
175 : :
176 [ + + ]: 361373 : l_gpu_option=ABI_GPU_DISABLED; if(present(gpu_option)) l_gpu_option=gpu_option
177 [ - + ]: 361373 : l_ndat=1; if(present(ndat)) l_ndat=ndat
178 : : ! matlu%mband = mband
179 : : ! matlu%dmftbandf = dmftbandf
180 : : ! matlu%dmftbandi = dmftbandi
181 : : ! matlu%nkpt = nkpt
182 : : ! matlu%mbandc = 0
183 [ + + ]: 980224 : do iatom=1,natom
184 : :
185 : 618851 : lpawu = lpawu_natom(iatom)
186 : 618851 : matlu(iatom)%lpawu = lpawu
187 : 618851 : matlu(iatom)%nspinor = nspinor
188 : 618851 : matlu(iatom)%nsppol = nsppol
189 : 618851 : matlu(iatom)%ndat = l_ndat
190 : 618851 : matlu(iatom)%gpu_option = l_gpu_option
191 [ + + ]: 618851 : if (lpawu == -1) cycle
192 : 369523 : ndim = (2*lpawu+1) * nspinor
193 [ - + + - : 1847615 : ABI_MALLOC(matlu(iatom)%mat,(ndim,ndim,nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
194 [ + - ]: 730896 : if(l_gpu_option==ABI_GPU_DISABLED) then
195 [ + + + + : 19895802 : matlu(iatom)%mat(:,:,:) = czero
+ + ]
196 [ # # ]: 0 : else if(l_gpu_option==ABI_GPU_OPENMP) then
197 [ # # # # : 0 : matlu(iatom)%mat(:,:,:) = czero
# # ]
198 : : #ifdef HAVE_OPENMP_OFFLOAD
199 : : call matlu_ompgpu_map_alloc(matlu(iatom)%mat)
200 : : #endif
201 : 0 : call gpu_set_to_zero_complex(matlu(iatom)%mat, int(nsppol,c_size_t)*ndim*ndim)
202 : : end if
203 : :
204 : : end do ! iatom
205 : :
206 : : #if defined HAVE_OPENMP_OFFLOAD
207 : : contains
208 : : ! mat is located in an array of structs,
209 : : ! which syntax is loosely supported in OpenMP directives
210 : : subroutine matlu_ompgpu_map_alloc(mat)
211 : : complex(dp),intent(inout) :: mat(:,:,:)
212 : : !$OMP TARGET ENTER DATA MAP(alloc:mat)
213 : : end subroutine matlu_ompgpu_map_alloc
214 : : #endif
215 : 361373 : end subroutine init_matlu
216 : : !!***
217 : :
218 : : !!****f* m_matlu/zero_matlu
219 : : !! NAME
220 : : !! zero_matlu
221 : : !!
222 : : !! FUNCTION
223 : : !! Set the elements of matlu to 0.
224 : : !!
225 : : !! INPUTS
226 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
227 : : !! natom = number of atoms
228 : : !! onlynondiag = set all the off-diagonal elements to 0
229 : : !! onlyimag = set the imaginary part to 0
230 : : !!
231 : : !! OUTPUT
232 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
233 : : !! err = maximal off-diagonal/imaginary element that is neglected
234 : : !!
235 : : !! SOURCE
236 : :
237 : 18975 : subroutine zero_matlu(matlu,natom,onlynondiag,onlyimag,err)
238 : :
239 : : !Arguments ------------------------------------
240 : : integer, intent(in) :: natom
241 : : type(matlu_type), intent(inout) :: matlu(natom)
242 : : integer, optional, intent(in) :: onlyimag,onlynondiag
243 : : real(dp), optional, intent(out) :: err
244 : : !Local variables-------------------------------
245 : : integer :: iatom,im,im1,isppol
246 : : integer :: lpawu,ndim,nspinor,nsppol,tndim
247 : : real(dp) :: err_
248 : : !*********************************************************************
249 : :
250 : 18975 : nspinor = matlu(1)%nspinor
251 : 18975 : nsppol = matlu(1)%nsppol
252 : :
253 [ - + ]: 18975 : if (present(err)) err = zero
254 : :
255 [ + + ]: 42949 : do iatom=1,natom
256 : 23974 : lpawu = matlu(iatom)%lpawu
257 [ + + ]: 23974 : if (lpawu == -1) cycle
258 [ + + ]: 38142 : if (present(onlynondiag)) then
259 : 17929 : ndim = 2*lpawu + 1
260 : 17929 : tndim = nspinor * ndim
261 [ + + ]: 52527 : do isppol=1,nsppol
262 [ + + ]: 225685 : do im1=1,tndim
263 [ + + ]: 1084130 : do im=1,tndim
264 [ + + ]: 1049532 : if (im /= im1) then
265 [ - + ]: 703216 : if (present(err)) then
266 : 0 : err_ = abs(matlu(iatom)%mat(im,im1,isppol))
267 [ # # ]: 0 : if (err_ > err) err = err_
268 : : end if
269 : 703216 : matlu(iatom)%mat(im,im1,isppol) = czero
270 : : end if ! im/=im1
271 : : end do ! im
272 : : end do ! im1
273 : : end do ! isppol
274 [ - + ]: 1238 : else if (present(onlyimag)) then
275 [ # # # # : 0 : if (present(err)) err = maxval(abs(aimag(matlu(iatom)%mat(:,:,:))))
# # # # #
# # # # #
# # ]
276 [ # # # # : 0 : matlu(iatom)%mat(:,:,:) = cmplx(dble(matlu(iatom)%mat(:,:,:)),zero,kind=dp)
# # ]
277 : : else
278 [ + + + + : 61355 : matlu(iatom)%mat(:,:,:) = czero
+ + ]
279 : : end if ! onlynondiag
280 : : end do ! iatom
281 : :
282 : 18975 : end subroutine zero_matlu
283 : : !!***
284 : :
285 : : !!****f* m_matlu/destroy_matlu
286 : : !! NAME
287 : : !! destroy_matlu
288 : : !!
289 : : !! FUNCTION
290 : : !! Deallocate matlu
291 : : !!
292 : : !! INPUTS
293 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
294 : : !! natom = number of atoms
295 : : !!
296 : : !! OUTPUT
297 : : !!
298 : : !! SOURCE
299 : :
300 : 361373 : subroutine destroy_matlu(matlu,natom)
301 : :
302 : : !Arguments ------------------------------------
303 : : integer, intent(in) :: natom
304 : : type(matlu_type),target, intent(inout) :: matlu(natom)
305 : : !Local variables-------------------------------
306 : : integer :: iatom
307 : : ! *********************************************************************
308 : :
309 [ + + ]: 980224 : do iatom=1,natom
310 : : if(matlu(iatom)%gpu_option==ABI_GPU_OPENMP) then
311 : : #ifdef HAVE_OPENMP_OFFLOAD
312 : : call matlu_ompgpu_map_delete(matlu(iatom)%mat)
313 : : #endif
314 : : end if
315 [ + + ]: 980224 : ABI_SFREE(matlu(iatom)%mat)
316 : : end do ! iatom
317 : :
318 : : #if defined HAVE_OPENMP_OFFLOAD
319 : : contains
320 : : ! mat is located in an array of structs,
321 : : ! which syntax is loosely supported in OpenMP directives
322 : : subroutine matlu_ompgpu_map_delete(mat)
323 : : complex(dp),intent(inout) :: mat(:,:,:)
324 : : !$OMP TARGET EXIT DATA MAP(delete:mat)
325 : : end subroutine matlu_ompgpu_map_delete
326 : : #endif
327 : 361373 : end subroutine destroy_matlu
328 : : !!***
329 : :
330 : : !!****f* m_matlu/copy_matlu
331 : : !! NAME
332 : : !! copy_matlu
333 : : !!
334 : : !! FUNCTION
335 : : !! Copy mat1 into mat2
336 : : !!
337 : : !! INPUTS
338 : : !! mat1 <type(matlu_type)>= density matrix nmat1 in the local orbital basis and related variables
339 : : !! natom = number of atoms
340 : : !! opt_diag = if present, only copy the diagonal elements (the off-diagonal elements are not set to 0)
341 : : !! opt_non_diag = if present, only copy the off-diagonal elements
342 : : !! opt_re = if present, only copy the real part
343 : : !!
344 : : !! OUTPUT
345 : : !! mat2 <type(matlu_type)>= density matrix nmat2 in the local orbital basis and related variables
346 : : !!
347 : : !! SOURCE
348 : :
349 : 170189 : subroutine copy_matlu(mat1,mat2,natom,opt_diag,opt_non_diag,opt_re)
350 : :
351 : : !Arguments ------------------------------------
352 : : integer, intent(in) :: natom
353 : : type(matlu_type), intent(in) :: mat1(natom)
354 : : type(matlu_type), intent(inout) :: mat2(natom) !vz_i
355 : : integer, optional, intent(in) :: opt_diag,opt_non_diag,opt_re
356 : : !Local variables-------------------------------
357 : : integer :: iatom,im,im1,isppol,lpawu,ndim,nspinor,nsppol
358 : : ! *********************************************************************
359 : :
360 : 170189 : nspinor = mat1(1)%nspinor
361 : 170189 : nsppol = mat1(1)%nsppol
362 : :
363 [ + + ]: 442556 : do iatom=1,natom
364 : :
365 : 272367 : lpawu = mat1(iatom)%lpawu
366 [ + + ]: 272367 : if (lpawu == -1) cycle
367 : 174413 : ndim = (2*lpawu+1) * nspinor
368 : :
369 : : !if both matrix are on GPU and no opt is provided, perform copy on GPU
370 : : !Other cases can be handled through OpenMP kernels on GPU but no use case exists.
371 : : if(mat1(iatom)%gpu_option==ABI_GPU_OPENMP .and. mat2(iatom)%gpu_option==ABI_GPU_OPENMP &
372 [ - + - - : 174413 : & .and. .not. (present(opt_diag) .or. present(opt_non_diag) .or. present(opt_re))) then
- - - - ]
373 : : #ifdef HAVE_OPENMP_OFFLOAD
374 : : call gpu_copy_complex(mat2(iatom)%mat, mat1(iatom)%mat, int(nsppol, c_size_t)*ndim*ndim)
375 : : #endif
376 : : cycle
377 : : end if
378 : :
379 : :
380 : : #ifdef HAVE_OPENMP_OFFLOAD
381 : : !$OMP TARGET UPDATE FROM(mat1(iatom)%mat) IF(mat1(iatom)%gpu_option==ABI_GPU_OPENMP)
382 : : !$OMP TARGET UPDATE FROM(mat2(iatom)%mat) IF(mat2(iatom)%gpu_option==ABI_GPU_OPENMP)
383 : : #endif
384 [ - + ]: 344602 : if (present(opt_diag)) then
385 [ # # ]: 0 : do isppol=1,nsppol
386 [ # # ]: 0 : do im=1,ndim
387 : 0 : mat2(iatom)%mat(im,im,isppol) = mat1(iatom)%mat(im,im,isppol)
388 : : end do ! im
389 : : end do ! isppol
390 [ - + ]: 174413 : else if (present(opt_non_diag)) then
391 [ # # ]: 0 : do isppol=1,nsppol
392 [ # # ]: 0 : do im1=1,ndim
393 [ # # ]: 0 : do im=1,ndim
394 [ # # ]: 0 : if (im /= im1) mat2(iatom)%mat(im,im1,isppol) = mat1(iatom)%mat(im,im1,isppol)
395 : : end do ! im
396 : : end do ! im1
397 : : end do ! isppol
398 [ + + ]: 174413 : else if (present(opt_re)) then
399 [ + + + + : 6438 : mat2(iatom)%mat(:,:,:) = cmplx(dble(mat1(iatom)%mat(:,:,:)),zero,kind=dp)
+ + ]
400 : : else
401 [ + + + + : 10115265 : mat2(iatom)%mat(:,:,:) = mat1(iatom)%mat(:,:,:)
+ + ]
402 : : end if ! opt
403 : : #ifdef HAVE_OPENMP_OFFLOAD
404 : : !$OMP TARGET UPDATE TO(mat2(iatom)%mat) IF(mat2(iatom)%gpu_option==ABI_GPU_OPENMP)
405 : : #endif
406 : :
407 : : end do ! iatom
408 : :
409 : : ! do iatom=1,natom
410 : : ! lpawu=nmat1(iatom)%lpawu
411 : : ! if(lpawu.ne.-1) then
412 : : ! nmat2(iatom)%mat=nmat1(iatom)%mat
413 : : ! endif
414 : : ! enddo
415 : :
416 : 170189 : end subroutine copy_matlu
417 : : !!***
418 : :
419 : : !!****f* m_matlu/copy_matlu_from_ndat
420 : : !! NAME
421 : : !! copy_matlu_from_ndat
422 : : !!
423 : : !! FUNCTION
424 : : !! Copy matlu1 into matlu2
425 : : !!
426 : : !! INPUTS
427 : : !! maltu1 <type(matlu_type)>= density matrix matlu1 in the local orbital basis and related variables
428 : : !! natom = number of atoms
429 : : !!
430 : : !! OUTPUT
431 : : !! maltu2 <type(matlu_type)>= density matrix matlu2 in the local orbital basis and related variables
432 : : !!
433 : : !! SOURCE
434 : :
435 : 0 : subroutine copy_matlu_from_ndat(mat1,mat2,natom,ndat,idat,opt_diag,opt_non_diag,opt_re)
436 : :
437 : : !Arguments ------------------------------------
438 : : !type
439 : : integer, intent(in) :: natom,ndat,idat
440 : : type(matlu_type),intent(in) :: mat1(natom)
441 : : type(matlu_type),intent(inout) :: mat2(natom) !vz_i
442 : : integer, optional, intent(in) :: opt_diag,opt_non_diag,opt_re
443 : :
444 : : !Local variables-------------------------------
445 : : integer :: iatom,isppol,im1,im,ndim,nsppol,nspinor,lpawu
446 : : ! *********************************************************************
447 : :
448 [ # # ]: 0 : ABI_CHECK(mat1(1)%nsppol==mat2(1)%nsppol*ndat, "bad ndat value")
449 : 0 : nspinor = mat1(1)%nspinor
450 : 0 : nsppol = mat2(1)%nsppol
451 [ # # ]: 0 : do iatom=1,natom
452 : 0 : lpawu = mat1(iatom)%lpawu
453 [ # # ]: 0 : if (lpawu == -1) cycle
454 : 0 : ndim = (2*lpawu+1) * nspinor
455 : :
456 [ # # ]: 0 : if (present(opt_diag)) then
457 [ # # ]: 0 : do isppol=1,nsppol
458 [ # # ]: 0 : do im=1,ndim
459 : 0 : mat2(iatom)%mat(im,im,isppol) = mat1(iatom)%mat(im,im,(isppol-1)*ndat+idat)
460 : : end do ! im
461 : : end do ! isppol
462 [ # # ]: 0 : else if (present(opt_non_diag)) then
463 [ # # ]: 0 : do isppol=1,nsppol
464 [ # # ]: 0 : do im1=1,ndim
465 [ # # ]: 0 : do im=1,ndim
466 [ # # ]: 0 : if (im /= im1) mat2(iatom)%mat(im,im1,isppol) = mat1(iatom)%mat(im,im1,(isppol-1)*ndat+idat)
467 : : end do ! im
468 : : end do ! im1
469 : : end do ! isppol
470 : :
471 [ # # ]: 0 : else if (present(opt_re)) then
472 [ # # ]: 0 : do isppol=1,nsppol
473 [ # # # # ]: 0 : mat2(iatom)%mat(:,:,isppol) = cmplx(dble(mat1(iatom)%mat(:,:,(isppol-1)*ndat+idat)),zero,kind=dp)
474 : : end do ! isppol
475 : : else
476 [ # # ]: 0 : do isppol=1,nsppol
477 [ # # # # ]: 0 : mat2(iatom)%mat(:,:,isppol) = mat1(iatom)%mat(:,:,(isppol-1)*ndat+idat)
478 : : end do ! isppol
479 : : end if ! opt
480 : : enddo ! iatom
481 : :
482 : 0 : end subroutine copy_matlu_from_ndat
483 : : !!***
484 : :
485 : : !!****f* m_matlu/copy_matlu_to_ndat
486 : : !! NAME
487 : : !! copy_matlu_to_ndat
488 : : !!
489 : : !! FUNCTION
490 : : !! Copy matlu1 into matlu2
491 : : !!
492 : : !! INPUTS
493 : : !! maltu1 <type(matlu_type)>= density matrix matlu1 in the local orbital basis and related variables
494 : : !! natom = number of atoms
495 : : !!
496 : : !! OUTPUT
497 : : !! maltu2 <type(matlu_type)>= density matrix matlu2 in the local orbital basis and related variables
498 : : !!
499 : : !! SOURCE
500 : :
501 : 0 : subroutine copy_matlu_to_ndat(mat1,mat2,natom,ndat,idat,opt_diag,opt_non_diag,opt_re)
502 : :
503 : : !Arguments ------------------------------------
504 : : !type
505 : : integer, intent(in) :: natom,ndat,idat
506 : : type(matlu_type),intent(in) :: mat1(natom)
507 : : type(matlu_type),intent(inout) :: mat2(natom) !vz_i
508 : : integer, optional, intent(in) :: opt_diag,opt_non_diag,opt_re
509 : :
510 : : !Local variables-------------------------------
511 : : integer :: iatom,isppol,im1,im,ndim,nsppol,nspinor,lpawu
512 : : ! *********************************************************************
513 : :
514 : :
515 [ # # ]: 0 : ABI_CHECK(mat1(1)%nsppol*ndat==mat2(1)%nsppol, "bad ndat value")
516 : 0 : nspinor = mat1(1)%nspinor
517 : 0 : nsppol = mat1(1)%nsppol
518 [ # # ]: 0 : do iatom=1,natom
519 : 0 : lpawu = mat1(iatom)%lpawu
520 [ # # ]: 0 : if (lpawu == -1) cycle
521 : 0 : ndim = (2*lpawu+1) * nspinor
522 : :
523 [ # # ]: 0 : if (present(opt_diag)) then
524 [ # # ]: 0 : do isppol=1,nsppol
525 [ # # ]: 0 : do im=1,ndim
526 : 0 : mat2(iatom)%mat(im,im,(isppol-1)*ndat+idat) = mat1(iatom)%mat(im,im,isppol)
527 : : end do ! im
528 : : end do ! isppol
529 [ # # ]: 0 : else if (present(opt_non_diag)) then
530 [ # # ]: 0 : do isppol=1,nsppol
531 [ # # ]: 0 : do im1=1,ndim
532 [ # # ]: 0 : do im=1,ndim
533 [ # # ]: 0 : if (im /= im1) mat2(iatom)%mat(im,im1,(isppol-1)*ndat+idat) = mat1(iatom)%mat(im,im1,isppol)
534 : : end do ! im
535 : : end do ! im1
536 : : end do ! isppol
537 : :
538 [ # # ]: 0 : else if (present(opt_re)) then
539 [ # # ]: 0 : do isppol=1,nsppol
540 [ # # # # ]: 0 : mat2(iatom)%mat(:,:,(isppol-1)*ndat+idat) = cmplx(dble(mat1(iatom)%mat(:,:,isppol)),zero,kind=dp)
541 : : end do ! isppol
542 : : else
543 [ # # ]: 0 : do isppol=1,nsppol
544 [ # # # # ]: 0 : mat2(iatom)%mat(:,:,(isppol-1)*ndat+idat) = mat1(iatom)%mat(:,:,isppol)
545 : : end do ! isppol
546 : : end if ! opt
547 : : enddo ! iatom
548 : :
549 : 0 : end subroutine copy_matlu_to_ndat
550 : : !!***
551 : :
552 : : !!****f* m_matlu/print_matlu
553 : : !! NAME
554 : : !! print_matlu
555 : : !!
556 : : !! FUNCTION
557 : : !! Print matlu
558 : : !!
559 : : !! INPUTS
560 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
561 : : !! natom= number of atoms
562 : : !! prtopt= /=0 print matlu
563 : : !! >=5 print matlu in n,mx,my,mz representation
564 : : !! opt_diag= 0 print non diagonal matrix (real or complex according to nspinor)
565 : : !! -1 print non diagonal complex matrix
566 : : !! >=1 print diagonal matrix (real or complex according to nspinor)
567 : : !! opt_ab_out= 0 print matrix on std_out
568 : : !! /=0 print matrix on ab_out
569 : : !! opt_exp= write in exponent format if present
570 : : !! argout= output unit
571 : : !! compl= 1 print complex matrix
572 : : !! /=1 only print complex matrix if nspinor=2
573 : : !!
574 : : !! OUTPUT
575 : : !!
576 : : !! SOURCE
577 : :
578 : 4835 : subroutine print_matlu(matlu,natom,prtopt,opt_diag,opt_ab_out,opt_exp,argout,compl)
579 : :
580 : : !Arguments ------------------------------------
581 : : integer, intent(in):: natom,prtopt
582 : : type(matlu_type), intent(in) :: matlu(natom)
583 : : integer, optional, intent(in) :: opt_diag,opt_ab_out,opt_exp,argout,compl
584 : : !Local variables-------------------------------
585 : : integer :: arg_out,iatom,im,im1,ispinor,ispinor1,isppol,lpawu
586 : : integer :: ndim,nspinor,nsppol,optab_out,optdiag
587 : : logical :: testcmplx,testcmplx_
588 : 4835 : complex(dp), allocatable :: mat_nmrep(:,:)
589 : : character(len=500) :: message
590 : : character(len=4) :: mode_paral,tag_at
591 : : character(len=9), parameter :: dspinm(2,2) = RESHAPE((/"n ","mx ","my ","mz "/),(/2,2/))
592 : : ! *********************************************************************
593 : :
594 : 4835 : arg_out = ab_out
595 : 4835 : mode_paral = 'COLL'
596 : 4835 : optab_out = 0
597 : 4835 : optdiag = 0
598 : :
599 [ + + ]: 4835 : if (present(opt_diag)) optdiag = opt_diag
600 [ + + ]: 4835 : if (present(opt_ab_out)) optab_out = opt_ab_out
601 [ - + ]: 4835 : if (optab_out == 0) arg_out = std_out
602 : :
603 [ - + ]: 4835 : if (present(argout)) then
604 : 0 : arg_out = argout
605 : 0 : mode_paral = 'PERS'
606 : : end if
607 : :
608 : 4835 : nspinor = matlu(1)%nspinor
609 : 4835 : nsppol = matlu(1)%nsppol
610 : 4835 : testcmplx_ = (nspinor == 2)
611 [ + + + + : 4835 : if (present(compl)) testcmplx_ = (nspinor == 2) .or. (compl == 1)
+ - ]
612 : :
613 [ + + ]: 19929 : do iatom=1,natom
614 : :
615 : 15094 : lpawu = matlu(iatom)%lpawu
616 [ + + ]: 15094 : if (lpawu == -1) cycle
617 : 5495 : ndim = 2*lpawu + 1
618 : :
619 : 5495 : write(tag_at,'(i4)') iatom
620 : 5495 : write(message,'(3a)') ch10,' -------> For Correlated Atom ',adjustl(tag_at)
621 : 5495 : call wrtout(arg_out,message,mode_paral)
622 : :
623 : 5495 : testcmplx = testcmplx_
624 [ + + + + : 255524 : if (maxval(abs(aimag(matlu(iatom)%mat(:,:,:)))) > tol5) testcmplx = .true.
+ + + + +
- + - +
+ ]
625 : :
626 : : !do isppol=1,nsppol
627 : : ! if (present(opt_ab_out) .and. nsppol == 2) then
628 : : ! noccspin = zero
629 : : ! do im=1,ndim
630 : : ! noccspin = noccspin + REAL(matlu(iatom)%mat(im,im,isppol))
631 : : ! end do
632 : : !write(message,fmt='(7x,a,i3,a,f10.5)') ". Occ. for lpawu and for spin",isppol," =",noccspin
633 : : !call wrtout(arg_out, message,mode_paral)
634 : : ! end if
635 : : !end do ! isppol
636 : :
637 [ + + ]: 18405 : do isppol=1,nsppol
638 [ + + ]: 8075 : if (nspinor == 1) then
639 : 7996 : write(message,'(a,10x,a,1x,i1)') ch10,'-- polarization spin component',isppol
640 : 7996 : call wrtout(arg_out,message,mode_paral)
641 : : end if ! nspinor=1
642 [ + + ]: 16229 : do ispinor=1,nspinor
643 [ + + ]: 24541 : do ispinor1=1,nspinor
644 [ + + ]: 8312 : if (nspinor == 2) then
645 : 316 : write(message,'(a,10x,a,i1,1x,i1)') ch10,'-- spin components ',ispinor,ispinor1
646 : 316 : call wrtout(arg_out,message,mode_paral)
647 : : end if
648 [ + - ]: 16466 : if (optdiag <= 0) then
649 [ + + ]: 47940 : do im1=1,ndim
650 [ + + ]: 39628 : if (optdiag == 0) then
651 [ + + + + ]: 34926 : if ((.not. testcmplx) .and. (abs(prtopt) > 0)) then
652 [ + + ]: 21653 : if (present(opt_exp)) then
653 [ + - + + ]: 6328 : write(message,'(5x,20e24.14)') (dble(matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
654 : : ! call wrtout(arg_out, message,mode_paral)
655 : : ! write(message,'(5x,20e20.14)') (REAL(sqrt(matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1))),m=1,2*lpawu+1)
656 : : ! call wrtout(arg_out, message,mode_paral)
657 : : ! write(message,'(5x,20e20.14)') (REAL(1.d0/sqrt(matlu(iatom)%mat(m,m,isppol,ispinor,ispinor1))),m=1,2*lpawu+1)
658 : : else
659 [ + - + + ]: 122804 : write(message,'(5x,20f10.5)') (dble(matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
660 : : end if ! opt_exp
661 [ + - ]: 13201 : else if (testcmplx .and. (abs(prtopt) > 0)) then
662 [ + + ]: 13201 : if (present(opt_exp)) then
663 [ - + ]: 2170 : if (opt_exp == 2) then
664 [ # # # # ]: 0 : write(message,'(5x,14(2e18.10,1x))') ((matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
665 : : else
666 [ + - + + ]: 13464 : write(message,'(5x,14(2e14.4,2x))') ((matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
667 : : end if ! opt_exp=2
668 : : else
669 [ + - + + ]: 72100 : write(message,'(5x,14(2f9.5,2x))') ((matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
670 : : end if ! opt_exp
671 : : !& write(message,'(5x,14(2f15.11,2x))')((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
672 : : end if ! testcmplx
673 [ + - ]: 4702 : else if (optdiag == -1) then
674 [ + - + + ]: 28076 : write(message,'(5x,14(2f10.5,2x))') ((matlu(iatom)%mat(im1+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
675 : : end if ! optdiag
676 : 47940 : call wrtout(arg_out,message,mode_paral)
677 : : end do ! im1
678 : : else if (optdiag >= 1) then
679 [ # # # # ]: 0 : if ((.not. testcmplx) .and. (abs(prtopt) > 0)) write(message,'(5x,20f10.5)') &
680 [ # # # # ]: 0 : & (dble(matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
681 [ # # ]: 0 : if (testcmplx .and. (abs(prtopt) > 0)) write(message,'(5x,14(2f9.5,2x))') &
682 [ # # # # ]: 0 : & ((matlu(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor1-1)*ndim,isppol)),im=1,ndim)
683 : : ! write(std_out,'(5x,14(2f9.5,2x))')((matlu(iatom)%mat(m1,m,isppol,ispinor,ispinor1)),m=1,2*lpawu+1)
684 : 0 : call wrtout(arg_out,message,mode_paral)
685 : : end if ! optdiag
686 : : end do ! ispinor1
687 : : end do ! ispinor
688 [ + + - + ]: 23169 : if (nspinor == 2 .and. prtopt >= 5) then
689 [ # # # # : 0 : ABI_MALLOC(mat_nmrep,(2*ndim,2*ndim)) ! Put matlu in n,mx,my,mz representation
# # # # #
# # # # #
# # # # ]
690 [ # # ]: 0 : do im=1,ndim
691 [ # # ]: 0 : do im1=1,ndim
692 : 0 : mat_nmrep(im1,im) = matlu(iatom)%mat(im1,im,isppol) + matlu(iatom)%mat(im1+ndim,im+ndim,isppol) ! n
693 : 0 : mat_nmrep(im1+ndim,im+ndim) = matlu(iatom)%mat(im1,im,isppol) - matlu(iatom)%mat(im1+ndim,im+ndim,isppol) ! mz
694 : 0 : mat_nmrep(im1+ndim,im) = matlu(iatom)%mat(im1,im+ndim,isppol) + matlu(iatom)%mat(im1+ndim,im,isppol) ! mx
695 : 0 : mat_nmrep(im1,im+ndim) = (matlu(iatom)%mat(im1,im+ndim,isppol)-matlu(iatom)%mat(im+ndim,im,isppol)) * j_dpc ! my
696 : : end do ! im1
697 : : end do ! im
698 [ # # ]: 0 : do ispinor=1,nspinor
699 [ # # ]: 0 : do ispinor1=1,nspinor
700 : 0 : write(message,'(a,10x,2a)') ch10,'-- spin components',dspinm(ispinor1,ispinor)
701 : 0 : call wrtout(arg_out,message,mode_paral)
702 [ # # ]: 0 : do im1=1,ndim
703 [ # # # # ]: 0 : write(message,'(5x,14(2f9.5,2x))') ((mat_nmrep(im1+(ispinor1-1)*ndim,im+(ispinor-1)*ndim)),im=1,ndim)
704 : 0 : call wrtout(arg_out,message,mode_paral)
705 : : end do ! im1
706 : : end do ! ispinor1
707 : : end do ! ispinor
708 [ # # ]: 0 : ABI_FREE(mat_nmrep)
709 : : end if ! nspinor=2 and prtopt >=5
710 : : end do ! isppol
711 : : ! if(nsppol==1.and.nspinor==1) then
712 : : ! write(message,'(a,10x,a,i3,a)') ch10,'-- polarization spin component',isppol+1,' is identical'
713 : : ! call wrtout(arg_out, message,mode_paral)
714 : : ! endif
715 : : end do ! iatom
716 : :
717 : 4835 : end subroutine print_matlu
718 : : !!***
719 : :
720 : : !!****f* m_matlu/sym_matlu
721 : : !! NAME
722 : : !! sym_matlu
723 : : !!
724 : : !! FUNCTION
725 : : !! Symmetrize local quantity.
726 : : !!
727 : : !! COPYRIGHT
728 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
729 : : !! This file is distributed under the terms of the
730 : : !! GNU General Public License, see ~abinit/COPYING
731 : : !! or http://www.gnu.org/copyleft/gpl.txt .
732 : : !!
733 : : !! INPUTS
734 : : !! gloc(natom) <type(matlu_type)>= density matrix in the local orbital basis and related variables
735 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
736 : : !!
737 : : !! OUTPUT
738 : : !! gloc(natom) <type(matlu_type)>= density matrix symmetrized in the local orbital basis and related variables
739 : : !!
740 : : !! SIDE EFFECTS
741 : : !!
742 : : !! NOTES
743 : : !!
744 : : !! SOURCE
745 : :
746 : 87632 : subroutine sym_matlu(gloc,paw_dmft)
747 : :
748 : : !Arguments ------------------------------------
749 : : type(paw_dmft_type), target, intent(in) :: paw_dmft
750 : : type(matlu_type), target, intent(inout) :: gloc(paw_dmft%natom)
751 : : !Local variables-------------------------------
752 : : integer :: at_indx,iatom,irot,isppol,lpawu,m1,m2,mu,natom
753 : : integer :: ndim,ndim_max,nspinor,nsppol,nsym,nu,gpu_option
754 : 87632 : complex(dp), target, allocatable :: gloc_tmp(:,:,:),gloc_tmp2(:,:,:)
755 : 87632 : complex(dp), target, allocatable :: gloc_tmp3(:,:,:,:),gloc_tmp4(:,:,:,:)
756 : 87632 : type(matlu_type), allocatable, target :: gloc_nmrep(:),glocsym(:)
757 : 87632 : complex(dp), ABI_CONTIGUOUS pointer :: zarot(:,:,:,:),gloc_mat(:,:,:),glocsym_mat(:,:,:)
758 : 87632 : real(dp), ABI_CONTIGUOUS pointer :: symrec_cart(:,:,:)
759 : : complex(dp) :: ratio
760 : :
761 : 87632 : natom = paw_dmft%natom
762 : 87632 : ndim_max = 2*paw_dmft%maxlpawu + 1
763 : 87632 : nspinor = paw_dmft%nspinor
764 : 87632 : nsppol = gloc(1)%nsppol
765 : 87632 : nsym = paw_dmft%nsym
766 : 87632 : gpu_option = gloc(1)%gpu_option
767 : 87632 : zarot => paw_dmft%zarot
768 : 87632 : ratio = dcmplx(1.0_dp/nsym,0.0_dp)
769 : :
770 : : !zarot => paw_dmft%zarot(:,1:ndim,irot,lpawu+1)
771 [ + - + - : 402075 : ABI_MALLOC(glocsym,(natom))
+ - - + -
+ + + ]
772 : : #ifdef HAVE_OPENMP_OFFLOAD
773 : : !$OMP TARGET ENTER DATA MAP(alloc:zarot) IF(gpu_option==ABI_GPU_OPENMP)
774 : : !$OMP TARGET UPDATE TO(zarot) IF(gpu_option==ABI_GPU_OPENMP)
775 : : #else
776 : : ABI_UNUSED((/m1,m2/))
777 : : #endif
778 : :
779 : : !========= Case nspinor ==1 ========================
780 : :
781 [ + + ]: 87632 : if (nspinor == 1) then
782 : :
783 [ - + + - : 434015 : ABI_MALLOC(gloc_tmp,(ndim_max,ndim_max*nsppol,nsym))
- + + - -
+ + - + -
+ - - + -
+ ]
784 : : #ifdef HAVE_OPENMP_OFFLOAD
785 : : !$OMP TARGET ENTER DATA MAP(alloc:gloc_tmp) IF(gpu_option==ABI_GPU_OPENMP)
786 : : #endif
787 : :
788 : 86803 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),glocsym(:),gpu_option=gpu_option)
789 : :
790 [ + + ]: 225153 : do iatom=1,natom
791 : :
792 : 138350 : lpawu = gloc(iatom)%lpawu
793 [ + + ]: 138350 : if (lpawu == -1) cycle
794 : 88683 : ndim = 2*lpawu + 1
795 : 88683 : glocsym_mat => glocsym(iatom)%mat
796 : :
797 [ - + + - : 443415 : ABI_MALLOC(gloc_tmp2,(ndim,ndim*nsppol,nsym))
- + + - -
+ + - + -
+ - - + -
+ ]
798 : : #ifdef HAVE_OPENMP_OFFLOAD
799 : : !$OMP TARGET ENTER DATA MAP(alloc:gloc_tmp2) IF(gpu_option==ABI_GPU_OPENMP)
800 : : #endif
801 : :
802 : :
803 [ + - ]: 88683 : if(gpu_option==ABI_GPU_DISABLED) then
804 : :
805 [ + + ]: 4241643 : do irot=1,nsym
806 : 4152960 : at_indx = paw_dmft%indsym(irot,iatom)
807 : 4152960 : gloc_mat => gloc(at_indx)%mat
808 : :
809 [ + + ]: 11637504 : do isppol=1,nsppol
810 : : call abi_xgemm("n","n",ndim,ndim,ndim,cone,gloc_mat(:,:,isppol),ndim,&
811 : 11637504 : & zarot(:,1:ndim,irot,lpawu+1),ndim_max,czero,gloc_tmp(:,1+ndim*(isppol-1):ndim*isppol,irot),ndim_max)
812 : : end do ! isppol
813 : :
814 : : call abi_xgemm("t","n",ndim,ndim*nsppol,ndim,cone,zarot(:,1:ndim,irot,lpawu+1),ndim_max,&
815 : 4152960 : & gloc_tmp(:,:,irot),ndim_max,czero,gloc_tmp2(:,:,irot),ndim)
816 : :
817 [ + + ]: 11726187 : do isppol=1,nsppol
818 [ + - + + : 461281152 : glocsym_mat(:,:,isppol) = glocsym_mat(:,:,isppol) + gloc_tmp2(:,1+ndim*(isppol-1):ndim*isppol,irot)
+ + + + +
+ ]
819 : : end do ! isppol
820 : :
821 : : end do ! irot
822 : :
823 [ + + + + : 4984541 : glocsym_mat(:,:,:) = glocsym_mat(:,:,:) / dble(nsym)
+ + ]
824 : :
825 : : else if(gpu_option==ABI_GPU_OPENMP) then
826 : : #ifdef HAVE_OPENMP_OFFLOAD
827 : :
828 : : do irot=1,nsym
829 : : at_indx = paw_dmft%indsym(irot,iatom)
830 : : gloc_mat => gloc(at_indx)%mat
831 : :
832 : : !$OMP TARGET DATA USE_DEVICE_ADDR(gloc_tmp,zarot,gloc_mat)
833 : : call abi_gpu_xgemm_strided(2,"n","n",ndim,ndim,ndim,cone,&
834 : : & c_loc(gloc_mat(:,:,:)),ndim,ndim*ndim,&
835 : : & c_loc(zarot(:,1:ndim,irot,lpawu+1)),ndim_max,0,czero,&
836 : : & c_loc(gloc_tmp(:,:,irot)),ndim_max,ndim*ndim,nsppol,async=.true.,stream_id=irot)
837 : : !$OMP END TARGET DATA
838 : : end do ! irot
839 : : call gpu_device_synchronize()
840 : :
841 : : !$OMP TARGET DATA USE_DEVICE_ADDR(gloc_tmp,zarot,gloc_tmp2)
842 : : call abi_gpu_xgemm_strided(2,"t","n",ndim,ndim*nsppol,ndim,cone,&
843 : : & c_loc(zarot(:,:,:,lpawu+1)),ndim_max,ndim_max*ndim_max,&
844 : : & c_loc(gloc_tmp(:,:,:)),ndim_max,ndim_max*ndim_max*nsppol,czero,&
845 : : & c_loc(gloc_tmp2(:,:,:)),ndim,ndim*ndim*nsppol,nsym)
846 : : !$OMP END TARGET DATA
847 : :
848 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:glocsym_mat,gloc_tmp2) PRIVATE(isppol)
849 : : do isppol=1,nsppol
850 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(m1,m2,irot)
851 : : do m2=1,ndim
852 : : do m1=1,ndim
853 : : do irot=1,nsym
854 : : glocsym_mat(m1,m2,isppol) = glocsym_mat(m1,m2,isppol) + gloc_tmp2(m1,m2+ndim*(isppol-1),irot)
855 : : end do ! irot
856 : : end do ! m1
857 : : end do ! m2
858 : : end do ! isppol
859 : :
860 : : call abi_xscal(ndim*ndim*nsppol, ratio, glocsym_mat, 1, gpu_option=gpu_option)
861 : :
862 : : #endif
863 : : end if
864 : :
865 : : #ifdef HAVE_OPENMP_OFFLOAD
866 : : !$OMP TARGET EXIT DATA MAP(delete:gloc_tmp2) IF(gpu_option==ABI_GPU_OPENMP)
867 : : #endif
868 [ - + ]: 225153 : ABI_FREE(gloc_tmp2)
869 : :
870 : : end do ! iatom
871 : :
872 : : !== Put glocsym into gloc
873 : 86803 : call copy_matlu(glocsym(:),gloc(:),natom)
874 : :
875 : : #ifdef HAVE_OPENMP_OFFLOAD
876 : : !$OMP TARGET EXIT DATA MAP(delete:gloc_tmp) IF(gpu_option==ABI_GPU_OPENMP)
877 : : #endif
878 [ - + ]: 86803 : ABI_FREE(gloc_tmp)
879 : : !========= Case nspinor ==2 ========================
880 : :
881 : : else
882 : :
883 : 829 : symrec_cart => paw_dmft%symrec_cart(:,:,:)
884 : : #ifdef HAVE_OPENMP_OFFLOAD
885 : : !$OMP TARGET ENTER DATA MAP(to:symrec_cart) IF(gpu_option==ABI_GPU_OPENMP)
886 : : #endif
887 : : !== Allocate temporary arrays
888 [ - + + - : 2487 : ABI_MALLOC(gloc_nmrep,(natom))
- + - + +
+ ]
889 : 829 : call init_matlu(natom,1,4*nsppol,paw_dmft%lpawu(:),glocsym(:),gpu_option=gpu_option)
890 : 829 : call init_matlu(natom,1,4*nsppol,paw_dmft%lpawu(:),gloc_nmrep(:),gpu_option=gpu_option)
891 : :
892 : : ! Put gloc into gloc_nmrep (density and magnetization representation)
893 : : ! gloc_nmrep(iatom)%mat(:,:,i) = n,mx,my,mz for i=1,2,3,4 respectively
894 : 829 : call chg_repr_matlu(gloc(:),gloc_nmrep(:),natom,1,1)
895 : :
896 : : !== Do the sum over symmetrized density matrix (in n,m repr)
897 [ + + ]: 1658 : do iatom=1,natom
898 : :
899 : 829 : lpawu = gloc(iatom)%lpawu
900 [ - + ]: 829 : if (lpawu == -1) cycle
901 : 829 : ndim = 2*lpawu + 1
902 : 829 : glocsym_mat => glocsym(iatom)%mat
903 : :
904 : :
905 [ - + + - : 5803 : ABI_MALLOC(gloc_tmp3,(ndim,ndim,4*nsppol,nsym))
- + + - +
- - + + -
- + + - +
- + - - +
- + ]
906 [ - + + - : 4145 : ABI_MALLOC(gloc_tmp4,(ndim,ndim,4*nsppol,nsym))
- + + - -
+ + - - +
+ - - + -
+ ]
907 : :
908 [ + - ]: 829 : if(gpu_option==ABI_GPU_DISABLED) then
909 [ + + ]: 20725 : do irot=1,nsym
910 : :
911 : 19896 : at_indx = paw_dmft%indsym(irot,iatom)
912 : 19896 : gloc_mat => gloc_nmrep(at_indx)%mat
913 : :
914 [ + + ]: 40621 : do isppol=1,nsppol
915 [ + + ]: 119376 : do mu=1,4 ! Symmetrize density and magnetization
916 : :
917 : : call abi_xgemm("n","n",ndim,ndim,ndim,cone,gloc_mat(:,:,mu+(isppol-1)*4),ndim, &
918 : 99480 : & zarot(:,1:ndim,irot,lpawu+1),ndim_max,czero,gloc_tmp3(:,:,mu+(isppol-1)*4,irot),ndim)
919 : :
920 : : end do ! mu
921 : : end do ! isppol
922 : : end do ! irot
923 : :
924 [ + + ]: 20725 : do irot=1,nsym
925 : : call abi_zgemm_2dd("t","n",ndim,ndim*4*nsppol,ndim,cone,zarot(:,1:ndim,irot,lpawu+1),ndim_max,&
926 : 20725 : & gloc_tmp3(:,:,:,irot),ndim,czero,gloc_tmp4(:,:,:,irot),ndim)
927 : :
928 : : end do ! irot
929 : :
930 [ + + ]: 20725 : do irot=1,nsym
931 [ + + ]: 40621 : do isppol=1,nsppol
932 [ + - + + : 2288040 : glocsym_mat(:,:,1+(isppol-1)*4) = glocsym_mat(:,:,1+(isppol-1)*4) + gloc_tmp4(:,:,1+(isppol-1)*4,irot)
+ + + + +
+ ]
933 : : end do ! isppol
934 : : end do ! irot
935 : :
936 : : ! Symmetrize magnetization
937 : :
938 [ + + ]: 20725 : do irot=1,nsym
939 [ + + ]: 40621 : do isppol=1,nsppol
940 [ + + ]: 99480 : do nu=2,4
941 [ + + ]: 258648 : do mu=2,4
942 : : glocsym_mat(:,:,mu+(isppol-1)*4) = glocsym_mat(:,:,mu+(isppol-1)*4) + &
943 [ + - + + : 20472984 : & symrec_cart(mu-1,nu-1,irot)*gloc_tmp4(:,:,nu+(isppol-1)*4,irot)
+ + + + +
+ ]
944 : : end do ! mu
945 : : end do ! nu
946 : : end do ! isppol
947 : : end do ! irot
948 : :
949 : : ! == Normalize sum
950 [ + + + + : 189841 : glocsym_mat(:,:,:) = glocsym_mat(:,:,:) / dble(nsym)
+ + ]
951 : :
952 : : else if(gpu_option==ABI_GPU_OPENMP) then
953 : : #ifdef HAVE_OPENMP_OFFLOAD
954 : : !$OMP TARGET ENTER DATA MAP(alloc:gloc_tmp3,gloc_tmp4)
955 : :
956 : : do irot=1,nsym
957 : :
958 : : at_indx = paw_dmft%indsym(irot,iatom)
959 : : gloc_mat => gloc_nmrep(at_indx)%mat
960 : :
961 : : !$OMP TARGET DATA USE_DEVICE_ADDR(gloc_tmp3,zarot,gloc_mat)
962 : : call abi_gpu_xgemm_strided(2,"n","n",ndim,ndim,ndim,cone,&
963 : : & c_loc(gloc_mat(:,:,:)),ndim,ndim*ndim,&
964 : : & c_loc(zarot(:,1:ndim,irot,lpawu+1)),ndim_max,0,czero,&
965 : : & c_loc(gloc_tmp3(:,:,:,irot)),ndim,ndim*ndim,4*nsppol,async=.true.,stream_id=irot)
966 : : !$OMP END TARGET DATA
967 : : end do ! irot
968 : : call gpu_device_synchronize()
969 : :
970 : :
971 : : !$OMP TARGET DATA USE_DEVICE_ADDR(gloc_tmp3,zarot,gloc_tmp4)
972 : : call abi_gpu_xgemm_strided(2,"t","n",ndim,ndim*4*nsppol,ndim,cone,&
973 : : & c_loc(zarot(:,:,:,lpawu+1)),ndim_max,ndim_max*ndim_max,&
974 : : & c_loc(gloc_tmp3(:,:,:,:)),ndim,ndim*ndim*4*nsppol,czero,&
975 : : & c_loc(gloc_tmp4(:,:,:,:)),ndim,ndim*ndim*4*nsppol,nsym)
976 : : !$OMP END TARGET DATA
977 : :
978 : :
979 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:glocsym_mat,gloc_tmp4) PRIVATE(isppol)
980 : : do isppol=1,nsppol
981 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(m1,m2,irot)
982 : : do m2=1,ndim
983 : : do m1=1,ndim
984 : : do irot=1,nsym
985 : : glocsym_mat(m1,m2,1+(isppol-1)*4) = glocsym_mat(m1,m2,1+(isppol-1)*4) + gloc_tmp4(m1,m2,1+(isppol-1)*4,irot)
986 : : end do ! irot
987 : : end do ! m1
988 : : end do ! m2
989 : : end do ! isppol
990 : :
991 : : ! Symmetrize magnetization
992 : :
993 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:glocsym_mat,gloc_tmp4,symrec_cart) PRIVATE(isppol,irot)
994 : : do isppol=1,nsppol
995 : : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(nu,mu,m2,m1)
996 : : do mu=2,4
997 : : do m2=1,ndim
998 : : do m1=1,ndim
999 : : do irot=1,nsym
1000 : : do nu=2,4
1001 : : glocsym_mat(m1,m2,mu+(isppol-1)*4) = glocsym_mat(m1,m2,mu+(isppol-1)*4) + &
1002 : : & symrec_cart(mu-1,nu-1,irot)*gloc_tmp4(m1,m2,nu+(isppol-1)*4,irot)
1003 : : end do ! m1
1004 : : end do ! m2
1005 : : end do ! mu
1006 : : end do ! nu
1007 : : end do ! isppol
1008 : : end do ! irot
1009 : :
1010 : : ! == Normalize sum
1011 : : call abi_xscal(ndim*ndim*4*nsppol, ratio, glocsym_mat, 1, gpu_option=gpu_option)
1012 : :
1013 : : !$OMP TARGET EXIT DATA MAP(delete:gloc_tmp3,gloc_tmp4)
1014 : : #endif
1015 : : end if ! gpu_option
1016 : :
1017 [ - + ]: 829 : ABI_FREE(gloc_tmp3)
1018 [ - + ]: 1658 : ABI_FREE(gloc_tmp4)
1019 : :
1020 : : end do ! iatom
1021 : :
1022 : : !== Compute back density matrix in upup dndn updn dnup representation
1023 : 829 : call chg_repr_matlu(gloc(:),glocsym(:),natom,-1,1)
1024 : :
1025 : 829 : call destroy_matlu(gloc_nmrep(:),natom)
1026 [ + - + + : 1658 : ABI_FREE(gloc_nmrep)
- + - + ]
1027 : :
1028 : : #ifdef HAVE_OPENMP_OFFLOAD
1029 : : !$OMP TARGET EXIT DATA MAP(delete:symrec_cart) IF(gpu_option==ABI_GPU_OPENMP)
1030 : : #endif
1031 : : !==============end of nspinor=2 case ===========
1032 : : end if ! nspinor
1033 : :
1034 : 87632 : call destroy_matlu(glocsym(:),natom)
1035 [ + - + + : 226811 : ABI_FREE(glocsym)
- + - + ]
1036 : : #ifdef HAVE_OPENMP_OFFLOAD
1037 : : !$OMP TARGET EXIT DATA MAP(delete:zarot) IF(gpu_option==ABI_GPU_OPENMP)
1038 : : #endif
1039 : :
1040 : : !mt2g(1)=1
1041 : : !mt2g(2)=2
1042 : : !mt2g(3)=4
1043 : : !mx2my2d=5
1044 : : !t2g=paw_dmft%dmftqmc_t2g
1045 : : !x2my2d=paw_dmft%dmftqmc_x2my2d
1046 : :
1047 : : ! DBG_ENTER("COLL")
1048 : :
1049 : : !ci=cone
1050 : : !nspinor=gloc(1)%nspinor
1051 : : !nsppol=gloc(1)%nsppol
1052 : : !natom=cryst_struc%natom
1053 : :
1054 : : !ABI_MALLOC(glocnm,(natom))
1055 : : !ABI_MALLOC(glocnms,(natom))
1056 : : !ABI_MALLOC(glocsym,(natom))
1057 : : !ABI_MALLOC(lpawu_natom,(natom))
1058 : : !lpawu_natom(1:natom)=gloc(1:natom)%lpawu ! If gloc(1:natom)%lpawu is directly used in the next three lines, warnings are generated by some compilers.
1059 : : !call init_matlu(natom,nspinor,nsppol,lpawu_natom,glocnm)
1060 : : !call init_matlu(natom,nspinor,nsppol,lpawu_natom,glocnms)
1061 : : !call init_matlu(natom,nspinor,nsppol,lpawu_natom,glocsym)
1062 : : !ABI_FREE(lpawu_natom)
1063 : :
1064 : : !========= Case nspinor ==1 ========================
1065 : :
1066 : : !if (nspinor==1) then
1067 : : ! ispinor=1
1068 : : ! ispinor1=1
1069 : : ! do iatom=1,cryst_struc%natom
1070 : : ! do isppol=1,nsppol
1071 : : ! if(gloc(iatom)%lpawu/=-1) then
1072 : : ! lpawu=gloc(iatom)%lpawu
1073 : : ! do m1=1, 2*lpawu+1
1074 : : ! do m2=1, 2*lpawu+1
1075 : : ! do irot=1,cryst_struc%nsym
1076 : : ! at_indx=cryst_struc%indsym(4,irot,iatom)
1077 : : ! do m3=1, 2*lpawu+1
1078 : : ! do m4=1, 2*lpawu+1
1079 : : ! if(t2g==1) then
1080 : : ! m1s=mt2g(m1)
1081 : : ! m2s=mt2g(m2)
1082 : : ! m3s=mt2g(m3)
1083 : : ! m4s=mt2g(m4)
1084 : : ! lpawu_zarot=2
1085 : : ! else if (x2my2d==1) then
1086 : : ! m1s=mx2my2d
1087 : : ! m2s=mx2my2d
1088 : : ! m3s=mx2my2d
1089 : : ! m4s=mx2my2d
1090 : : ! lpawu_zarot=2
1091 : : ! else
1092 : : ! m1s=m1
1093 : : ! m2s=m2
1094 : : ! m3s=m3
1095 : : ! m4s=m4
1096 : : ! lpawu_zarot=lpawu
1097 : : ! endif
1098 : : ! zarot2=pawang%zarot(m3s,m1s,lpawu_zarot+1,irot)*pawang%zarot(m4s,m2s,lpawu_zarot+1,irot)
1099 : : ! glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)=&
1100 : : !& glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)&
1101 : : !& +gloc(at_indx)%mat(m3,m4,isppol,ispinor,ispinor1)*zarot2
1102 : : ! end do ! m3
1103 : : ! end do ! m4
1104 : : ! end do ! irot
1105 : : ! glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)=&
1106 : : !& glocsym(iatom)%mat(m1,m2,isppol,ispinor,ispinor1)/real(cryst_struc%nsym,kind=dp)
1107 : : ! end do ! m2
1108 : : ! end do ! m1
1109 : : ! endif ! lpawu/=-1
1110 : : ! end do ! isppol
1111 : : ! end do ! iatom
1112 : : !== Put glocsym into gloc
1113 : : ! do iatom=1,cryst_struc%natom
1114 : : ! if(gloc(iatom)%lpawu/=-1) then
1115 : : ! gloc(iatom)%mat=glocsym(iatom)%mat
1116 : : ! gloc(iatom)%mat(:,:,1,:,:)=(glocsym(iatom)%mat(:,:,1,:,:) &
1117 : : !& + glocsym(iatom)%mat(:,:,2,:,:))/two
1118 : : ! gloc(iatom)%mat(:,:,2,:,:)= gloc(iatom)%mat(:,:,1,:,:)
1119 : : ! write(std_out,*) "WARNING: SYM non mag"
1120 : : ! write(ab_out,*) "WARNING: SYM non mag"
1121 : : ! endif
1122 : : ! end do ! iatom
1123 : :
1124 : : !========= Case nspinor ==2 ========================
1125 : :
1126 : : ! else if (nspinor==2) then
1127 : :
1128 : : !== Allocate temporary arrays
1129 : : ! do iatom=1,cryst_struc%natom
1130 : : ! if(gloc(iatom)%lpawu/=-1) then
1131 : : ! ndim=2*gloc(iatom)%lpawu+1
1132 : : ! ABI_FREE(glocnm(iatom)%mat)
1133 : : ! ABI_FREE(glocnms(iatom)%mat)
1134 : : ! ABI_FREE(glocsym(iatom)%mat)
1135 : : ! ABI_MALLOC(glocnm(iatom)%mat,(ndim,ndim,nsppol,4,1))
1136 : : ! ABI_MALLOC(glocnms(iatom)%mat,(ndim,ndim,nsppol,4,1))
1137 : : ! ABI_MALLOC(glocsym(iatom)%mat,(ndim,ndim,nsppol,2,2))
1138 : : ! endif
1139 : : ! enddo
1140 : : ! ABI_MALLOC(symrec_cart,(3,3,cryst_struc%nsym))
1141 : :
1142 : : !== Compute symrec_cart
1143 : : ! do irot=1,cryst_struc%nsym
1144 : : ! call symredcart(cryst_struc%gprimd,cryst_struc%rprimd,symrec_cart(:,:,irot),cryst_struc%symrec(:,:,irot))
1145 : : ! end do
1146 : :
1147 : : !== Compute density matrix in density and magnetization representation
1148 : : ! call chg_repr_matlu(gloc,glocnm,cryst_struc%natom,option=1,prtopt=1)
1149 : :
1150 : : !== Do the sum over symetrized density matrix (in n,m repr)
1151 : : ! isppol=1
1152 : : ! do iatom=1,cryst_struc%natom
1153 : : ! if(gloc(iatom)%lpawu/=-1) then
1154 : : ! lpawu=gloc(iatom)%lpawu
1155 : : ! ndim=2*gloc(iatom)%lpawu+1
1156 : : ! do m1=1, 2*lpawu+1
1157 : : ! do m2=1, 2*lpawu+1
1158 : : ! sumrho=czero
1159 : : ! rotmag=czero
1160 : : ! do irot=1,cryst_struc%nsym
1161 : : ! summag=czero
1162 : : ! at_indx=cryst_struc%indsym(4,irot,iatom)
1163 : : ! do m3=1, 2*lpawu+1
1164 : : ! do m4=1, 2*lpawu+1
1165 : : ! if(t2g==1) then
1166 : : ! m1s=mt2g(m1)
1167 : : ! m2s=mt2g(m2)
1168 : : ! m3s=mt2g(m3)
1169 : : ! m4s=mt2g(m4)
1170 : : ! lpawu_zarot=2
1171 : : !! else if (x2my2d==1) then
1172 : : ! m1s=mx2my2d
1173 : : ! m2s=mx2my2d
1174 : : ! m3s=mx2my2d
1175 : : ! m4s=mx2my2d
1176 : : ! lpawu_zarot=2
1177 : : ! else
1178 : : ! m1s=m1
1179 : : ! m2s=m2
1180 : : ! m3s=m3
1181 : : ! m4s=m4
1182 : : ! lpawu_zarot=lpawu
1183 : : ! endif
1184 : : ! zarot2=pawang%zarot(m3s,m2s,lpawu_zarot+1,irot)*pawang%zarot(m4s,m1s,lpawu_zarot+1,irot)
1185 : : ! sumrho=sumrho + glocnm(at_indx)%mat(m4,m3,isppol,1,1) * zarot2
1186 : : ! do mu=1,3
1187 : : ! summag(mu)=summag(mu) + glocnm(at_indx)%mat(m4,m3,isppol,mu+1,1) * zarot2
1188 : : ! enddo
1189 : : ! end do ! m3
1190 : : ! end do !m4
1191 : :
1192 : : ! == special case of magnetization
1193 : : ! do nu=1,3
1194 : : ! do mu=1,3
1195 : : ! rotmag(mu)=rotmag(mu)+symrec_cart(mu,nu,irot)*summag(nu)
1196 : : ! end do
1197 : : ! end do
1198 : : ! write(std_out,'(a,3i4,2x,3(2f10.5,2x))') "rotmag",irot,m1,m2,(rotmag(mu),mu=1,3)
1199 : : ! end do ! irot
1200 : :
1201 : : ! == Normalizes sum
1202 : : ! sumrho=sumrho/real(cryst_struc%nsym,kind=dp)
1203 : : ! sumrho=glocnm(isppol,1,iatom,m1,m2) ! test without sym
1204 : : ! glocnms(iatom)%mat(m1,m2,isppol,1,1)=sumrho
1205 : : ! do mu=1,3
1206 : : ! rotmag(mu)=rotmag(mu)/real(cryst_struc%nsym,kind=dp)
1207 : : ! rotmag(mu)=glocnm(isppol,mu+1,iatom,m1,m2) ! test without sym
1208 : : ! glocnms(iatom)%mat(m1,m2,isppol,mu+1,1)=rotmag(mu)
1209 : : ! enddo
1210 : : ! end do ! m2
1211 : : ! end do ! m1
1212 : : ! endif ! lpawu/=-1
1213 : : ! end do ! iatom
1214 : :
1215 : : !== Compute back density matrix in upup dndn updn dnup representation
1216 : : ! call chg_repr_matlu(glocsym,glocnms,cryst_struc%natom,option=-1,prtopt=1)
1217 : :
1218 : : !== Put glocsym into gloc
1219 : : ! do iatom=1,cryst_struc%natom
1220 : : ! if(gloc(iatom)%lpawu/=-1) then
1221 : : ! gloc(iatom)%mat=glocsym(iatom)%mat
1222 : : ! gloc(iatom)%mat(:,:,1,:,:)=(glocsym(iatom)%mat(:,:,1,:,:) &
1223 : : !& + glocsym(iatom)%mat(:,:,2,:,:))/two
1224 : : ! gloc(iatom)%mat(:,:,2,:,:)= gloc(iatom)%mat(:,:,1,:,:)
1225 : : ! write(std_out,*) "WARNING: SYM non mag"
1226 : : ! write(ab_out,*) "WARNING: SYM non mag"
1227 : : ! endif
1228 : : ! end do ! iatom
1229 : :
1230 : : ! ABI_FREE(symrec_cart)
1231 : : ! endif
1232 : :
1233 : : ! call destroy_matlu(glocnm,cryst_struc%natom)
1234 : : ! call destroy_matlu(glocnms,cryst_struc%natom)
1235 : : ! call destroy_matlu(glocsym,cryst_struc%natom)
1236 : : ! ABI_FREE(glocnm)
1237 : : ! ABI_FREE(glocnms)
1238 : : ! ABI_FREE(glocsym)
1239 : :
1240 [ - + ]: 87632 : end subroutine sym_matlu
1241 : : !!***
1242 : :
1243 : : !!****f* m_matlu/inverse_matlu
1244 : : !! NAME
1245 : : !! inverse_matlu
1246 : : !!
1247 : : !! FUNCTION
1248 : : !! Inverse local quantity.
1249 : : !!
1250 : : !! COPYRIGHT
1251 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
1252 : : !! This file is distributed under the terms of the
1253 : : !! GNU General Public License, see ~abinit/COPYING
1254 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1255 : : !!
1256 : : !! INPUTS
1257 : : !! matlu(natom) :: input quantity to inverse
1258 : : !! natom=number of atoms in cell.
1259 : : !!
1260 : : !! OUTPUT
1261 : : !! matlu(natom) :: inverse of input matrix
1262 : : !!
1263 : : !! SIDE EFFECTS
1264 : : !!
1265 : : !! NOTES
1266 : : !!
1267 : : !! SOURCE
1268 : :
1269 : 48249 : subroutine inverse_matlu(matlu,natom)
1270 : :
1271 : : !Arguments ------------------------------------
1272 : : integer, intent(in) :: natom
1273 : : type(matlu_type), intent(inout) :: matlu(natom)
1274 : : !Local variables-------------------------------
1275 : : integer :: iatom,isppol,lpawu,ndim,nspinor,nsppol
1276 : : !************************************************************************
1277 : :
1278 : 48249 : nspinor = matlu(1)%nspinor
1279 : 48249 : nsppol = matlu(1)%nsppol
1280 : :
1281 [ + + ]: 124066 : do iatom=1,natom
1282 : 75817 : lpawu = matlu(iatom)%lpawu
1283 [ + + ]: 75817 : if (lpawu == -1) cycle
1284 : 49401 : ndim = nspinor * (2*lpawu+1)
1285 [ + + ]: 188712 : do isppol=1,nsppol
1286 : 166879 : call xginv(matlu(iatom)%mat(:,:,isppol),ndim)
1287 : : end do ! isppol
1288 : : end do ! iatom
1289 : :
1290 : : !if(prtopt>0) then
1291 : : !endif
1292 : : !ABI_MALLOC(gathermatlu,(natom))
1293 : : !do iatom=1,natom
1294 : : ! if(matlu(iatom)%lpawu.ne.-1) then
1295 : : ! tndim=nsppol*nspinor*(2*matlu(iatom)%lpawu+1)
1296 : : ! ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
1297 : : ! gathermatlu(iatom)%value=czero
1298 : : ! endif
1299 : : !enddo
1300 : :
1301 : : !call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
1302 : : !do iatom=1,natom
1303 : : ! if(matlu(iatom)%lpawu.ne.-1) then
1304 : : ! tndim=nsppol*nspinor*(2*matlu(iatom)%lpawu+1)
1305 : : !call matcginv_dpc(gathermatlu(iatom)%value,tndim,tndim)
1306 : : ! call xginv(gathermatlu(iatom)%value,tndim)
1307 : : ! endif
1308 : : !enddo
1309 : : !call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
1310 : :
1311 : : !do iatom=1,natom
1312 : : ! if(matlu(iatom)%lpawu.ne.-1) then
1313 : : ! ABI_FREE(gathermatlu(iatom)%value)
1314 : : ! endif
1315 : : !enddo
1316 : : !ABI_FREE(gathermatlu)
1317 : :
1318 : 48249 : end subroutine inverse_matlu
1319 : : !!***
1320 : :
1321 : : !!****f* m_matlu/diff_matlu
1322 : : !! NAME
1323 : : !! diff_matlu
1324 : : !!
1325 : : !! FUNCTION
1326 : : !!
1327 : : !! INPUTS
1328 : : !! char1 = character describing matlu1
1329 : : !! char2 = character describing matlu2
1330 : : !! matlu1(natom) <type(matlu_type)>= density matrix 1 in the local orbital basis and related variables
1331 : : !! matlu2(natom) <type(matlu_type)>= density matrix 2 in the local orbital basis and related variables
1332 : : !! natom = number of atoms
1333 : : !! option =1 if diff > toldiff , stop
1334 : : !! 0 print diff and toldiff
1335 : : !! else do not test and do not print
1336 : : !! toldiff = threshold for the difference between matlu1 and matlu2
1337 : : !! zero_or_one = useful when comparing the overlap of Wannier functions at one kpt with the identity
1338 : : !!
1339 : : !! OUTPUT
1340 : : !! ierr = 0 if diff < toldiff
1341 : : !! = -1 otherwise
1342 : : !!
1343 : : !! SOURCE
1344 : :
1345 : 926 : subroutine diff_matlu(char1,char2,matlu1,matlu2,natom,option,toldiff,ierr,zero_or_one)
1346 : :
1347 : : !Arguments ------------------------------------
1348 : : integer, intent(in) :: natom,option
1349 : : type(matlu_type), intent(in) :: matlu1(natom),matlu2(natom)
1350 : : character(len=*), intent(in) :: char1,char2
1351 : : real(dp), intent(in) :: toldiff
1352 : : integer, optional, intent(out) :: ierr
1353 : : integer, optional, intent(in) :: zero_or_one
1354 : : !Local variables-------------------------------
1355 : : integer :: iatom,idiff,lpawu,nspinor,nsppol
1356 : : real(dp) :: matludiff
1357 : : character(len=500) :: message
1358 : : ! *********************************************************************
1359 : :
1360 [ + + ]: 926 : if (option /= 1 .and. option /= 0) return
1361 : :
1362 : 918 : idiff = 0
1363 : 918 : matludiff = zero
1364 : 918 : nspinor = matlu1(1)%nspinor
1365 : 918 : nsppol = matlu1(1)%nsppol
1366 : :
1367 [ + + ]: 3984 : do iatom=1,natom
1368 : 3066 : lpawu = matlu1(iatom)%lpawu
1369 [ + + ]: 3066 : if (lpawu == -1) cycle
1370 [ + + + + : 45413 : matludiff = matludiff + sum(abs(matlu1(iatom)%mat(:,:,:)-matlu2(iatom)%mat(:,:,:)))
+ + ]
1371 : 3984 : idiff = idiff + (2*lpawu+1)**2
1372 : : end do ! iatom
1373 : 918 : idiff = idiff * (nspinor**2) * nsppol
1374 : :
1375 [ + + ]: 918 : if (.not. present(zero_or_one)) matludiff = matludiff / dble(idiff)
1376 : :
1377 [ + + ]: 918 : if (matludiff < toldiff) then
1378 : 845 : write(message,'(5a,6x,3a,4x,e12.4,a,e12.4)') ch10,' ** Differences between ',trim(char1),' and ',&
1379 : 1690 : & ch10,trim(char2),' are small enough:',ch10,matludiff,' is lower than',toldiff
1380 : 845 : call wrtout(std_out,message,'COLL')
1381 [ + + ]: 845 : if (present(ierr)) ierr = 0
1382 : : else
1383 : 73 : write(message,'(5a,3x,3a,3x,e12.4,a,e12.4)') ch10,'Differences between ',trim(char1),' and ',&
1384 : 146 : & ch10,trim(char2),' is too large:',ch10,matludiff,' is larger than',toldiff
1385 : 73 : ABI_WARNING(message)
1386 : : ! write(message,'(8a,4x,e12.4,a,e12.4)') ch10," Matrix for ",trim(char1)
1387 : 73 : write(message,'(a,3x,a)') ch10,trim(char1)
1388 : 73 : call wrtout(std_out,message,'COLL')
1389 : 73 : call print_matlu(matlu1(:),natom,1,opt_diag=-1)
1390 : 73 : write(message,'(a,3x,a)') ch10,trim(char2)
1391 : 73 : call wrtout(std_out,message,'COLL')
1392 : 73 : call print_matlu(matlu2(:),natom,1,opt_diag=-1)
1393 [ - + - - ]: 73 : if (present(zero_or_one) .and. (mod(matludiff,one) < toldiff)) then
1394 : 0 : write(message,'(a,3x,a)') ch10," The norm is not identity for this k-point &
1395 : 0 : & but is compatible with a high symmetry point"
1396 : 0 : call wrtout(std_out,message,'COLL')
1397 : : else if (present(zero_or_one)) then
1398 : 0 : write(message,'(a,3x,a)') ch10," The norm is not identity for this k-point but might be compatible &
1399 : 0 : & with a high symmetry point: it should be checked"
1400 : 0 : call wrtout(std_out,message,'COLL')
1401 [ - + ]: 73 : else if (option == 1) then
1402 : 0 : call flush_unit(std_out)
1403 : 0 : write(message,'(5a,6x,3a,4x,e12.4,a,e12.4)') ch10,' ** Differences between ',trim(char1),' and ',&
1404 : 0 : & ch10,trim(char2),' are too high:',ch10,matludiff,' is greater than',toldiff
1405 : 0 : ABI_ERROR(message)
1406 : : end if ! zero_or_one
1407 [ + + ]: 73 : if (present(ierr)) ierr = -1
1408 : : end if ! matludiff < toldiff
1409 : :
1410 : : end subroutine diff_matlu
1411 : : !!***
1412 : :
1413 : : !!****f* m_matlu/add_matlu
1414 : : !! NAME
1415 : : !! add_matlu
1416 : : !!
1417 : : !! FUNCTION
1418 : : !!
1419 : : !! INPUTS
1420 : : !! matlu1 <type(matlu_type)>= density matrix matlu1 in the local orbital basis and related variables
1421 : : !! matlu2 <type(matlu_type)>= density matrix matlu2 in the local orbital basis and related variables
1422 : : !! natom = number of atoms
1423 : : !! sign_matlu2= 1 add matlu1 and matlu2
1424 : : !! -1 substract matlu2 to matlu1
1425 : : !!
1426 : : !! OUTPUT
1427 : : !! matlu3 <type(matlu_type)>= density matrix matlu3, sum/substract matlu1 and matlu2
1428 : : !!
1429 : : !! SOURCE
1430 : :
1431 : 974058 : subroutine add_matlu(matlu1,matlu2,matlu3,natom,sign_matlu2,idat,ndat)
1432 : :
1433 : : !Arguments ------------------------------------
1434 : : integer, intent(in) :: natom,sign_matlu2
1435 : : integer, optional,intent(in) :: idat,ndat
1436 : : type(matlu_type), intent(in) :: matlu1(natom),matlu2(natom)
1437 : : type(matlu_type), intent(inout) :: matlu3(natom) !vz_i
1438 : : !Local variables-------------------------------
1439 : : integer :: iatom,lpawu,isppol
1440 : : ! *********************************************************************
1441 : :
1442 [ + - ]: 974058 : if(present(idat) .and. present(ndat)) then
1443 [ # # ]: 0 : do iatom=1,natom
1444 : 0 : lpawu = matlu1(iatom)%lpawu
1445 [ # # ]: 0 : if (lpawu == -1) cycle
1446 [ # # ]: 0 : do isppol=1,matlu1(iatom)%nsppol
1447 [ # # ]: 0 : if (sign_matlu2 == 1) then
1448 [ # # # # ]: 0 : matlu3(iatom)%mat(:,:,(isppol-1)*ndat+idat) = matlu1(iatom)%mat(:,:,isppol) + matlu2(iatom)%mat(:,:,isppol)
1449 [ # # ]: 0 : else if (sign_matlu2 == -1) then
1450 [ # # # # ]: 0 : matlu3(iatom)%mat(:,:,(isppol-1)*ndat+idat) = matlu1(iatom)%mat(:,:,isppol) - matlu2(iatom)%mat(:,:,isppol)
1451 : : end if
1452 : : end do ! isppol
1453 : : end do ! iatom
1454 : : else
1455 [ + + ]: 3085443 : do iatom=1,natom
1456 : 2111385 : lpawu = matlu1(iatom)%lpawu
1457 [ + + ]: 2111385 : if (lpawu == -1) cycle
1458 [ + + ]: 1968768 : if (sign_matlu2 == 1) then
1459 [ + + + + : 489568 : matlu3(iatom)%mat(:,:,:) = matlu1(iatom)%mat(:,:,:) + matlu2(iatom)%mat(:,:,:)
+ + ]
1460 [ + - ]: 986294 : else if (sign_matlu2 == -1) then
1461 [ + + + + : 50089479 : matlu3(iatom)%mat(:,:,:) = matlu1(iatom)%mat(:,:,:) - matlu2(iatom)%mat(:,:,:)
+ + ]
1462 : : end if
1463 : : end do ! iatom
1464 : : end if
1465 : :
1466 : 974058 : end subroutine add_matlu
1467 : : !!***
1468 : :
1469 : : !!****f* m_matlu/chg_repr_matlu
1470 : : !! NAME
1471 : : !! chg_repr_matlu
1472 : : !!
1473 : : !! FUNCTION
1474 : : !! Change representation of density matrix (useful for nspinor=2)
1475 : : !!
1476 : : !! COPYRIGHT
1477 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
1478 : : !! This file is distributed under the terms of the
1479 : : !! GNU General Public License, see ~abinit/COPYING
1480 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1481 : : !!
1482 : : !! INPUTS
1483 : : !! glocspsp(natom) :: density matrix in the spin spin representation
1484 : : !! glocnm(natom) :: density matrix in the magnetization representation
1485 : : !! natom=number of atoms in cell.
1486 : : !! option= 1 glocspsp is input, glocnm is computed
1487 : : !! = -1 glocspsp is computed, glocnm is input
1488 : : !! prtopt= abs(prtopt) >= 3 : print in magnetization representation
1489 : : !!
1490 : : !! OUTPUT
1491 : : !! glocspsp(natom) :: density matrix in the spin spin representation
1492 : : !! glocnm(natom) :: density matrix in the magnetization representation
1493 : : !!
1494 : : !! SIDE EFFECTS
1495 : : !!
1496 : : !! NOTES
1497 : : !!
1498 : : !! SOURCE
1499 : :
1500 : 1658 : subroutine chg_repr_matlu(glocspsp,glocnm,natom,option,prtopt)
1501 : :
1502 : : !Arguments ------------------------------------
1503 : : integer, intent(in) :: natom,option,prtopt
1504 : : type(matlu_type), intent(inout), target :: glocnm(natom),glocspsp(natom)
1505 : : !Local variables-------------------------------
1506 : : integer :: iatom,lpawu,m1,m2,mu,ndim,nsppol,isppol,gpu_option
1507 : 1658 : complex(dp), ABI_CONTIGUOUS pointer :: glocnm_mat(:,:,:),glocspsp_mat(:,:,:)
1508 : : character(len=500) :: message
1509 : :
1510 : : ! DBG_ENTER("COLL")
1511 : 0 : ABI_CHECK(glocspsp(1)%nsppol*4==glocnm(1)%nsppol, "Mismatch in nsppol between glocspsp and glocnm")
1512 [ - + ]: 1658 : ABI_CHECK(glocspsp(1)%gpu_option==glocnm(1)%gpu_option, "Mismatch in gpu_option between glocspsp and glocnm")
1513 : 1658 : nsppol=glocspsp(1)%nsppol
1514 : 1658 : gpu_option=glocspsp(1)%gpu_option
1515 : :
1516 : : !== Compute density matrix in density magnetization representation
1517 [ + + ]: 1658 : if (option == 1) then
1518 [ + + ]: 1658 : do iatom=1,natom
1519 : 829 : lpawu = glocspsp(iatom)%lpawu
1520 [ + - ]: 829 : if (lpawu == -1) cycle
1521 : 829 : ndim = 2*lpawu + 1
1522 : 829 : glocnm_mat => glocnm(iatom)%mat
1523 : 829 : glocspsp_mat => glocspsp(iatom)%mat
1524 : :
1525 : : #ifdef HAVE_OPENMP_OFFLOAD
1526 : : !$OMP TARGET TEAMS DISTRIBUTE PRIVATE(isppol) MAP(to:glocnm_mat,glocspsp_mat) &
1527 : : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1528 : : #endif
1529 [ + + ]: 1658 : do isppol=1,nsppol
1530 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(m1,m2)
1531 [ + + ]: 7461 : do m2=1,ndim
1532 [ + + ]: 47253 : do m1=1,ndim
1533 : 40621 : glocnm_mat(m1,m2,1+(isppol-1)*4) = glocspsp_mat(m1,m2, isppol) + glocspsp_mat(m1+ndim,m2+ndim,isppol)
1534 : 40621 : glocnm_mat(m1,m2,4+(isppol-1)*4) = glocspsp_mat(m1,m2, isppol) - glocspsp_mat(m1+ndim,m2+ndim,isppol)
1535 : 40621 : glocnm_mat(m1,m2,2+(isppol-1)*4) = glocspsp_mat(m1,m2+ndim,isppol) + glocspsp_mat(m1+ndim,m2,isppol)
1536 : 46424 : glocnm_mat(m1,m2,3+(isppol-1)*4) = (glocspsp_mat(m1,m2+ndim,isppol) - glocspsp_mat(m1+ndim,m2,isppol)) * j_dpc
1537 : : end do ! m1
1538 : : end do ! m2
1539 : : end do ! isppol
1540 : :
1541 [ - + ]: 1658 : if (abs(prtopt) >= 3) then
1542 : : #ifdef HAVE_OPENMP_OFFLOAD
1543 : : !$OMP TARGET UPDATE FROM(glocnm) IF(gpu_option==ABI_GPU_OPENMP)
1544 : : #endif
1545 : 0 : write(message,'(a)') " -- in n, m repr "
1546 : 0 : call wrtout(std_out,message,'COLL')
1547 [ # # ]: 0 : do mu=1,4
1548 [ # # ]: 0 : do m1=1,ndim
1549 [ # # # # ]: 0 : write(message,'(8x,(14(2f9.5,2x)))') (glocnm(iatom)%mat(m1,m2,isppol+mu-1),m2=1,ndim)
1550 : 0 : call wrtout(std_out,message,'COLL')
1551 : : end do ! m1
1552 : 0 : write(message,'(a)') ch10
1553 : 0 : call wrtout(std_out,message,'COLL')
1554 : : end do ! mu
1555 : : end if ! prtopt >3
1556 : : end do ! iatom
1557 : :
1558 : : !== Compute back density matrix in upup dndn updn dnup representation
1559 [ + - ]: 829 : else if (option == -1) then
1560 : :
1561 [ + + ]: 1658 : do iatom=1,natom
1562 : 829 : lpawu = glocnm(iatom)%lpawu
1563 [ + - ]: 829 : if (lpawu == -1) cycle
1564 : 829 : ndim = 2*lpawu + 1
1565 : 829 : glocnm_mat => glocnm(iatom)%mat
1566 : 829 : glocspsp_mat => glocspsp(iatom)%mat
1567 : :
1568 : : #ifdef HAVE_OPENMP_OFFLOAD
1569 : : !$OMP TARGET TEAMS DISTRIBUTE PRIVATE(isppol) MAP(to:glocnm_mat,glocspsp_mat) &
1570 : : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1571 : : #endif
1572 [ + + ]: 2487 : do isppol=1,nsppol
1573 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(m1,m2)
1574 [ + + ]: 7461 : do m2=1,ndim
1575 [ + + ]: 47253 : do m1=1,ndim
1576 : : glocspsp_mat(m1,m2,isppol) = &
1577 : 40621 : & half * (glocnm_mat(m1,m2,1+(isppol-1)*4)+glocnm_mat(m1,m2,4+(isppol-1)*4))
1578 : : glocspsp_mat(m1+ndim,m2+ndim,isppol) = &
1579 : 40621 : & half * (glocnm_mat(m1,m2,1+(isppol-1)*4)-glocnm_mat(m1,m2,4+(isppol-1)*4))
1580 : : glocspsp_mat(m1,m2+ndim,isppol) = &
1581 : 40621 : & half * (glocnm_mat(m1,m2,2+(isppol-1)*4)-j_dpc*glocnm_mat(m1,m2,3+(isppol-1)*4))
1582 : : glocspsp_mat(m1+ndim,m2,isppol) = &
1583 : 46424 : & half * (glocnm_mat(m1,m2,2+(isppol-1)*4)+j_dpc*glocnm_mat(m1,m2,3+(isppol-1)*4))
1584 : : end do ! m1
1585 : : end do ! m2
1586 : : end do ! isppol
1587 : : !if (abs(prtopt) > 6) then
1588 : : ! write(message,'(a)') " -- in spin spin repr "
1589 : : ! call wrtout(std_out,message,'COLL')
1590 : : ! do mu=1,4
1591 : : ! do m1=1,ndim
1592 : : ! write(message,'(8x,14(2f9.5,2x))') (glocspsp(iatom)%mat(m1,m2,isppol,mu,1),m2=1,ndim)
1593 : : ! call wrtout(std_out, message,'COLL')
1594 : : ! end do ! m1
1595 : : ! write(message,'(a)') ch10
1596 : : ! call wrtout(std_out,message,'COLL')
1597 : : ! end do
1598 : : !end if ! prtopt>6
1599 : : end do ! iatom
1600 : : else
1601 : 0 : message = "stop in chg_repr_matlu"
1602 : 0 : ABI_ERROR(message)
1603 : : end if ! option
1604 : :
1605 : : ! DBG_EXIT("COLL")
1606 : :
1607 [ - + ]: 1658 : end subroutine chg_repr_matlu
1608 : : !!***
1609 : :
1610 : : !!****f* m_matlu/trace_matlu
1611 : : !! NAME
1612 : : !! trace_matlu
1613 : : !!
1614 : : !! FUNCTION
1615 : : !! Compute the trace of the matlu matrix
1616 : : !!
1617 : : !! INPUTS
1618 : : !! matlu(natom) <type(matlu_type)>= density matrix in the
1619 : : !! local orbital basis and related variables
1620 : : !! natom = number of atoms
1621 : : !! itau = flag for print
1622 : : !! not present (default) : occupations from G(iw)
1623 : : !! = 1 : occupations from G(tau)
1624 : : !! = -1 : occupations from G0(tau)
1625 : : !! = 4 : trace of matlu
1626 : : !!
1627 : : !! OUTPUT
1628 : : !! trace_loc(nsppol+1,natom)= trace for each atom and each polarization,
1629 : : !! trace_loc(iatom,nsppol+1) is
1630 : : !! the full trace over all polarizations
1631 : : !! trace= trace over all correlated atoms
1632 : : !!
1633 : : !! SOURCE
1634 : :
1635 : 779 : subroutine trace_matlu(matlu,natom,trace_loc,itau,trace)
1636 : :
1637 : : !Arguments ------------------------------------
1638 : : integer, intent(in) :: natom
1639 : : type(matlu_type), intent(in) :: matlu(natom)
1640 : : real(dp), target, optional, intent(inout) :: trace_loc(matlu(1)%nsppol+1,natom)
1641 : : integer, optional, intent(in) :: itau
1642 : : complex(dp), optional, intent(out) :: trace
1643 : : !Local variables-------------------------------
1644 : : integer :: iatom,im,isppol,lpawu,ndim,nspinor,nsppol
1645 : : complex(dp) :: trace_tmp,trace_tmp2
1646 : : real(dp), ABI_CONTIGUOUS pointer :: traceloc(:,:) => null()
1647 : : character(len=4) :: tag
1648 : : character(len=12) :: tag_nb_elec
1649 : : character(len=500) :: message
1650 : : ! *********************************************************************
1651 : :
1652 : 780 : nspinor = matlu(1)%nspinor
1653 : 780 : nsppol = matlu(1)%nsppol
1654 : :
1655 [ + + ]: 780 : if (present(trace_loc)) then
1656 : 779 : traceloc => trace_loc(:,:)
1657 : : else
1658 [ + - - + : 4 : ABI_MALLOC(traceloc,(nsppol+1,natom))
+ - - + +
- + - + -
- + - + ]
1659 : : end if
1660 : :
1661 : 780 : trace_tmp = czero
1662 [ + + + + ]: 9893 : traceloc(:,:) = zero
1663 : :
1664 [ + + ]: 3489 : do iatom=1,natom
1665 : :
1666 : 2709 : lpawu = matlu(iatom)%lpawu
1667 [ + + ]: 2709 : if (lpawu == -1) cycle
1668 : :
1669 : 912 : ndim = nspinor * (2*lpawu+1)
1670 : 912 : write(tag,'(i4)') iatom
1671 : 912 : write(message,'(3a)') ch10,' -------> For Correlated Atom ',adjustl(tag)
1672 [ + + ]: 912 : if (.not. present(itau)) then
1673 : 911 : call wrtout(std_out,message,'COLL')
1674 : : end if
1675 : :
1676 : : if (present(itau)) then
1677 [ - + ]: 1 : if (itau > 0) then
1678 : 1 : call wrtout(std_out,message,'COLL')
1679 : : end if
1680 : : end if ! present(itau)
1681 : :
1682 [ + + ]: 2293 : do isppol=1,nsppol
1683 : : trace_tmp2 = czero
1684 [ + + ]: 7809 : do im=1,ndim
1685 : 7809 : trace_tmp2 = trace_tmp2 + matlu(iatom)%mat(im,im,isppol)
1686 : : end do ! im
1687 : 1381 : trace_tmp = trace_tmp + trace_tmp2
1688 : 1381 : traceloc(isppol,iatom) = dble(trace_tmp2)
1689 : 2293 : traceloc(nsppol+1,iatom) = traceloc(nsppol+1,iatom) + traceloc(isppol,iatom)
1690 : : end do ! isppol
1691 [ + + ]: 912 : if (nsppol == 1 .and. nspinor == 1) traceloc(nsppol+1,iatom) = traceloc(nsppol+1,iatom) * two
1692 : 912 : write(tag_nb_elec,'(f12.6)') traceloc(nsppol+1,iatom)
1693 : 912 : tag_nb_elec = adjustl(tag_nb_elec)
1694 [ + + ]: 912 : if (.not. present(itau)) then
1695 : 911 : write(message,'(8x,2a)') 'Nb of Corr. elec. from G(iw) is: ',tag_nb_elec
1696 : 911 : call wrtout(std_out,message,'COLL')
1697 : : end if ! not present(itau)
1698 : 780 : if (present(itau)) then
1699 [ + - ]: 1 : if (itau == 1) then
1700 : 1 : write(message,'(8x,2a)') 'Nb of Corr. elec. from G(tau=0-) is: ',tag_nb_elec
1701 : 1 : call wrtout(std_out,message,'COLL')
1702 [ # # ]: 0 : else if (itau == -1) then
1703 : 0 : write(message,'(8x,2a)') 'Nb: Sum of the values of G0(tau=0-) is: ',tag_nb_elec
1704 : 0 : call wrtout(std_out,message,'COLL')
1705 [ # # ]: 0 : else if (itau == 4) then
1706 : 0 : write(message,'(8x,2a)') 'Trace of matlu matrix is: ',tag_nb_elec
1707 : 0 : call wrtout(std_out,message,'COLL')
1708 : : end if ! itau
1709 : : end if ! present(itau)
1710 : : end do ! iatom
1711 : :
1712 [ - + ]: 780 : if (present(trace)) then
1713 [ # # ]: 0 : if (nsppol == 1 .and. nspinor == 1) trace_tmp = trace_tmp * two
1714 : 0 : trace = trace_tmp
1715 : : end if ! present(trace)
1716 : :
1717 [ + + ]: 780 : if (nsppol > 1 .and. (.not. present(trace))) then
1718 [ + + ]: 1323 : do iatom=1,natom
1719 : 986 : lpawu = matlu(iatom)%lpawu
1720 [ + + ]: 986 : if (lpawu == -1) cycle
1721 : : !! MAG
1722 : : ! if(nsppol>1.and.present(itau)) then
1723 : : ! if(itau==1) then
1724 : 469 : write(tag_nb_elec,'(f12.6)') traceloc(2,iatom) - traceloc(1,iatom)
1725 : 469 : write(message,'(8x,2a)') 'DMFT Corr. Elec. Mag.: ',adjustl(tag_nb_elec)
1726 : 1323 : call wrtout(std_out,message,'COLL')
1727 : : ! endif
1728 : : end do ! iatom
1729 : : end if ! nsppol>1
1730 : :
1731 [ + + ]: 780 : if (.not. present(trace_loc)) then
1732 [ - + ]: 1 : ABI_FREE(traceloc)
1733 : : end if
1734 : 780 : traceloc => null()
1735 : :
1736 [ + + ]: 780 : end subroutine trace_matlu
1737 : : !!***
1738 : :
1739 : : !!****f* m_matlu/gather_matlu
1740 : : !! NAME
1741 : : !! gather_matlu
1742 : : !!
1743 : : !! FUNCTION
1744 : : !! Create new array from matlu
1745 : : !!
1746 : : !! COPYRIGHT
1747 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
1748 : : !! This file is distributed under the terms of the
1749 : : !! GNU General Public License, see ~abinit/COPYING
1750 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1751 : : !!
1752 : : !! INPUTS
1753 : : !! gloc(natom) <type(matlu_type)> = density matrix in the spin spin representation
1754 : : !! gatherloc(natom) <type(coeff2c_type)> = density matrix where spin and angular momentum are gathered in the same index
1755 : : !! natom=number of atoms in cell.
1756 : : !! option= 1 go from gloc to gathergloc
1757 : : !! option= -1 go from gathergloc to gloc
1758 : : !! prtopt= option to define level of printing
1759 : : !!
1760 : : !! OUTPUT
1761 : : !! gloc(natom) <type(matlu_type)> = density matrix in the spin spin representation
1762 : : !! gatherloc(natom) <type(coeff2c_type)> = density matrix where spin and angular momentum are gathered in the same index
1763 : : !!
1764 : : !! SIDE EFFECTS
1765 : : !!
1766 : : !! SOURCE
1767 : 0 : subroutine gather_matlu(gloc,gathergloc,natom,option,prtopt)
1768 : :
1769 : : use defs_wvltypes
1770 : : use m_crystal, only : crystal_t
1771 : :
1772 : : ! type matlus_type
1773 : : ! SEQUENCE
1774 : : ! complex(dp), pointer :: mat(:,:)
1775 : : ! end type matlus_type
1776 : :
1777 : : !Arguments ------------------------------------
1778 : : !scalars
1779 : : integer,intent(in) :: natom,option,prtopt
1780 : : type(coeff2c_type), intent(inout) :: gathergloc(natom)
1781 : : type(matlu_type),intent(inout) :: gloc(natom)
1782 : : !Local variables-------------------------------
1783 : : !scalars
1784 : : integer :: iatom,im1,im2,ispinor,ispinor1,isppol,isppol1
1785 : : integer :: jc1,jc2,ml1,ml2,ndim,nspinor,nsppol,tndim
1786 : : character(len=500) :: message
1787 : :
1788 : : ! DBG_ENTER("COLL")
1789 : 0 : nsppol=gloc(1)%nsppol
1790 : 0 : nspinor=gloc(1)%nspinor
1791 : :
1792 [ # # ]: 0 : do iatom=1,natom
1793 [ # # ]: 0 : if(gloc(iatom)%lpawu.ne.-1) then
1794 : : !==-------------------------------------
1795 : :
1796 : 0 : ndim=2*gloc(iatom)%lpawu+1
1797 : 0 : tndim=nsppol*nspinor*ndim
1798 : :
1799 : : !== Put norm into array "gathergloc"
1800 : 0 : jc1=0
1801 [ # # ]: 0 : do isppol=1,nsppol
1802 [ # # ]: 0 : do ispinor=1,nspinor
1803 [ # # ]: 0 : do ml1=1,ndim
1804 : 0 : jc1=jc1+1
1805 : 0 : jc2=0
1806 [ # # ]: 0 : do isppol1=1,nsppol
1807 [ # # ]: 0 : do ispinor1=1,nspinor
1808 [ # # ]: 0 : do ml2=1,ndim
1809 : 0 : jc2=jc2+1
1810 [ # # ]: 0 : if(option==1) then
1811 [ # # ]: 0 : if(isppol==isppol1) then
1812 : 0 : gathergloc(iatom)%value(jc1,jc2)=gloc(iatom)%mat(ml1+(ispinor-1)*ndim,ml2+(ispinor1-1)*ndim,isppol)
1813 : : endif
1814 [ # # ]: 0 : else if(option==-1) then
1815 [ # # ]: 0 : if(isppol==isppol1) then
1816 : 0 : gloc(iatom)%mat(ml1+(ispinor-1)*ndim,ml2+(ispinor1-1)*ndim,isppol)=gathergloc(iatom)%value(jc1,jc2)
1817 : : endif
1818 : : endif
1819 : : enddo
1820 : : enddo ! ispinor1
1821 : : enddo ! isppol1
1822 : : enddo
1823 : : enddo !ispinor
1824 : : enddo ! isppol
1825 : : endif
1826 : : enddo ! iatom
1827 [ # # # # ]: 0 : if(option==1.and.prtopt==3) then
1828 [ # # ]: 0 : do iatom=1,natom
1829 [ # # ]: 0 : if(gloc(iatom)%lpawu.ne.-1) then
1830 : 0 : tndim=nsppol*nspinor*(2*gloc(iatom)%lpawu+1)
1831 : 0 : write(message,'(2a,i5)') ch10,' (gathermatlu:) For atom', iatom
1832 : 0 : call wrtout(std_out,message,'COLL')
1833 [ # # ]: 0 : do im1=1,tndim
1834 : : write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
1835 [ # # # # ]: 0 : & (gathergloc(iatom)%value(im1,im2),im2=1,tndim)
1836 : 0 : call wrtout(std_out,message,'COLL')
1837 : : end do
1838 : : endif
1839 : : enddo ! iatom
1840 [ # # # # ]: 0 : else if(option==-1.and.prtopt==3) then
1841 : 0 : call print_matlu(gloc,natom,prtopt)
1842 : : endif
1843 : :
1844 : :
1845 : :
1846 : : ! DBG_EXIT("COLL")
1847 : :
1848 : 0 : end subroutine gather_matlu
1849 : : !!***
1850 : :
1851 : : !!****f* m_matlu/diag_matlu
1852 : : !! NAME
1853 : : !! diag_matlu
1854 : : !!
1855 : : !! FUNCTION
1856 : : !! Diagonalize hermitian matlu matrix
1857 : : !!
1858 : : !! COPYRIGHT
1859 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
1860 : : !! This file is distributed under the terms of the
1861 : : !! GNU General Public License, see ~abinit/COPYING
1862 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1863 : : !!
1864 : : !! INPUTS
1865 : : !! matlu(natom) :: input quantity to diagonalize (careful, matlu must be hermitian)
1866 : : !! natom=number of atoms
1867 : : !! prtopt: >=3 : print eigenvectors
1868 : : !! >=4 : print matlu before diagonalization
1869 : : !! nsppol_imp= if 1, one can diagonalize with the same matrix the Up
1870 : : !! and Dn matlu matrix. It is convenient because one can thus have the
1871 : : !! same interaction matrix for up and dn spins. Default is nsppol.
1872 : : !! checkstop= if true (default), print the matrix for spin down in the diagonalization basis of spin up
1873 : : !! (useful when nsppol=2 and nsppol_imp=1)
1874 : : !! optreal= diagonalize the real matrix if max(imag(matlu)) < 1e-6
1875 : : !! test= if 8 or 10, use the block diagonalization algorithm (8 for real and 10 for complex)
1876 : : !!
1877 : : !! OUTPUT
1878 : : !! matlu_diag(natom) :: diagonalized density matrix
1879 : : !! eigvectmatlu(natom) = Eigenvectors corresponding to the diagonalization
1880 : : !!
1881 : : !! SIDE EFFECTS
1882 : : !!
1883 : : !! NOTES
1884 : : !!
1885 : : !! SOURCE
1886 : :
1887 : 189 : subroutine diag_matlu(matlu,matlu_diag,natom,prtopt,eigvectmatlu,nsppol_imp,checkstop,opt_real,test)
1888 : :
1889 : : !Arguments ------------------------------------
1890 : : integer, intent(in) :: natom,prtopt
1891 : : type(matlu_type), intent(in) :: matlu(natom)
1892 : : type(matlu_type), intent(inout) :: eigvectmatlu(natom),matlu_diag(natom) !vz_i
1893 : : integer, optional, intent(in) :: nsppol_imp,opt_real,test
1894 : : logical, optional, intent(in) :: checkstop
1895 : : !Local variables-------------------------------
1896 : : integer :: iatom,im1,im2,info,isppol,lpawu,lwork,lworkr
1897 : : integer :: nspinor,nsppol,nsppolimp,optreal,tndim
1898 : : logical :: blockdiag,blockdiagc,checkstop_in,print_temp_mat2
1899 : : character(len=4) :: tag
1900 : : character(len=500) :: message
1901 : 189 : real(dp), allocatable :: eig(:),rwork(:),valuer(:,:),work(:)!,valuer2(:,:)
1902 : : !real(dp),allocatable :: valuer3(:,:),valuer4(:,:)
1903 : : ! real(dp),allocatable :: eigvec(:,:)
1904 : 189 : complex(dp), allocatable :: temp_mat(:,:),zwork(:)
1905 : : !debug complex(dp),allocatable :: temp_mat3(:,:)
1906 : : !************************************************************************
1907 : :
1908 : 189 : blockdiag = .false.
1909 : 189 : blockdiagc = .false.
1910 : 189 : checkstop_in = .true.
1911 : 189 : nspinor = matlu(1)%nspinor
1912 : 189 : nsppol = matlu(1)%nsppol
1913 : 189 : nsppolimp = nsppol
1914 : 189 : optreal = 0
1915 : :
1916 [ + + ]: 189 : if (present(nsppol_imp)) nsppolimp = nsppol_imp
1917 [ + + ]: 189 : if (present(checkstop)) checkstop_in = checkstop
1918 [ + + ]: 189 : if (present(test)) blockdiag = (test == 8 )
1919 : 175 : if (present(test)) blockdiagc = ( test == 10)
1920 [ + + ]: 189 : if (present(opt_real)) optreal = opt_real
1921 : :
1922 : 189 : call zero_matlu(matlu_diag(:),natom)
1923 : 189 : call copy_matlu(matlu(:),eigvectmatlu(:),natom)
1924 : :
1925 : : !donotdiag=.true.
1926 : : !donotdiag=.false.
1927 : : ! ===========================
1928 : : ! Check is diagonalization is necessary and how
1929 : : ! ===========================
1930 : : !do isppol=1,matlu(1)%nsppol
1931 : : ! do iatom=1,natom
1932 : : ! if(matlu(iatom)%lpawu.ne.-1) then
1933 : : ! tndim=(2*matlu(iatom)%lpawu+1)
1934 : : ! do im1=1,tndim
1935 : : ! do im2=1,tndim
1936 : : ! do ispinor=1,nspinor
1937 : : ! do ispinor1=1,nspinor
1938 : : ! if(abs(matlu(iatom)%mat(im1,im2,isppol,ispinor,ispinor1))>tol8.and.&
1939 : : !& (im1/=im2.or.ispinor/=ispinor1)) then
1940 : : ! if matrix is diagonal: do not diagonalize
1941 : : ! donotdiag=.false.
1942 : : ! exit
1943 : : ! endif
1944 : : ! enddo
1945 : : ! enddo
1946 : : ! enddo
1947 : : ! enddo
1948 : : ! endif
1949 : : ! enddo
1950 : : !enddo
1951 : :
1952 : : !if(donotdiag) then
1953 : : ! do isppol=1,matlu(1)%nsppol
1954 : : ! do iatom=1,natom
1955 : : ! if(matlu(iatom)%lpawu.ne.-1) then
1956 : : ! tndim=(2*matlu(iatom)%lpawu+1)
1957 : : ! eigvectmatlu(iatom,isppol)%value(:,:)=czero
1958 : : ! do im1=1,tndim
1959 : : ! eigvectmatlu(iatom,isppol)%value(im1,im1)=cone
1960 : : ! enddo
1961 : : ! endif
1962 : : ! enddo
1963 : : ! enddo
1964 : : ! call copy_matlu(matlu,matlu_diag,natom)
1965 : : ! write(message,'(a)') " Diagonalisation of matlu will not be performed"
1966 : : ! call wrtout(std_out,message,'COLL')
1967 : : ! return
1968 : : !endif
1969 : :
1970 : : ! For nsppol=2, and if nsppolimp=1, the eigenvectors are computed for isppol=1, and applied through
1971 : : ! rotate_matlu to isppol=2. It is the reason why the sum below is only from 1 to nsppolimp !
1972 : :
1973 [ + + ]: 875 : do iatom=1,natom
1974 : :
1975 : 686 : lpawu = matlu(iatom)%lpawu
1976 [ + + ]: 686 : if (lpawu == -1) cycle
1977 : 237 : write(tag,'(i4)') iatom
1978 : 237 : tndim = nspinor * (2*lpawu+1)
1979 : 237 : lwork = 2*tndim - 1
1980 : 237 : lworkr = tndim * (tndim+2) * 2
1981 [ - + + - : 711 : ABI_MALLOC(eig,(tndim))
+ - + - -
+ - + ]
1982 : :
1983 : : ! ===========================
1984 : : ! Define gathermatlu
1985 : : ! ===========================
1986 : : !ABI_MALLOC(gathermatlu,(natom))
1987 : : !do iatom=1,natom
1988 : : ! if(matlu(iatom)%lpawu.ne.-1) then
1989 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
1990 : : ! ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
1991 : : ! gathermatlu(iatom)%value=czero
1992 : : ! endif
1993 : : !enddo
1994 : : !if(nsppol==1.and.nspinor==2) then
1995 : : ! call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
1996 : : !else if((nsppol==2.or.nsppol==1).and.nspinor==1) then
1997 : : ! do iatom=1,natom
1998 : : ! if(matlu(iatom)%lpawu.ne.-1) then
1999 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2000 : : ! do im1=1,tndim
2001 : : ! do im2=1,tndim
2002 : : ! gathermatlu(iatom)%value(im1,im2)=matlu(iatom)%mat(im1,im2,isppol,1,1)
2003 : : ! enddo
2004 : : ! enddo
2005 : : ! endif
2006 : : ! enddo
2007 : : !endif
2008 : :
2009 : : ! ===========================
2010 : : ! Diagonalize
2011 : : ! ===========================
2012 [ + + ]: 599 : do isppol=1,nsppolimp
2013 : : !debug allocate(temp_mat2(tndim,tndim))
2014 : : !debug temp_mat2=zero
2015 : : ! ABI_MALLOC(valuer2,(tndim,tndim))
2016 : : ! ABI_MALLOC(valuer3,(tndim,tndim))
2017 : : ! ABI_MALLOC(valuer4,(tndim,tndim))
2018 : : ! valuer2=zero
2019 : : ! valuer3=zero
2020 : : ! valuer4=zero
2021 [ + + ]: 362 : if (prtopt >= 4) then
2022 : 14 : write(message,'(a,i4,a,i4)') " BEFORE DIAGONALIZATION for atom",iatom," and isppol",isppol
2023 : 14 : call wrtout(std_out,message,'COLL')
2024 [ + + ]: 112 : do im1=1,tndim
2025 [ + - + + ]: 784 : write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))') (matlu(iatom)%mat(im1,im2,isppol),im2=1,tndim)
2026 : 112 : call wrtout(std_out,message,'COLL')
2027 : : end do ! im1
2028 : : end if ! prtopt>=4
2029 : : !debug temp_mat2(:,:)=gathermatlu(iatom)%value(:,:)
2030 : : ! write(std_out,*)"diag"
2031 : :
2032 [ + + + + : 11186 : if (optreal == 1 .and. maxval(abs(aimag(matlu(iatom)%mat(:,:,isppol)))) < tol6 ) then
+ + + - -
+ - - + +
+ + ]
2033 : 142 : write(message,'(a,2x,a,e9.3,a)') ch10,"Imaginary part of Local Hamiltonian is lower than ",&
2034 : 284 : & tol6,": the real matrix is used"
2035 : 142 : call wrtout(std_out,message,'COLL')
2036 [ - + + - : 568 : ABI_MALLOC(valuer,(tndim,tndim))
- + + - +
- + - - +
- + ]
2037 [ + + + + ]: 3786 : valuer(:,:) = dble(matlu(iatom)%mat(:,:,isppol))
2038 : : ! write(message,'(a)') ch10
2039 : : ! call wrtout(std_out,message,'COLL')
2040 : : ! write(message,'(a,i4,a,i4)') "BEFORE valuer for atom",iatom," and isppol",isppol
2041 : : ! call wrtout(std_out,message,'COLL')
2042 : : ! do im1=1,tndim
2043 : : ! write(message,'(2(1x,18(1x,"(",f20.15,",",f20.15,")")))')&
2044 : : !& (valuer(im1,im2),im2=1,tndim)
2045 : : ! call wrtout(std_out,message,'COLL')
2046 : : ! end do
2047 : : ! do im1=1,tndim
2048 : : ! valuer(im1,im1)=real(im1,kind=dp)*0.00000000001_dp+valuer(im1,im1)
2049 : : ! enddo
2050 : : ! write(message,'(a)') ch10
2051 : : ! call wrtout(std_out,message,'COLL')
2052 : : ! write(message,'(a,i4,a,i4)') "BEFORE valuer for atom",iatom," and isppol",isppol
2053 : : ! call wrtout(std_out,message,'COLL')
2054 : : ! do im1=1,tndim
2055 : : ! write(message,'(2(1x,18(1x,f20.15,f20.15)))')&
2056 : : !& (valuer(im1,im2),im2=1,tndim)
2057 : : ! call wrtout(std_out,message,'COLL')
2058 : : ! end do
2059 : : !call dsyev('v','u',tndim,valuer,tndim,eig,work,lworkr,info)
2060 [ - + ]: 142 : if (blockdiag .or. blockdiagc) then
2061 : 0 : call blockdiago_fordsyev(valuer(:,:),tndim,eig(:))
2062 : : else
2063 [ - + + - : 426 : ABI_MALLOC(work,(lworkr))
+ - + - -
+ - + ]
2064 : 142 : call dsyev('v','u',tndim,valuer(:,:),tndim,eig(:),work(:),lworkr,info)
2065 [ - + ]: 142 : ABI_FREE(work)
2066 : : end if ! blockdiag
2067 : : !! For reproductibility
2068 : : ! ! valuer2: eigenvector for the perturb matrix
2069 : : ! valuer2=real(gathermatlu(iatom)%value,kind=dp)
2070 : : ! do im1=1,tndim
2071 : : ! valuer2(im1,im1)=float(im1)*0.00000000001+valuer2(im1,im1)
2072 : : ! enddo
2073 : : ! call dsyev('v','u',tndim,valuer2,tndim,eig,work,lworkr,info)
2074 : : ! write(message,'(a)') ch10
2075 : : ! call wrtout(std_out,message,'COLL')
2076 : : ! write(message,'(a,i4,a,i4)') " valuer2 for atom",iatom," and isppol",isppol
2077 : : ! call wrtout(std_out,message,'COLL')
2078 : : ! do im1=1,tndim
2079 : : ! write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
2080 : : !& (valuer2(im1,im2),im2=1,tndim)
2081 : : ! call wrtout(std_out,message,'COLL')
2082 : : ! end do
2083 : : ! call dgemm('n','n',tndim,tndim,tndim,cone,valuer,tndim,&
2084 : : !& valuer2,tndim,czero,valuer3 ,tndim)
2085 : : ! call dgemm('c','n',tndim,tndim,tndim,cone,valuer2,tndim,&
2086 : : !& valuer3 ,tndim,czero,valuer4,tndim)
2087 : : ! ! valuer4: compute unpert matrix in the basis of the
2088 : : ! ! perturb basis
2089 : : ! write(message,'(a)') ch10
2090 : : ! call wrtout(std_out,message,'COLL')
2091 : : ! write(message,'(a,i4,a,i4)') "BEFORE valuer4 for atom",iatom," and isppol",isppol
2092 : : ! call wrtout(std_out,message,'COLL')
2093 : : ! do im1=1,tndim
2094 : : ! write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
2095 : : !& (valuer4(im1,im2),im2=1,tndim)
2096 : : ! call wrtout(std_out,message,'COLL')
2097 : : ! end do
2098 : : ! call dsyev('v','u',tndim,valuer4,tndim,eig,work,lworkr,info)
2099 : : ! ! valuer4: Diago valuer4 (nearly diag)
2100 : : ! write(message,'(a)') ch10
2101 : : ! call wrtout(std_out,message,'COLL')
2102 : : ! write(message,'(a,i4,a,i4)') "AFTER valuer4 for atom",iatom," and isppol",isppol
2103 : : ! call wrtout(std_out,message,'COLL')
2104 : : ! do im1=1,tndim
2105 : : ! write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
2106 : : !& (valuer4(im1,im2),im2=1,tndim)
2107 : : ! call wrtout(std_out,message,'COLL')
2108 : : ! end do
2109 : : ! call dgemm('n','n',tndim,tndim,tndim,cone,valuer2,tndim,&
2110 : : !& valuer4,tndim,czero,valuer ,tndim)
2111 : : !write(6,*) "INFO",info
2112 [ + + + + ]: 3786 : eigvectmatlu(iatom)%mat(:,:,isppol) = cmplx(valuer(:,:),zero,kind=dp)
2113 [ - + ]: 142 : ABI_FREE(valuer)
2114 : : ! write(message,'(a,i4,a,i4)') "AFTER valuer for atom",iatom," and isppol",isppol
2115 : : ! call wrtout(std_out,message,'COLL')
2116 : : ! do im1=1,tndim
2117 : : ! write(message,'(2(1x,18(1x,"(",f20.15,",",f20.15,")")))')&
2118 : : !& (valuer(im1,im2),im2=1,tndim)
2119 : : ! call wrtout(std_out,message,'COLL')
2120 : : ! end do
2121 : : else
2122 : :
2123 [ - + ]: 220 : if (blockdiagc) then
2124 : 0 : write(message,'(a,a,a)') ch10, " == The local Hamiltonian in Ylm basis is complex.&
2125 : 0 : & The complex matrix is used for the diagonalisation. Printing real and imaginary part of rotation matrix: "
2126 : 0 : call wrtout(std_out,message,'COLL')
2127 : :
2128 [ # # # # ]: 0 : eigvectmatlu(iatom)%mat(:,:,isppol) = matlu(iatom)%mat(:,:,isppol)
2129 : :
2130 : 0 : call blockdiago_forzheev(eigvectmatlu(iatom)%mat(:,:,isppol),tndim,eig(:))
2131 : :
2132 : : !ABI_MALLOC(zwork,(lwork))
2133 : : !ABI_MALLOC(rwork,(3*tndim-2))
2134 : : !call zheev('v','u',tndim,eigvectmatlu(iatom)%mat(:,:,isppol),tndim,eig(:),zwork(:),lwork,rwork(:),info)
2135 : : !ABI_FREE(zwork)
2136 : : !ABI_FREE(rwork)
2137 : : else
2138 [ + + + + : 7400 : if (optreal == 1 .and. maxval(abs(aimag(matlu(iatom)%mat(:,:,isppol)))) > tol8 ) then
+ + + - -
+ - - + +
+ - ]
2139 : 2 : write(message,'(a)') " Local hamiltonian in correlated basis is complex"
2140 : 2 : ABI_COMMENT(message)
2141 : : end if
2142 : : !eigvectmatlu(iatom)%mat(:,:,isppol) = matlu(iatom)%mat(:,:,isppol)
2143 [ - + + - : 660 : ABI_MALLOC(zwork,(lwork))
+ - + - -
+ - + ]
2144 [ + - - + : 660 : ABI_MALLOC(rwork,(3*tndim-2))
+ - + - +
- - + -
+ ]
2145 : 220 : call zheev('v','u',tndim,eigvectmatlu(iatom)%mat(:,:,isppol),tndim,eig(:),zwork(:),lwork,rwork(:),info)
2146 [ - + ]: 220 : ABI_FREE(zwork)
2147 [ - + ]: 220 : ABI_FREE(rwork)
2148 : : endif !blockdiag
2149 : : end if ! present(optreal)
2150 [ + + ]: 362 : if (prtopt >= 3) then
2151 [ + - ]: 56 : write(message,'(3a,i1)') " EIGENVECTORS for atom ",trim(adjustl(tag))," and isppol ",isppol
2152 : 56 : call wrtout(std_out,message,'COLL')
2153 [ + + ]: 302 : do im1=1,tndim
2154 : : !write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))') (eigvectmatlu(iatom)%mat(im1,im2,isppol),im2=1,tndim)
2155 [ + - + + ]: 1684 : write(message,'(12(1x,18(1x,f6.3,1x,f6.3)))') (eigvectmatlu(iatom)%mat(im1,im2,isppol),im2=1,tndim)
2156 : 302 : call wrtout(std_out,message,'COLL')
2157 : : end do ! im1
2158 : : ! do im1=1,tndim
2159 : : ! xcheck=czero
2160 : : ! do im3=1,tndim
2161 : : ! do im2=1,tndim
2162 : : ! xcheck=xcheck+gathermatlu(iatom)%value(im1,im2)*conjg(gathermatlu(iatom)%value(im2,im3))
2163 : : ! end do
2164 : : ! end do
2165 : : ! write(6,*) "check",im3,im1,xcheck
2166 : : ! end do
2167 : : end if ! prtopt>=3
2168 : : ! write(std_out,*) "eig",eig
2169 : : ! ===========================
2170 : : ! Put eigenvalue in matlu_diag
2171 : : ! ===========================
2172 [ + + ]: 2092 : do im1=1,tndim
2173 : 2092 : matlu_diag(iatom)%mat(im1,im1,isppol) = cmplx(eig(im1),zero,kind=dp)
2174 : : end do ! im1
2175 : : !if(prtopt>=2) then
2176 : : ! write(message,'(a,12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
2177 : : !& ch10,(eig(im1),im1=1,tndim)
2178 : : ! call wrtout(std_out,message,'COLL')
2179 : : !call wrtout(std_out,message,'COLL')
2180 : : !write(std_out,*) "EIG", eig
2181 : : !endif
2182 : : ! ABI_FREE(valuer2)
2183 : : ! ABI_FREE(valuer3)
2184 : : ! ABI_FREE(valuer4)
2185 : : ! endif
2186 : : ! enddo
2187 : : ! ===========================
2188 : : ! Keep eigenvectors gathermatlu
2189 : : ! ===========================
2190 : : !if (present(eigvectmatlu)) then
2191 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2192 : : ! eigvectmatlu(iatom,isppol)%value(:,:)=gathermatlu(iatom)%value(:,:)
2193 : : ! write(std_out,*) "eigvect in diag_matlu"
2194 : : ! do im1=1,tndim
2195 : : ! write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
2196 : : !& (gathermatlu(iatom)%value(im1,im2),im2=1,tndim)
2197 : : ! call wrtout(std_out,message,'COLL')
2198 : : ! end do
2199 : : !endif
2200 : :
2201 [ - + ]: 599 : if (nsppolimp == 1 .and. nsppol == 2) then
2202 : : ! ==================================================================
2203 : : ! If necessary rotate levels for this other spin, assuming the same
2204 : : ! rotation matrix: it has to be checked afterwards that the matrix is
2205 : : ! diagonal
2206 : : ! ===================================================================
2207 : :
2208 : : ! input matrix: gathermatlu
2209 : : ! rotation matrix: eigvectmatlu
2210 : : ! intermediate matrix: temp_mat
2211 : : ! result matrix: temp_mat2
2212 : : !do im1=1,tndim
2213 : : ! do im2=1,tndim
2214 : : ! gathermatlu(iatom)%value(im1,im2)=matlu(iatom)%mat(im1,im2,2,1,1)
2215 : : ! enddo
2216 : : !enddo
2217 : :
2218 [ # # ]: 0 : if (prtopt >= 3) then
2219 : 0 : write(message,'(a,i4,a,i4)') " MATLU for atom",iatom," inside if nsppolimp==1"
2220 : 0 : call wrtout(std_out,message,'COLL')
2221 [ # # ]: 0 : do im1=1,tndim
2222 [ # # # # ]: 0 : write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))') (matlu(iatom)%mat(im1,im2,2),im2=1,tndim)
2223 : 0 : call wrtout(std_out,message,'COLL')
2224 : : end do ! im1
2225 : : end if ! prtopt>=3
2226 : :
2227 [ # # # # : 0 : ABI_MALLOC(temp_mat,(tndim,tndim))
# # # # #
# # # # #
# # ]
2228 : :
2229 : : call abi_xgemm('n','n',tndim,tndim,tndim,cone,matlu(iatom)%mat(:,:,2),tndim,&
2230 : 0 : & eigvectmatlu(iatom)%mat(:,:,1),tndim,czero,temp_mat(:,:),tndim)
2231 : :
2232 : : call abi_xgemm('c','n',tndim,tndim,tndim,cone,eigvectmatlu(iatom)%mat(:,:,1),tndim,&
2233 : 0 : & temp_mat(:,:),tndim,czero,matlu_diag(iatom)%mat(:,:,2),tndim)
2234 : :
2235 [ # # ]: 0 : ABI_FREE(temp_mat)
2236 : :
2237 [ # # # # ]: 0 : eigvectmatlu(iatom)%mat(:,:,2) = eigvectmatlu(iatom)%mat(:,:,1)
2238 : : print_temp_mat2 = .false.
2239 [ # # ]: 0 : do im2=1,tndim
2240 [ # # ]: 0 : do im1=1,tndim
2241 [ # # # # ]: 0 : if (im1 /= im2 .and. abs(matlu_diag(iatom)%mat(im1,im2,2)) > tol5) then
2242 : 0 : write(message,'(3a,i4,2f16.4)') ch10,'diag_matlu= Matrix for spin number 2 obtained with', &
2243 : 0 : & ' eigenvectors from diagonalization for spin nb 1 is non diagonal for atom:',iatom,&
2244 : 0 : & abs(matlu_diag(iatom)%mat(im1,im2,2)),tol5
2245 : 0 : call wrtout(std_out,message,'COLL')
2246 [ # # # # ]: 0 : if (abs(matlu_diag(iatom)%mat(im1,im2,2)) > tol1 .or. checkstop_in) print_temp_mat2 = .true.
2247 : : end if
2248 : : end do ! im1
2249 : : end do ! im2
2250 : :
2251 [ # # # # ]: 0 : if (print_temp_mat2 .and. prtopt >= 3) then
2252 : 0 : write(message,'(a)') " temp_mat2"
2253 : 0 : call wrtout(std_out,message,'COLL')
2254 [ # # ]: 0 : do im1=1,tndim
2255 [ # # # # ]: 0 : write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))') (matlu_diag(iatom)%mat(im1,im2,2),im2=1,tndim)
2256 : 0 : call wrtout(std_out,message,'COLL')
2257 : : end do ! im1
2258 [ # # ]: 0 : if (iatom == 2) ABI_ERROR("iatom==2")
2259 : : end if ! print_temp_mat2
2260 : : end if ! nsppol_imp=1 and nsppol=2
2261 : : end do ! isppol
2262 : :
2263 [ - + ]: 875 : ABI_FREE(eig)
2264 : :
2265 : : !! for check only
2266 : : !debug if(matlu(iatom)%lpawu.ne.-1) then
2267 : : !debug allocate(temp_mat(tndim,tndim))
2268 : : !debug allocate(temp_mat3(tndim,tndim))
2269 : : !debug do im1=1,tndim
2270 : : !debug do im2=1,tndim
2271 : : !debug! rot_mat(iatom,isppol)%value(im1,im2)=rot_mat_orig(iatom,isppol)%value(im1,im2)
2272 : : !debug temp_mat3(im1,im2)=conjg(gathermatlu(iatom)%value(im2,im1))
2273 : : !debug enddo
2274 : : !debug enddo
2275 : : !debug temp_mat(:,:)=czero
2276 : : !debug! input matrix: temp_mat2
2277 : : !debug! rotation matrix: gathermatlu
2278 : : !debug! intermediate matrix: temp_mat
2279 : : !debug! result matrix: temp_mat2
2280 : : !debug call zgemm('n','c',tndim,tndim,tndim,cone,temp_mat2 ,tndim,&
2281 : : !debug& temp_mat3,tndim,czero,temp_mat ,tndim)
2282 : : !debug call zgemm('n','n',tndim,tndim,tndim,cone,temp_mat3,tndim,&
2283 : : !debug& temp_mat ,tndim,czero,temp_mat2,tndim)
2284 : : !debug! call zgemm('n','c',tndim,tndim,tndim,cone,temp_mat2 ,tndim,&
2285 : : !debug!& gathermatlu(iatom)%value,tndim,czero,temp_mat ,tndim)
2286 : : !debug! call zgemm('n','n',tndim,tndim,tndim,cone,gathermatlu(iatom)%value,tndim,&
2287 : : !debug!& temp_mat ,tndim,czero,temp_mat2,tndim)
2288 : : !debug write(std_out,*) "result"
2289 : : !debug do im1=1,tndim
2290 : : !debug write(message,'(12(1x,18(1x,"(",f9.3,",",f9.3,")")))')&
2291 : : !debug& (temp_mat2(im1,im2),im2=1,tndim)
2292 : : !debug call wrtout(std_out,message,'COLL')
2293 : : !debug end do
2294 : : !debug deallocate(temp_mat)
2295 : : !debug deallocate(temp_mat3)
2296 : : !debug endif ! lpawu
2297 : :
2298 : : end do ! iatom
2299 : : ! End loop over atoms
2300 : : ! ===========================
2301 : :
2302 : 0 : end subroutine diag_matlu
2303 : : !!***
2304 : :
2305 : : !!****f* m_matlu/rotate_matlu
2306 : : !! NAME
2307 : : !! rotate_matlu
2308 : : !!
2309 : : !! FUNCTION
2310 : : !! Rotate matlu matrix
2311 : : !!
2312 : : !! COPYRIGHT
2313 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2314 : : !! This file is distributed under the terms of the
2315 : : !! GNU General Public License, see ~abinit/COPYING
2316 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2317 : : !!
2318 : : !! INPUTS
2319 : : !! matlu_inp(natom) :: input quantity to rotate
2320 : : !! rot_mat(natom) :: Rotation matrix (usually from diag_matlu)
2321 : : !! natom=number of atoms in cell.
2322 : : !! inverse= 1: rot_mat^H * matlu * rot_mat (from original basis to diagonal basis)
2323 : : !! /=1: rot_mat * matlu * rot_mat^H (from diagonal basis to original basis)
2324 : : !!
2325 : : !! OUTPUT
2326 : : !!
2327 : : !! SIDE EFFECTS
2328 : : !!
2329 : : !! NOTES
2330 : : !!
2331 : : !! SOURCE
2332 : :
2333 : 21879 : subroutine rotate_matlu(matlu_inp,rot_mat,natom,inverse)
2334 : :
2335 : : !Arguments ------------------------------------
2336 : : integer, intent(in) :: inverse,natom
2337 : : type(matlu_type), intent(inout) :: matlu_inp(natom)
2338 : : type(matlu_type), intent(in) :: rot_mat(natom)
2339 : : !Local variables-------------------------------
2340 : : integer :: iatom,isppol,lpawu,nspinor,nsppol,tndim
2341 : 21879 : complex(dp), allocatable :: temp_mat(:,:)
2342 : : character(len=1) :: c1,c2
2343 : : !************************************************************************
2344 : :
2345 : 21879 : nspinor = matlu_inp(1)%nspinor
2346 : 21879 : nsppol = matlu_inp(1)%nsppol
2347 : :
2348 [ + + ]: 21879 : if (inverse == 1) then
2349 : 10581 : c1 = "n" ; c2 = "c"
2350 : : else
2351 : 11298 : c1 = "c" ; c2 = "n"
2352 : : end if ! inverse
2353 : :
2354 [ + + ]: 69726 : do iatom=1,natom
2355 : 47847 : lpawu = matlu_inp(iatom)%lpawu
2356 [ + + ]: 47847 : if (lpawu == -1) cycle
2357 : 23031 : tndim = nspinor * (2*lpawu+1)
2358 [ - + + - : 92124 : ABI_MALLOC(temp_mat,(tndim,tndim))
- + + - +
- + - - +
+ - ]
2359 [ + + ]: 54190 : do isppol=1,nsppol
2360 : : call abi_xgemm('n',c1,tndim,tndim,tndim,cone,matlu_inp(iatom)%mat(:,:,isppol),tndim,&
2361 : 31159 : & rot_mat(iatom)%mat(:,:,isppol),tndim,czero,temp_mat(:,:),tndim)
2362 : : call abi_xgemm(c2,'n',tndim,tndim,tndim,cone,rot_mat(iatom)%mat(:,:,isppol),tndim,&
2363 : 54190 : & temp_mat(:,:),tndim,czero,matlu_inp(iatom)%mat(:,:,isppol),tndim)
2364 : : end do ! isppol
2365 [ - + ]: 69726 : ABI_FREE(temp_mat)
2366 : : end do ! iatom
2367 : :
2368 : : !do isppol=1,nsppol
2369 : :
2370 : : ! ===========================
2371 : : ! Define gathermatlu and rot_mat_orig and allocate
2372 : : ! ===========================
2373 : : !ABI_MALLOC(rot_mat_orig,(natom))
2374 : : !ABI_MALLOC(gathermatlu,(natom))
2375 : : !do iatom=1,natom
2376 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2377 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2378 : : ! ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
2379 : : ! gathermatlu(iatom)%value=czero
2380 : : ! ABI_MALLOC(rot_mat_orig(iatom,isppol)%value,(tndim,tndim))
2381 : : ! rot_mat_orig(iatom,isppol)%value(:,:)=rot_mat(iatom,isppol)%value(:,:)
2382 : : ! ABI_MALLOC(rot_mat_orig(iatom)%value,(tndim,tndim))
2383 : : ! rot_mat_orig(iatom)%value(:,:)=rot_mat(iatom,isppol)%value(:,:)
2384 : : ! endif
2385 : : !enddo
2386 : : !if(nsppol==1.and.nspinor==2) then
2387 : : ! call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
2388 : : !else if((nsppol==2.or.nsppol==1).and.nspinor==1) then
2389 : : ! do iatom=1,natom
2390 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2391 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2392 : : ! do im1=1,tndim
2393 : : ! do im2=1,tndim
2394 : : ! gathermatlu(iatom)%value(im1,im2)=matlu(iatom)%mat(im1,im2,isppol,1,1)
2395 : : ! enddo
2396 : : ! enddo
2397 : : ! endif
2398 : : ! enddo
2399 : : !endif
2400 : : ! write(std_out,*) "gathermatlu in rotate matlu"
2401 : : ! do im1=1,tndim
2402 : : ! write(message,'(12(1x,18(1x,"(",e17.10,",",e17.10,")")))')&
2403 : : ! (gathermatlu(1)%value(im1,im2),im2=1,tndim)
2404 : : ! call wrtout(std_out,message,'COLL')
2405 : : ! end do
2406 : :
2407 : : ! ===========================
2408 : : ! If necessary, invert rot_mat
2409 : : ! ===========================
2410 : : !if(inverse==1) then
2411 : : ! do iatom=1,natom
2412 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2413 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2414 : : ! do im1=1,tndim
2415 : : ! do im2=1,tndim
2416 : : ! rot_mat(iatom,isppol)%value(im1,im2)=conjg(rot_mat_orig(iatom,isppol)%value(im2,im1))
2417 : : ! rot_mat(iatom,isppol)%value(im1,im2)=conjg(rot_mat_orig(iatom)%value(im2,im1))
2418 : : ! enddo
2419 : : ! enddo
2420 : : ! endif ! lpawu
2421 : : ! enddo ! iatom
2422 : : !endif
2423 : : ! write(std_out,*) "rot_mat_orig "
2424 : : ! do im1=1,tndim
2425 : : ! write(message,'(12(1x,18(1x,"(",e18.10,",",e18.10,")")))')&
2426 : : ! & (rot_mat_orig(1)%value(im1,im2),im2=1,tndim)
2427 : : ! call wrtout(std_out,message,'COLL')
2428 : : ! end do
2429 : : ! write(std_out,*) "rot_mat "
2430 : : ! do im1=1,tndim
2431 : : ! write(message,'(12(1x,18(1x,"(",e18.10,",",e18.10,")")))')&
2432 : : ! & (rot_mat(1,1)%value(im1,im2),im2=1,tndim)
2433 : : ! call wrtout(std_out,message,'COLL')
2434 : : ! end do
2435 : :
2436 : : ! ===========================
2437 : : ! Rotate
2438 : : ! ===========================
2439 : : !ABI_MALLOC(temp_mat,(tndim,tndim))
2440 : : !do iatom=1,natom
2441 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2442 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2443 : : ! temp_mat(:,:)=czero
2444 : : ! input matrix: gathermatlu
2445 : : ! rotation matrix: rot_mat
2446 : : ! intermediate matrix: temp_mat
2447 : : ! result matrix: gathermatlu
2448 : : ! temp_mat = gathermatlu * conjg(rot_mat)
2449 : : ! call zgemm('n','c',tndim,tndim,tndim,cone,gathermatlu(iatom)%value
2450 : : ! ,tndim,&
2451 : : !& rot_mat(iatom,isppol)%value,tndim,czero,temp_mat
2452 : : !,tndim)
2453 : : ! gathermatlu = rot_mat * temp_mat = rot_mat * gathermatlu *
2454 : : ! conjg(rot_mat)
2455 : : ! call
2456 : : ! zgemm('n','n',tndim,tndim,tndim,cone,rot_mat(iatom,isppol)%value,tndim,&
2457 : : !& temp_mat
2458 : : !,tndim,czero,gathermatlu(iatom)%value,tndim)
2459 : : ! endif ! lpawu
2460 : : ! enddo ! iatom
2461 : : !do iatom=1,natom
2462 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2463 : : ! write(std_out,*) "temp_mat in rotate_matlu 2"
2464 : : ! do im1=1,tndim
2465 : : ! write(message,'(12(1x,18(1x,"(",f17.10,",",f17.10,")")))')&
2466 : : !& (temp_mat(im1,im2),im2=1,tndim)
2467 : : ! call wrtout(std_out,message,'COLL')
2468 : : ! end do
2469 : : ! endif
2470 : : !enddo
2471 : : !do iatom=1,natom
2472 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2473 : : ! write(std_out,*) "gathermatlu in rotate_matlu 2"
2474 : : ! do im1=1,tndim
2475 : : ! write(message,'(12(1x,18(1x,"(",f17.10,",",f17.10,")")))')&
2476 : : !& (gathermatlu(iatom)%value(im1,im2),im2=1,tndim)
2477 : : ! call wrtout(std_out,message,'COLL')
2478 : : ! end do
2479 : : ! endif
2480 : : !enddo
2481 : : ! ABI_FREE(temp_mat)
2482 : : !ABI_ERROR("Aborting now")
2483 : :
2484 : : ! Choose inverse rotation: reconstruct correct rot_mat from rot_mat_orig
2485 : : ! ========================================================================
2486 : : !if(inverse==1) then
2487 : : ! do iatom=1,natom
2488 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2489 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2490 : : ! do im1=1,tndim
2491 : : ! do im2=1,tndim
2492 : : ! !
2493 : : ! rot_mat(iatom,isppol)%value(im1,im2)=rot_mat_orig(iatom,isppol)%value(im1,im2)
2494 : : ! rot_mat(iatom,isppol)%value(im1,im2)=rot_mat_orig(iatom)%value(im1,im2)
2495 : : ! enddo
2496 : : ! enddo
2497 : : ! endif ! lpawu
2498 : : ! enddo ! iatom
2499 : : !endif
2500 : :
2501 : : ! ===========================
2502 : : ! Put data into matlu(iatom)
2503 : : ! ===========================
2504 : : !if(nsppol==1.and.nspinor==2) then
2505 : : ! call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
2506 : : !else if((nsppol==2.or.nsppol==1).and.nspinor==1) then
2507 : : ! do iatom=1,natom
2508 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2509 : : ! tndim=nspinor*(2*matlu(iatom)%lpawu+1)
2510 : : ! do im1=1,tndim
2511 : : ! do im2=1,tndim
2512 : : ! matlu(iatom)%mat(im1,im2,isppol,1,1)=
2513 : : ! gathermatlu(iatom)%value(im1,im2)
2514 : : ! enddo
2515 : : ! enddo
2516 : : ! endif
2517 : : ! enddo
2518 : : !endif ! test nsppol/nspinor
2519 : : ! ===========================
2520 : : ! Deallocations
2521 : : ! ===========================
2522 : : !do iatom=1,natom
2523 : : ! if(matlu(iatom)%lpawu.ne.-1) then
2524 : : ! ABI_FREE(gathermatlu(iatom)%value)
2525 : : ! ABI_FREE(rot_mat_orig(iatom,isppol)%value)
2526 : : ! ABI_FREE(rot_mat_orig(iatom)%value)
2527 : : ! endif
2528 : : !enddo
2529 : : !ABI_FREE(gathermatlu)
2530 : : !ABI_FREE(rot_mat_orig)
2531 : : !enddo ! isppol
2532 : :
2533 : 21879 : end subroutine rotate_matlu
2534 : : !!***
2535 : :
2536 : : !!****f* m_matlu/shift_matlu
2537 : : !! NAME
2538 : : !! shift_matlu
2539 : : !!
2540 : : !! FUNCTION
2541 : : !! Add/subtract a scalar to the diagonal part
2542 : : !!
2543 : : !! COPYRIGHT
2544 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2545 : : !! This file is distributed under the terms of the
2546 : : !! GNU General Public License, see ~abinit/COPYING
2547 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2548 : : !!
2549 : : !! INPUTS
2550 : : !! matlu(natom) :: input quantity to shift
2551 : : !! natom=number of atoms in cell.
2552 : : !! shift= shift of the diagonal part.
2553 : : !! signe= 1 (default) : the shift is substracted
2554 : : !! = -1 : the shift is added
2555 : : !!
2556 : : !! OUTPUT
2557 : : !! matlu(natom) :: shifted matrix
2558 : : !!
2559 : : !! SIDE EFFECTS
2560 : : !!
2561 : : !! NOTES
2562 : : !!
2563 : : !! SOURCE
2564 : :
2565 : 18679 : subroutine shift_matlu(matlu,natom,shift,signe)
2566 : :
2567 : : !Arguments ------------------------------------
2568 : : integer, intent(in) :: natom
2569 : : type(matlu_type), intent(inout) :: matlu(natom)
2570 : : complex(dp), intent(in) :: shift(natom)
2571 : : integer, optional, intent(in) :: signe
2572 : : !Local variables-------------------------------
2573 : : integer :: iatom,im,lpawu,ndim,nspinor,nsppol,signe_used
2574 : : ! character(len=500) :: message
2575 : : !************************************************************************
2576 : :
2577 : 18679 : nspinor = matlu(1)%nspinor
2578 : 18679 : nsppol = matlu(1)%nsppol
2579 : 18679 : signe_used = 1
2580 : :
2581 [ + + ]: 18679 : if (present(signe)) then
2582 [ - + ]: 682 : if (signe == -1) signe_used = -1
2583 : : end if ! present(signe)
2584 : :
2585 [ + + ]: 41462 : do iatom=1,natom
2586 : 22783 : lpawu = matlu(iatom)%lpawu
2587 [ + + ]: 22783 : if (lpawu == -1) cycle
2588 : 18799 : ndim = nspinor * (2*lpawu+1)
2589 [ + + ]: 132501 : do im=1,ndim
2590 [ + + + + ]: 300836 : matlu(iatom)%mat(im,im,:) = matlu(iatom)%mat(im,im,:) + merge(-shift(iatom),shift(iatom),signe_used==1)
2591 : : end do ! im
2592 : : end do ! iatom
2593 : :
2594 : 18679 : end subroutine shift_matlu
2595 : : !!***
2596 : :
2597 : : !!****f* m_matlu/checkreal_matlu
2598 : : !! NAME
2599 : : !! checkreal_matlu
2600 : : !!
2601 : : !! FUNCTION
2602 : : !! Check that matlu is real and diagonal with given precision
2603 : : !!
2604 : : !! COPYRIGHT
2605 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2606 : : !! This file is distributed under the terms of the
2607 : : !! GNU General Public License, see ~abinit/COPYING
2608 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2609 : : !!
2610 : : !! INPUTS
2611 : : !! dmft_solv :: impurity solver
2612 : : !! matlu(natom) :: input quantity to check
2613 : : !! natom=number of atoms in cell.
2614 : : !! tol : threshold. Print a warning if max(abs(imag(off diagonal elements))) > tol or
2615 : : !! max(abs(off diagonal elements)) > tol, and throws an error if max(abs(imag(diagonal elements))) > tol
2616 : : !!
2617 : : !! SIDE EFFECTS
2618 : : !!
2619 : : !! NOTES
2620 : : !!
2621 : : !! SOURCE
2622 : :
2623 : 0 : subroutine checkreal_matlu(matlu,natom,tol)
2624 : :
2625 : : !Arguments ------------------------------------
2626 : : real(dp), intent(in) :: tol
2627 : : integer, intent(in) :: natom
2628 : : type(matlu_type), intent(in) :: matlu(natom)
2629 : : !Local variables-------------------------------
2630 : : integer :: iatom,im,im1,isppol,lpawu,ndim,nspinor,nsppol
2631 : : character(len=500) :: message
2632 : : real(dp) :: elem,maximag,maximagdiag,maxoffdiag
2633 : : !************************************************************************
2634 : :
2635 : 0 : maximag = zero
2636 : 0 : maximagdiag = zero
2637 : 0 : maxoffdiag = zero
2638 : 0 : nspinor = matlu(1)%nspinor
2639 : 0 : nsppol = matlu(1)%nsppol
2640 : :
2641 [ # # ]: 0 : do iatom=1,natom
2642 : 0 : lpawu = matlu(iatom)%lpawu
2643 [ # # ]: 0 : if (lpawu == -1) cycle
2644 : 0 : ndim = nspinor * (2*lpawu+1)
2645 [ # # ]: 0 : do isppol=1,nsppol
2646 [ # # ]: 0 : do im1=1,ndim
2647 [ # # ]: 0 : do im=1,ndim
2648 : 0 : elem = abs(aimag(matlu(iatom)%mat(im,im1,isppol)))
2649 [ # # ]: 0 : if (elem > maximag) maximag = elem
2650 [ # # ]: 0 : if (im == im1) then
2651 [ # # ]: 0 : if (elem > maximagdiag) maximagdiag = elem
2652 : : else
2653 : 0 : elem = abs(matlu(iatom)%mat(im,im1,isppol))
2654 [ # # ]: 0 : if (elem > maxoffdiag) maxoffdiag = elem
2655 : : end if ! im/=im1
2656 : : end do ! im
2657 : : end do ! im1
2658 : : end do ! isppol
2659 : : end do ! iatom
2660 : :
2661 [ # # ]: 0 : if (maximagdiag > tol) then
2662 : 0 : write(message,'(3x,2a,e12.4,a,e12.4,2a)') ch10,&
2663 : 0 : & ' Diagonal part of the occupation matrix is complex: the imaginary part ',&
2664 : 0 : & maximagdiag,' is larger than',tol,ch10, &
2665 : 0 : & "The calculation cannot handle it : check that your calculation is meaningful"
2666 : 0 : ABI_ERROR(message)
2667 : : end if ! maximagdiag > tol
2668 [ # # ]: 0 : if (maximag > tol) then
2669 : 0 : write(message,'(3x,2a,e12.4,a,e12.4,2a)') ch10,&
2670 : 0 : & ' The off diagonal occupation matrix is complex: the imaginary part ',maximag,' is larger than',tol,ch10,&
2671 : 0 : & "Check that your calculation is meaningful"
2672 : 0 : ABI_WARNING(message)
2673 : : end if ! maximag > tol
2674 [ # # ]: 0 : if (maxoffdiag > tol) then
2675 : 0 : write(message,'(3x,2a,e12.4,a,e12.4,6a)') ch10,&
2676 : 0 : & ' Occupation matrix is non diagonal : the maximum off-diag part ',maxoffdiag,' is larger than',tol,ch10,&
2677 : 0 : & "The corresponding non diagonal elements will be neglected in the Weiss/Hybridization functions",ch10,&
2678 : 0 : & "(Except if dmft_solv=8,9 where these elements are taken into account)",ch10,"This is an approximation."
2679 : 0 : ABI_WARNING(message)
2680 : : else
2681 : 0 : write(message,'(3x,2a,e12.4,a,e12.4,2a)') ch10,' Occupation matrix is diagonal : the off-diag part ',&
2682 : 0 : & maxoffdiag,' is lower than',tol
2683 : 0 : ABI_COMMENT(message)
2684 : : end if ! maxoffdiag > tol
2685 : :
2686 : 0 : end subroutine checkreal_matlu
2687 : : !!***
2688 : :
2689 : : !!****f* m_matlu/checkdiag_matlu
2690 : : !! NAME
2691 : : !! checkdiag_matlu
2692 : : !!
2693 : : !! FUNCTION
2694 : : !! Check that matlu is diagonal in the orbital index with given precision
2695 : : !!
2696 : : !! COPYRIGHT
2697 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2698 : : !! This file is distributed under the terms of the
2699 : : !! GNU General Public License, see ~abinit/COPYING
2700 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2701 : : !!
2702 : : !! INPUTS
2703 : : !! matlu(natom) :: input quantity to check
2704 : : !! natom=number of atoms in cell.
2705 : : !! tol : precision
2706 : : !!
2707 : : !! OUTPUT
2708 : : !! nondiag= true if max(abs(off diagonal elements)) > tol
2709 : : !!
2710 : : !! SIDE EFFECTS
2711 : : !!
2712 : : !! NOTES
2713 : : !!
2714 : : !! SOURCE
2715 : :
2716 : 140 : subroutine checkdiag_matlu(matlu,natom,tol,nondiag)
2717 : :
2718 : : !Arguments ------------------------------------
2719 : : real(dp), intent(in) :: tol
2720 : : integer, intent(in) :: natom
2721 : : logical, intent(out) :: nondiag
2722 : : type(matlu_type), intent(in) :: matlu(natom)
2723 : : !Local variables-------------------------------
2724 : : integer :: iatom,im,im1,isppol,lpawu,ndim,nsppol,nspinor
2725 : : !************************************************************************
2726 : :
2727 : 140 : nondiag = .false.
2728 : 140 : nspinor = matlu(1)%nspinor
2729 : 140 : nsppol = matlu(1)%nsppol
2730 : :
2731 [ + + ]: 572 : do iatom=1,natom
2732 : 432 : lpawu = matlu(iatom)%lpawu
2733 [ + + ]: 432 : if (lpawu == -1) cycle
2734 : 164 : ndim = nspinor * (2*lpawu+1)
2735 [ + + ]: 562 : do isppol=1,nsppol
2736 [ + + ]: 1926 : do im1=1,ndim
2737 [ + + ]: 7966 : do im=1,ndim
2738 : : ! if(im/=im1) write(std_out,*) "im,im1",im,im1,matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor)
2739 : : ! if(present(nondiag).eqv..false.) then
2740 : : ! if(im/=im1.and.(abs(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor))>tol)) then
2741 : : ! write(message,'(5i5)') im,im1,isppol,ispinor,ispinor
2742 : : ! call wrtout(std_out,message,'COLL')
2743 : : ! write(message,'(a,3e16.5)')" checkdiag_matlu: Warning ",matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor),tol
2744 : : ! call wrtout(std_out,message,'COLL')
2745 : : ! if(.not.present(opt)) ABI_ERROR("not present(opt)")
2746 : : ! if(matlu(1)%nspinor==1) ABI_ERROR("matlu%nspinor==1")
2747 : : ! endif
2748 : : ! endif
2749 : : ! if(present(nondiag)) then
2750 [ + + + + ]: 7708 : if ((im /= im1) .and. abs(dble(matlu(iatom)%mat(im,im1,isppol))) > tol) nondiag = .true.
2751 : : ! write(6,*) "NONDIAG", matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)
2752 : : !if(ispinor/=ispinor1.and.(abs(matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))>tol)) then
2753 : : ! write(message,'(a,3e16.5)')" checkdiag_matlu :i Warning ",matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1),tol
2754 : : ! call wrtout(std_out,message,'COLL')
2755 : : ! write(message,'(5i5)') im,im1,isppol,ispinor,ispinor
2756 : : ! call wrtout(std_out,message,'COLL')
2757 : : ! if(matlu(1)%nspinor==1) ABI_ERROR("matlu%nspinor==1")
2758 : : !endif
2759 : : end do ! im
2760 : : end do ! im1
2761 : : end do ! isppol
2762 : : end do ! iatom
2763 : :
2764 : 140 : end subroutine checkdiag_matlu
2765 : : !!***
2766 : :
2767 : : !!****f* m_matlu/prod_matlu
2768 : : !! NAME
2769 : : !! prod_matlu
2770 : : !!
2771 : : !! FUNCTION
2772 : : !! Matrix product of two matlus
2773 : : !!
2774 : : !! COPYRIGHT
2775 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2776 : : !! This file is distributed under the terms of the
2777 : : !! GNU General Public License, see ~abinit/COPYING
2778 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2779 : : !!
2780 : : !! INPUTS
2781 : : !! matlu1(natom) :: input quantity
2782 : : !! matlu2(natom) :: input quantity
2783 : : !!
2784 : : !! OUTPUT
2785 : : !! matlu3(natom) :: output quantity
2786 : : !!
2787 : : !! SIDE EFFECTS
2788 : : !!
2789 : : !! NOTES
2790 : : !!
2791 : : !! SOURCE
2792 : :
2793 : 68 : subroutine prod_matlu(matlu1,matlu2,matlu3,natom)
2794 : :
2795 : : !Arguments ------------------------------------
2796 : : integer, intent(in) :: natom
2797 : : type(matlu_type), intent(in) :: matlu1(natom),matlu2(natom)
2798 : : type(matlu_type), intent(inout) :: matlu3(natom)
2799 : : !Local variables-------------------------------
2800 : : integer :: iatom,isppol,lpawu,ndim,nspinor,nsppol
2801 : : !************************************************************************
2802 : :
2803 : 68 : nspinor = matlu1(1)%nspinor
2804 : 68 : nsppol = matlu1(1)%nsppol
2805 : :
2806 [ + + ]: 192 : do iatom=1,natom
2807 : 124 : lpawu = matlu1(iatom)%lpawu
2808 [ + + ]: 124 : if (lpawu == -1) cycle
2809 : 68 : ndim = nspinor * (2*lpawu+1)
2810 [ + + ]: 244 : do isppol=1,nsppol
2811 : : call abi_xgemm("n","n",ndim,ndim,ndim,cone,matlu1(iatom)%mat(:,:,isppol),ndim,&
2812 : 232 : & matlu2(iatom)%mat(:,:,isppol),ndim,czero,matlu3(iatom)%mat(:,:,isppol),ndim)
2813 : : end do ! isppol
2814 : : end do ! iatom
2815 : :
2816 : : !call zero_matlu(matlu3,natom)
2817 : : !do iatom=1,natom
2818 : : ! lpawu=matlu1(iatom)%lpawu
2819 : : ! if(lpawu.ne.-1) then
2820 : : ! do isppol=1,matlu1(1)%nsppol
2821 : : ! do ispinor1=1,matlu1(1)%nspinor
2822 : : ! do ispinor2=1,matlu1(1)%nspinor
2823 : : ! do ispinor3=1,matlu1(1)%nspinor
2824 : : ! do im1=1,2*lpawu+1
2825 : : ! do im2=1,2*lpawu+1
2826 : : ! do im3=1,2*lpawu+1
2827 : : ! matlu3(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)= &
2828 : : !& matlu3(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)+ &
2829 : : !& matlu1(iatom)%mat(im1,im3,isppol,ispinor1,ispinor3)*&
2830 : : !& matlu2(iatom)%mat(im3,im2,isppol,ispinor3,ispinor2)
2831 : : ! enddo ! im3
2832 : : ! enddo ! im2
2833 : : ! enddo ! im1
2834 : : ! enddo ! ispinor3
2835 : : ! enddo ! ispinor2
2836 : : ! enddo ! ispinor1
2837 : : ! enddo ! isppol
2838 : : ! endif ! lpawu
2839 : : !enddo ! iatom
2840 : :
2841 : 68 : end subroutine prod_matlu
2842 : : !!***
2843 : :
2844 : : !!****f* m_matlu/conjg_matlu
2845 : : !! NAME
2846 : : !! conjg_matlu
2847 : : !!
2848 : : !! FUNCTION
2849 : : !! conjugate of input matlu
2850 : : !!
2851 : : !! COPYRIGHT
2852 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2853 : : !! This file is distributed under the terms of the
2854 : : !! GNU General Public License, see ~abinit/COPYING
2855 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2856 : : !!
2857 : : !! INPUTS
2858 : : !! matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity
2859 : : !!
2860 : : !! SIDE EFFECTS
2861 : : !!
2862 : : !! NOTES
2863 : : !!
2864 : : !! SOURCE
2865 : 0 : subroutine conjg_matlu(matlu1,natom)
2866 : : use defs_wvltypes
2867 : :
2868 : : !Arguments ------------------------------------
2869 : : !scalars
2870 : : integer, intent(in) :: natom
2871 : : !arrays
2872 : : type(matlu_type), intent(inout) :: matlu1(natom)
2873 : : !Local variables-------------------------------
2874 : : !scalars
2875 : : integer :: iatom,im1,im2,ispinor2,ispinor1,isppol
2876 : : integer :: lpawu,ndim
2877 : : !arrays
2878 : : !************************************************************************
2879 [ # # ]: 0 : do iatom=1,natom
2880 : 0 : lpawu=matlu1(iatom)%lpawu
2881 [ # # ]: 0 : if(lpawu.ne.-1) then
2882 : 0 : ndim=2*lpawu+1
2883 [ # # ]: 0 : do isppol=1,matlu1(1)%nsppol
2884 [ # # ]: 0 : do ispinor1=1,matlu1(1)%nspinor
2885 [ # # ]: 0 : do ispinor2=1,matlu1(1)%nspinor
2886 [ # # ]: 0 : do im1=1,2*lpawu+1
2887 [ # # ]: 0 : do im2=1,2*lpawu+1
2888 : : matlu1(iatom)%mat(im1+(ispinor1-1)*ndim,im2+(ispinor2-1)*ndim,isppol)= &
2889 : 0 : & conjg(matlu1(iatom)%mat(im1+(ispinor1-1)*ndim,im2+(ispinor2-1)*ndim,isppol))
2890 : : enddo ! im2
2891 : : enddo ! im1
2892 : : enddo ! ispinor2
2893 : : enddo ! ispinor1
2894 : : enddo ! isppol
2895 : : endif ! lpawu
2896 : : enddo ! iatom
2897 : :
2898 : 0 : end subroutine conjg_matlu
2899 : : !!***
2900 : :
2901 : : !!****f* m_matlu/ln_matlu
2902 : : !! NAME
2903 : : !! ln_matlu
2904 : : !!
2905 : : !! FUNCTION
2906 : : !! Compute the logarithm of matlu (only if diagonal for the moment)
2907 : : !!
2908 : : !! COPYRIGHT
2909 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2910 : : !! This file is distributed under the terms of the
2911 : : !! GNU General Public License, see ~abinit/COPYING
2912 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2913 : : !!
2914 : : !! INPUTS
2915 : : !! matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity
2916 : :
2917 : : !! SIDE EFFECTS
2918 : : !!
2919 : : !! NOTES
2920 : : !!
2921 : : !! SOURCE
2922 : 0 : subroutine ln_matlu(matlu1,natom)
2923 : : use defs_wvltypes
2924 : :
2925 : : !Arguments ------------------------------------
2926 : : !scalars
2927 : : integer, intent(in) :: natom
2928 : : !arrays
2929 : : type(matlu_type), intent(inout) :: matlu1(natom)
2930 : : !Local variables-------------------------------
2931 : : !scalars
2932 : : integer :: iatom,im,ispinor,isppol
2933 : : integer :: lpawu,ndim
2934 : : character(len=500) :: message
2935 : : !arrays
2936 : : !************************************************************************
2937 : : !call checkdiag_matlu(matlu1,natom,tol8)
2938 [ # # ]: 0 : do iatom=1,natom
2939 : 0 : lpawu=matlu1(iatom)%lpawu
2940 [ # # ]: 0 : if(lpawu.ne.-1) then
2941 : 0 : ndim=2*lpawu+1
2942 [ # # ]: 0 : do isppol=1,matlu1(1)%nsppol
2943 [ # # ]: 0 : do ispinor=1,matlu1(1)%nspinor
2944 [ # # ]: 0 : do im=1,2*lpawu+1
2945 [ # # ]: 0 : if( real(matlu1(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol))<zero) then
2946 : 0 : write(message,'(2a,2es13.5,a)') ch10," ln_matlu: PROBLEM " &
2947 : 0 : & , matlu1(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol)
2948 : 0 : ABI_ERROR(message)
2949 : : endif
2950 : : matlu1(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol)= &
2951 : 0 : & log(matlu1(iatom)%mat(im+(ispinor-1)*ndim,im+(ispinor-1)*ndim,isppol))
2952 : : enddo ! im
2953 : : enddo ! ispinor
2954 : : enddo ! isppol
2955 : : endif ! lpawu
2956 : : enddo ! iatom
2957 : :
2958 : 0 : end subroutine ln_matlu
2959 : : !!***
2960 : :
2961 : : !!****f* m_matlu/slm2ylm_matlu
2962 : : !! NAME
2963 : : !! slm2ylm_matlu
2964 : : !!
2965 : : !! FUNCTION
2966 : : !! Transform mat from Slm to Ylm basis or vice versa
2967 : : !!
2968 : : !! COPYRIGHT
2969 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
2970 : : !! This file is distributed under the terms of the
2971 : : !! GNU General Public License, see ~abinit/COPYING
2972 : : !! or http://www.gnu.org/copyleft/gpl.txt .
2973 : : !!
2974 : : !! INPUTS
2975 : : !! matlu1(natom) :: input quantity
2976 : : !! natom :: number of atoms
2977 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
2978 : : !! option=1 go from Slm to Ylm basis
2979 : : !! option=2 go from Ylm to Slm basis
2980 : : !! SIDE EFFECTS
2981 : : !!
2982 : : !! NOTES
2983 : : !!
2984 : : !! SOURCE
2985 : :
2986 : 0 : subroutine slm2ylm_matlu(matlu,natom,paw_dmft,option,optprt)
2987 : :
2988 : : !Arguments ------------------------------------
2989 : : integer, intent(in) :: natom,option,optprt
2990 : : type(matlu_type), target, intent(inout) :: matlu(natom)
2991 : : type(paw_dmft_type), target, intent(in) :: paw_dmft
2992 : : !Local variables-------------------------------
2993 : : integer :: iatom,im1,im2,ispin,ispinor1,ispinor2,isppol
2994 : : integer :: lpawu,ndim,ndim_max,nspin,nspinor,nsppol
2995 : : complex(dp), pointer :: mat_out(:,:) => null(), slm2ylm(:,:) => null()
2996 : 0 : complex(dp), allocatable :: mat_inp(:,:),mat_tmp(:,:)
2997 : 0 : complex(dp), target, allocatable :: mat_tmp2(:,:)
2998 : : character(len=1) :: c1,c2
2999 : : character(len=500) :: message
3000 : : !************************************************************************
3001 : :
3002 : 0 : ndim_max = 2*paw_dmft%maxlpawu + 1
3003 : 0 : nspinor = paw_dmft%nspinor
3004 : 0 : nsppol = paw_dmft%nsppol
3005 : 0 : nspin = nsppol * (nspinor**2)
3006 : :
3007 [ # # ]: 0 : if (option == 1) then
3008 : 0 : c1 = "n" ; c2 = "c"
3009 [ # # ]: 0 : else if (option == 2) then
3010 : 0 : c1 = "c" ; c2 = "n"
3011 : : end if
3012 : :
3013 [ # # ]: 0 : do iatom=1,natom
3014 : :
3015 : 0 : lpawu = paw_dmft%lpawu(iatom)
3016 [ # # ]: 0 : if (lpawu == -1) cycle
3017 : 0 : ndim = 2*lpawu + 1
3018 : 0 : slm2ylm => paw_dmft%slm2ylm(:,:,lpawu+1)
3019 : :
3020 [ # # ]: 0 : if (optprt > 2) then
3021 : 0 : write(message,'(2a)') ch10,"SLM2YLM matrix"
3022 : 0 : call wrtout(std_out,message,'COLL')
3023 [ # # ]: 0 : do im1=1,ndim
3024 : 0 : write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') (slm2ylm(im1,im2),im2=1,ndim)
3025 : 0 : call wrtout(std_out,message,'COLL')
3026 : : end do ! im1
3027 : : end if ! optprt>2
3028 : :
3029 [ # # # # : 0 : ABI_MALLOC(mat_inp,(ndim,ndim))
# # # # #
# # # # #
# # ]
3030 [ # # # # : 0 : ABI_MALLOC(mat_tmp,(ndim,nspin*ndim))
# # # # #
# # # # #
# # ]
3031 [ # # # # : 0 : ABI_MALLOC(mat_tmp2,(ndim,nspin*ndim))
# # # # #
# # # ]
3032 : :
3033 : 0 : ispin = 0
3034 : :
3035 [ # # ]: 0 : do isppol=1,nsppol
3036 [ # # ]: 0 : do ispinor2=1,nspinor
3037 [ # # ]: 0 : do ispinor1=1,nspinor
3038 : :
3039 : 0 : ispin = ispin + 1
3040 : :
3041 : : ! Make copy here instead of creating a temporary when calling zgemm in order to please -fcheck
3042 [ # # # # ]: 0 : mat_inp(:,:) = matlu(iatom)%mat(1+(ispinor1-1)*ndim:ndim*ispinor1,1+(ispinor2-1)*ndim:ispinor2*ndim,isppol)
3043 : :
3044 [ # # ]: 0 : if (optprt > 2) then
3045 : 0 : write(message,'(2a,i2,a,i2,a,i2)') ch10,"SLM input matrix,&
3046 : 0 : & isppol=",isppol,", ispinor1=",ispinor1,", ispinor2=",ispinor2
3047 : 0 : call wrtout(std_out,message,'COLL')
3048 [ # # ]: 0 : do im1=1,ndim
3049 : : write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') &
3050 : 0 : & (mat_inp(im1,im2),im2=1,ndim)
3051 : 0 : call wrtout(std_out,message,'COLL')
3052 : : end do ! im1
3053 : : end if ! optprt>2
3054 : :
3055 : : call abi_xgemm("n",c2,ndim,ndim,ndim,cone,mat_inp(:,:),ndim,slm2ylm(:,1:ndim), &
3056 [ # # ]: 0 : & ndim_max,czero,mat_tmp(:,1+(ispin-1)*ndim:ndim*ispin),ndim)
3057 : :
3058 : : end do ! ispinor2
3059 : : end do ! ispinor1
3060 : : end do ! isppol
3061 : :
3062 : : call abi_xgemm(c1,"n",ndim,ndim*nspin,ndim,cone,slm2ylm(:,1:ndim),ndim_max,mat_tmp(:,:),ndim, &
3063 [ # # ]: 0 : & czero,mat_tmp2(:,:),ndim)
3064 : :
3065 : 0 : ispin = 0
3066 : :
3067 [ # # ]: 0 : do isppol=1,nsppol
3068 [ # # ]: 0 : do ispinor2=1,nspinor
3069 [ # # ]: 0 : do ispinor1=1,nspinor
3070 : :
3071 : 0 : ispin = ispin + 1
3072 : :
3073 : 0 : mat_out => mat_tmp2(:,1+(ispin-1)*ndim:ndim*ispin)
3074 : :
3075 [ # # # # : 0 : matlu(iatom)%mat(1+(ispinor1-1)*ndim:ndim*ispinor1,1+(ispinor2-1)*ndim:ispinor2*ndim,isppol) = mat_out(:,:)
# # # # #
# ]
3076 : :
3077 [ # # ]: 0 : if (optprt > 2) then
3078 : 0 : write(message,'(2a,i2,a,i2,a,i2)') ch10,"YLM output matrix, isppol=",isppol,", ispinor=",ispinor1,&
3079 : 0 : & ", ispinor2=",ispinor2
3080 : 0 : call wrtout(std_out,message,'COLL')
3081 [ # # ]: 0 : do im1=1,ndim
3082 : 0 : write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') (mat_out(im1,im2),im2=1,ndim)
3083 : 0 : call wrtout(std_out,message,'COLL')
3084 : : end do ! im1
3085 : : end if ! optprt>2
3086 : :
3087 : : end do ! ispinor1
3088 : : end do ! ispinor2
3089 : : end do ! isppol
3090 : :
3091 [ # # ]: 0 : ABI_FREE(mat_inp)
3092 [ # # ]: 0 : ABI_FREE(mat_tmp)
3093 [ # # ]: 0 : ABI_FREE(mat_tmp2)
3094 : :
3095 : : end do ! iatom
3096 : :
3097 : 0 : mat_out => null()
3098 : 0 : slm2ylm => null()
3099 : :
3100 : : !do iatom=1,natom
3101 : : ! lpawu=matlu(iatom)%lpawu
3102 : : ! if(lpawu.ne.-1) then
3103 : : ! ndim=2*lpawu+1
3104 : : ! ll=lpawu
3105 : : ! ABI_MALLOC(slm2ylm,(2*ll+1,2*ll+1))
3106 : : ! slm2ylm=czero
3107 : : ! do im=1,2*ll+1
3108 : : ! mm=im-ll-1;jm=-mm+ll+1
3109 : : ! onem=dble((-1)**mm)
3110 : : ! if (mm> 0) then
3111 : : ! slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
3112 : : ! slm2ylm(jm,im)= cmplx(invsqrt2, zero,kind=dp)
3113 : : ! end if
3114 : : ! if (mm==0) then
3115 : : ! slm2ylm(im,im)=cone
3116 : : ! end if
3117 : : ! if (mm< 0) then
3118 : : ! slm2ylm(im,im)= cmplx(zero, invsqrt2,kind=dp)
3119 : : ! slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
3120 : : ! end if
3121 : : ! end do
3122 : : ! if(optprt>2) then
3123 : : ! write(message,'(2a)') ch10,"SLM2YLM matrix"
3124 : : ! call wrtout(std_out,message,'COLL')
3125 : : ! do im1=1,ll*2+1
3126 : : ! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3127 : : !& (slm2ylm(im1,im2),im2=1,ll*2+1)
3128 : : ! call wrtout(std_out,message,'COLL')
3129 : : ! end do
3130 : : ! endif
3131 : : ! do isppol=1,matlu(1)%nsppol
3132 : : ! do ispinor=1,matlu(1)%nspinor
3133 : : ! do ispinor2=1,matlu(1)%nspinor
3134 : : ! ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1))
3135 : : ! ABI_MALLOC(mat_inp_c,(2*ll+1,2*ll+1))
3136 : : ! mat_inp_c(:,:) = matlu(iatom)%mat(1+(ispinor-1)*ndim:ndim+(ispinor-1)*ndim,1+(ispinor2-1)*ndim:ndim+(ispinor2-1)*ndim,isppol)
3137 : : ! mat_out_c=czero
3138 : :
3139 : : ! if(optprt>2) then
3140 : : ! write(message,'(2a, i2, a, i2, a, i2)') ch10,"SLM input matrix, isppol=", isppol, ", ispinor=", ispinor,&
3141 : : !& ", ispinor2=", ispinor2
3142 : : ! call wrtout(std_out,message,'COLL')
3143 : : ! do im1=1,ll*2+1
3144 : : ! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3145 : : !& (mat_inp_c(im1,im2),im2=1,ll*2+1)
3146 : : ! call wrtout(std_out,message,'COLL')
3147 : : ! end do
3148 : : ! endif
3149 : :
3150 : : ! do jm=1,2*ll+1
3151 : : ! do im=1,2*ll+1
3152 : : ! tmp2=czero
3153 : : ! do ii=1,2*ll+1
3154 : : ! do jj=1,2*ll+1
3155 : : ! if(option==1) then
3156 : : ! tmp2=tmp2+mat_inp_c(ii,jj)*(slm2ylm(im,ii))*CONJG(slm2ylm(jm,jj))
3157 : : ! else if(option==2) then
3158 : : ! tmp2=tmp2+mat_inp_c(ii,jj)*CONJG(slm2ylm(ii,im))*(slm2ylm(jj,jm))
3159 : : ! end if
3160 : : ! end do
3161 : : ! end do
3162 : : ! mat_out_c(im,jm)=tmp2
3163 : : ! end do
3164 : : ! end do
3165 : :
3166 : : ! if(optprt>2) then
3167 : : ! write(message,'(2a, i2, a, i2, a, i2)') ch10,"YLM output matrix, isppol=", isppol, ", ispinor=", ispinor,&
3168 : : !& ", ispinor2=", ispinor2
3169 : : ! call wrtout(std_out,message,'COLL')
3170 : : ! do im1=1,ll*2+1
3171 : : ! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3172 : : ! & (mat_out_c(im1,im2),im2=1,ll*2+1)
3173 : : ! call wrtout(std_out,message,'COLL')
3174 : : ! end do
3175 : : ! endif
3176 : :
3177 : : ! matlu(iatom)%mat(1+(ispinor-1)*ndim:ndim+(ispinor-1)*ndim,1+(ispinor2-1)*ndim:ndim+(ispinor2-1)*ndim,isppol)=mat_out_c(:,:)
3178 : : ! ABI_FREE(mat_out_c)
3179 : : ! ABI_FREE(mat_inp_c)
3180 : : ! enddo ! im
3181 : : ! enddo ! ispinor
3182 : : ! enddo ! isppol
3183 : : ! ABI_FREE(slm2ylm)
3184 : : ! endif ! lpawu
3185 : : !enddo ! iatom
3186 : :
3187 : 0 : end subroutine slm2ylm_matlu
3188 : : !!***
3189 : :
3190 : : !!****f* m_matlu/fac_matlu
3191 : : !! NAME
3192 : : !! fac_matlu
3193 : : !!
3194 : : !! FUNCTION
3195 : : !! Multiply matlu by a scalar
3196 : : !!
3197 : : !! COPYRIGHT
3198 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
3199 : : !! This file is distributed under the terms of the
3200 : : !! GNU General Public License, see ~abinit/COPYING
3201 : : !! or http://www.gnu.org/copyleft/gpl.txt .
3202 : : !!
3203 : : !! INPUTS
3204 : : !! matlu(natom) :: input quantity
3205 : : !! natom=number of atoms in cell.
3206 : : !! fac= factor
3207 : : !!
3208 : : !! OUTPUT
3209 : : !! matlu(natom) :: fac * matlu
3210 : : !!
3211 : : !! SIDE EFFECTS
3212 : : !!
3213 : : !! NOTES
3214 : : !!
3215 : : !! SOURCE
3216 : :
3217 : 19177 : subroutine fac_matlu(matlu,natom,fac)
3218 : :
3219 : : !Arguments ------------------------------------
3220 : : integer, intent(in) :: natom
3221 : : type(matlu_type), intent(inout) :: matlu(natom)
3222 : : complex(dp), intent(in) :: fac
3223 : : !Local variables-------------------------------
3224 : : integer :: iatom,lpawu
3225 : : ! character(len=500) :: message
3226 : : !************************************************************************
3227 : :
3228 [ + + ]: 40706 : do iatom=1,natom
3229 : 21529 : lpawu = matlu(iatom)%lpawu
3230 [ + + ]: 21529 : if (lpawu == -1) cycle
3231 [ + + + + : 1163524 : matlu(iatom)%mat(:,:,:) = fac * matlu(iatom)%mat(:,:,:)
+ + ]
3232 : : end do ! iatom
3233 : :
3234 : 19177 : end subroutine fac_matlu
3235 : : !!***
3236 : :
3237 : : !!****f* m_matlu/printplot_matlu
3238 : : !! NAME
3239 : : !! printplot_matlu
3240 : : !!
3241 : : !! FUNCTION
3242 : : !! Write matlu for a given frequency
3243 : : !!
3244 : : !! COPYRIGHT
3245 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
3246 : : !! This file is distributed under the terms of the
3247 : : !! GNU General Public License, see ~abinit/COPYING
3248 : : !! or http://www.gnu.org/copyleft/gpl.txt .
3249 : : !!
3250 : : !! INPUTS
3251 : : !! matlu(natom) :: input quantity to write
3252 : : !! natom=number of atoms in cell.
3253 : : !! freq :: frequency
3254 : : !! char1 :: name of the file on which to write
3255 : : !! units :: unit of the file
3256 : : !! imre :: if present, write real and imaginary parts in two different files
3257 : : !!
3258 : : !! OUTPUT
3259 : : !!
3260 : : !! SIDE EFFECTS
3261 : : !!
3262 : : !! NOTES
3263 : : !!
3264 : : !! SOURCE
3265 : :
3266 : 0 : subroutine printplot_matlu(matlu,natom,freq,char1,units,imre)
3267 : :
3268 : : !Arguments ------------------------------------
3269 : : !scalars
3270 : : integer, intent(in) :: natom,units
3271 : : integer, optional, intent(in) :: imre
3272 : : real(dp), intent(in) :: freq
3273 : : !arrays
3274 : : type(matlu_type), intent(in) :: matlu(natom)
3275 : : character(len=*), intent(in) :: char1
3276 : : !Local variables-------------------------------
3277 : : !scalars
3278 : : integer :: iatom,im,im1,ispinor,ispinor1,isppol
3279 : : integer :: lpawu,ndim,nspinor,nsppol,unitnb
3280 : : character(len=4) :: tag_at
3281 : : character(len=fnlen) :: tmpfil,tmpfilim,tmpfilre
3282 : : ! character(len=500) :: message
3283 : : !arrays
3284 : : !************************************************************************
3285 : :
3286 : 0 : nspinor = matlu(1)%nspinor
3287 : 0 : nsppol = matlu(1)%nsppol
3288 : :
3289 : : ! not yet tested and used
3290 [ # # ]: 0 : do iatom=1,natom
3291 : 0 : lpawu = matlu(iatom)%lpawu
3292 [ # # ]: 0 : if (lpawu == -1) cycle
3293 : 0 : ndim = 2*lpawu + 1
3294 : 0 : unitnb = units + iatom
3295 : 0 : call int2char4(iatom,tag_at)
3296 [ # # ]: 0 : if (present(imre)) then
3297 [ # # ]: 0 : tmpfilre = trim(char1)//tag_at//"re"
3298 [ # # ]: 0 : tmpfilim = trim(char1)//tag_at//"im"
3299 : 0 : open(unit=unitnb+10,file=trim(tmpfilre),status='unknown',form='formatted')
3300 : 0 : open(unit=unitnb+20,file=trim(tmpfilim),status='unknown',form='formatted')
3301 : 0 : write(unitnb+10,'(400e26.16)') freq,(((((dble(matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol)),&
3302 [ # # # # : 0 : & im=1,ndim),ispinor=1,nspinor),im1=1,ndim),ispinor1=1,nspinor),isppol=1,nsppol)
# # # # #
# # # # #
# # # # #
# ]
3303 : 0 : write(unitnb+20,'(400e26.16)') freq,(((((aimag(matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol)),&
3304 [ # # # # : 0 : & im=1,ndim),ispinor=1,nspinor),im1=1,ndim),ispinor1=1,nspinor),isppol=1,nsppol)
# # # # #
# # # # #
# # # # #
# ]
3305 : : else
3306 [ # # ]: 0 : tmpfil = trim(char1)//tag_at
3307 : 0 : open(unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
3308 : 0 : write(unitnb,'(400e26.16)') freq,(((((matlu(iatom)%mat(im+(ispinor-1)*ndim,im1+(ispinor1-1)*ndim,isppol),&
3309 [ # # # # : 0 : & im=1,ndim),ispinor=1,nspinor),im1=1,ndim),ispinor1=1,nspinor),isppol=1,nsppol)
# # # # #
# # # # #
# # # # #
# ]
3310 : : end if ! present(imre)
3311 : : end do ! iatom
3312 : :
3313 : 0 : end subroutine printplot_matlu
3314 : : !!***
3315 : :
3316 : : !!****f* m_matlu/identity_matlu
3317 : : !! NAME
3318 : : !! identity_matlu
3319 : : !!
3320 : : !! FUNCTION
3321 : : !! Set the diagonal elements to 1 (the off-diagonal are not set to 0)
3322 : : !!
3323 : : !! COPYRIGHT
3324 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
3325 : : !! This file is distributed under the terms of the
3326 : : !! GNU General Public License, see ~abinit/COPYING
3327 : : !! or http://www.gnu.org/copyleft/gpl.txt .
3328 : : !!
3329 : : !! INPUTS
3330 : : !! matlu(natom) :: input quantity
3331 : : !! natom=number of atoms in cell.
3332 : : !!
3333 : : !! OUTPUT
3334 : : !!
3335 : : !! SIDE EFFECTS
3336 : : !!
3337 : : !! NOTES
3338 : : !!
3339 : : !! SOURCE
3340 : :
3341 : 267 : subroutine identity_matlu(matlu,natom)
3342 : :
3343 : : !Arguments ------------------------------------
3344 : : integer, intent(in) :: natom
3345 : : type(matlu_type), intent(inout) :: matlu(natom)
3346 : : !Local variables-------------------------------
3347 : : integer :: iatom,im,lpawu,ndim,nspinor
3348 : : ! character(len=500) :: message
3349 : : !arrays
3350 : : !************************************************************************
3351 : :
3352 : 267 : nspinor = matlu(1)%nspinor
3353 : :
3354 [ + + ]: 1215 : do iatom=1,natom
3355 : 948 : lpawu = matlu(iatom)%lpawu
3356 [ + + ]: 948 : if (lpawu == -1) cycle
3357 : 303 : ndim = nspinor * (2*lpawu+1)
3358 [ + + ]: 1926 : do im=1,ndim
3359 [ + + ]: 4300 : matlu(iatom)%mat(im,im,:) = cone
3360 : : end do ! im
3361 : : end do ! iatom
3362 : :
3363 : 267 : end subroutine identity_matlu
3364 : : !!***
3365 : :
3366 : : !!***
3367 : : !!****f* m_matlu/magmomforb_matlu
3368 : : !! NAME
3369 : : !! magmomforb_matlu
3370 : : !!
3371 : : !! FUNCTION
3372 : : !! return the product of occupation matrix of dimension [(2*ll+1)]**4 in the Ylm basis
3373 : : !! with the matrix of orbital angular momentum element for the x,y and z direction.
3374 : : !! Option gives the direction of the magnetic moment x==1, y==2 and z==3.
3375 : : !!
3376 : : !!
3377 : : !! COPYRIGHT
3378 : : !! Copyright (C) 2005-2026 ABINIT group (FGendron)
3379 : : !! This file is distributed under the terms of the
3380 : : !! GNU General Public License, see ~abinit/COPYING
3381 : : !! or http://www.gnu.org/copyleft/gpl.txt .
3382 : : !!
3383 : : !! INPUTS
3384 : : !! matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity in Ylm basis
3385 : : !! natom :: number of atoms
3386 : : !! option = 1 :: x axis
3387 : : !! = 2 :: y axis
3388 : : !! = 3 :: z axis
3389 : : !! optptr > 2 :: print orbital angular matrix elements and resulting product
3390 : : !!
3391 : : !! OUTPUT
3392 : : !! matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: product
3393 : : !!
3394 : : !! SIDE EFFECTS
3395 : : !!
3396 : : !! NOTES
3397 : : !!
3398 : : !! SOURCE
3399 : 0 : subroutine magmomforb_matlu(matlu,mu,natom,option,optprt)
3400 : : use defs_wvltypes
3401 : :
3402 : : !Arguments ------------------------------------
3403 : : !scalars
3404 : : integer, intent(in) :: natom,option,optprt
3405 : : complex(dp), allocatable, intent(inout) :: mu(:)
3406 : : !arrays
3407 : : type(matlu_type), intent(inout) :: matlu(natom)
3408 : : !Local variables-------------------------------
3409 : : !scalars
3410 : : integer :: iatom,im,ispinor,isppol,ispinor2
3411 : : integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,im1,im2,tndim,ndim
3412 : : character(len=500) :: message
3413 : : real(dp) :: xj
3414 : : !arrays
3415 : 0 : complex(dp),allocatable :: mat_out_c(:,:)
3416 : : ! integer, allocatable :: ind_msml(:,:)
3417 : 0 : complex(dp), allocatable :: temp_mat(:,:)
3418 : 0 : type(coeff2c_type), allocatable :: gathermatlu(:)
3419 : 0 : type(coeff2c_type), allocatable :: muorb(:)
3420 : : !************************************************************************
3421 : :
3422 : : !=====================================
3423 : : ! Allocate matrices
3424 : : !=====================================
3425 : :
3426 [ # # # # : 0 : ABI_MALLOC(gathermatlu,(natom))
# # # # #
# # # #
# ]
3427 [ # # # # : 0 : ABI_MALLOC(muorb,(natom))
# # # # #
# ]
3428 [ # # ]: 0 : do iatom=1,natom
3429 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
3430 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
3431 [ # # # # : 0 : ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
# # # # #
# # # # #
# # # # ]
3432 [ # # # # ]: 0 : gathermatlu(iatom)%value=czero
3433 [ # # # # : 0 : ABI_MALLOC(muorb(iatom)%value,(tndim,tndim))
# # # # #
# # # #
# ]
3434 [ # # # # ]: 0 : muorb(iatom)%value=czero
3435 : : end if
3436 : : end do
3437 : :
3438 [ # # ]: 0 : do iatom=1,natom
3439 : 0 : lpawu=matlu(iatom)%lpawu
3440 [ # # ]: 0 : if(lpawu.ne.-1) then
3441 : :
3442 : 0 : ll=lpawu
3443 : 0 : lcor=lpawu
3444 : : !=====================================
3445 : : !build orbital angular momentum matrix along x axis
3446 : : !=====================================
3447 [ # # ]: 0 : if(option==1) then
3448 : :
3449 : : jc1=0
3450 [ # # ]: 0 : do ms1 =-1,1
3451 : 0 : xj=float(ms1)+half
3452 [ # # ]: 0 : do ml1 = -ll,ll
3453 : 0 : jc1=jc1+1
3454 [ # # ]: 0 : if(jc1 == 1) then
3455 : 0 : muorb(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
3456 : : endif
3457 [ # # # # ]: 0 : if(jc1 > 1 .and. jc1 < 2*(2*ll+1)) then
3458 : 0 : muorb(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
3459 : 0 : muorb(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
3460 [ # # ]: 0 : else if(jc1== 2*(2*ll+1)) then
3461 : 0 : muorb(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
3462 : : end if
3463 : : end do
3464 : : end do
3465 : :
3466 : : !=====================================
3467 : : !build orbital angular momentum matrix along y axis
3468 : : !=====================================
3469 [ # # ]: 0 : else if(option==2) then
3470 : :
3471 : : jc1=0
3472 [ # # ]: 0 : do ms1 =-1,1
3473 : 0 : xj=float(ms1)+half
3474 [ # # ]: 0 : do ml1 = -ll,ll
3475 : 0 : jc1=jc1+1
3476 [ # # ]: 0 : if(jc1 == 1) then
3477 : 0 : muorb(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
3478 : : endif
3479 [ # # # # ]: 0 : if(jc1 > 1 .and. jc1 < 2*(2*ll+1)) then
3480 : 0 : muorb(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
3481 : 0 : muorb(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp)
3482 [ # # ]: 0 : else if(jc1 == 2*(2*ll+1)) then
3483 : 0 : muorb(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp)
3484 : : end if
3485 : : end do
3486 : : end do
3487 : :
3488 : : !=====================================
3489 : : !build orbital angular momentum matrix along z axis
3490 : : !=====================================
3491 [ # # ]: 0 : else if(option==3) then
3492 : : jc1=0
3493 [ # # ]: 0 : do ms1=-1,1
3494 [ # # ]: 0 : do ml1=-ll,ll
3495 : 0 : jc1=jc1+1
3496 [ # # ]: 0 : if(jc1 < tndim+1) then
3497 : 0 : muorb(iatom)%value(jc1,jc1) = ml1
3498 : : endif
3499 : : end do
3500 : : end do
3501 : : end if
3502 : :
3503 [ # # ]: 0 : if(optprt>2) then
3504 : 0 : write(message,'(a,i4)') "Orbital angular momentum matrix elements in |m_l,m_s> basis for axis=", option
3505 : 0 : call wrtout(std_out,message,"COLL")
3506 [ # # ]: 0 : do im=1,2*(ll*2+1)
3507 [ # # # # ]: 0 : write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muorb(iatom)%value(im,jm),jm=1,2*(ll*2+1))
3508 : 0 : call wrtout(std_out,message,"COLL")
3509 : : end do
3510 : : end if
3511 : :
3512 : : end if !lpawu
3513 : : end do !atom
3514 : :
3515 : : !=====================================
3516 : : ! Reshape input Ylm matlu in one 14x14 matrix
3517 : : !=====================================
3518 : :
3519 : 0 : call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
3520 : :
3521 : : !!printing for debug
3522 : : ! write(std_out,*) "gathermatlu in magmomforb"
3523 : : ! do im1=1,tndim
3524 : : ! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3525 : : ! (gathermatlu(1)%value(im1,im2),im2=1,tndim)
3526 : : ! call wrtout(std_out,message,'COLL')
3527 : : ! end do
3528 : :
3529 : : !=====================================
3530 : : ! Matrix product of Occ and muorb
3531 : : !=====================================
3532 [ # # ]: 0 : do iatom=1,natom
3533 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
3534 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
3535 [ # # # # : 0 : ABI_MALLOC(temp_mat,(tndim,tndim))
# # # # #
# # # # #
# # ]
3536 : :
3537 : : call abi_zgemm('n','n',tndim,tndim,tndim,cone,&
3538 : : & gathermatlu(iatom)%value,tndim,&
3539 : : & muorb(iatom)%value,tndim,czero,&
3540 : 0 : & temp_mat,tndim)
3541 : :
3542 [ # # # # : 0 : gathermatlu(iatom)%value=temp_mat
# # # # #
# # # ]
3543 [ # # ]: 0 : ABI_FREE(temp_mat)
3544 : :
3545 : : !=====================================
3546 : : ! Trace of matrix product
3547 : : !=====================================
3548 : :
3549 [ # # ]: 0 : do im1=1,tndim
3550 [ # # ]: 0 : do im2=1,tndim
3551 [ # # ]: 0 : if(im1==im2) then
3552 : 0 : mu(iatom) = mu(iatom) + gathermatlu(iatom)%value(im1,im2)
3553 : : end if
3554 : : end do
3555 : : end do
3556 : :
3557 : :
3558 : : !=====================================
3559 : : ! Reshape product matrix into matlu format
3560 : : !=====================================
3561 : :
3562 : : !call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
3563 : :
3564 : :
3565 [ # # ]: 0 : if(optprt>2) then
3566 [ # # # # : 0 : ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1))
# # # # #
# # # # #
# # # # ]
3567 : 0 : ndim = 2*ll+1
3568 [ # # ]: 0 : do isppol=1,matlu(1)%nsppol
3569 [ # # ]: 0 : do ispinor=1,matlu(1)%nspinor
3570 [ # # ]: 0 : do ispinor2=1,matlu(1)%nspinor
3571 [ # # # # ]: 0 : mat_out_c(:,:) = matlu(iatom)%mat(1+(ispinor-1)*ndim:ndim+(ispinor-1)*ndim,1+(ispinor2-1)*ndim:ndim+(ispinor2-1)*ndim,isppol)
3572 : :
3573 : 0 : write(message,'(2a, i2, a, i2, a, i2)') ch10,"Orbital angular momentum matrix, isppol=", isppol, ", ispinor=",&
3574 : 0 : & ispinor,", ispinor2=", ispinor2
3575 : 0 : call wrtout(std_out,message,'COLL')
3576 [ # # ]: 0 : do im1=1,ll*2+1
3577 : : write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3578 : 0 : & (mat_out_c(im1,im2),im2=1,ll*2+1)
3579 : 0 : call wrtout(std_out,message,'COLL')
3580 : : end do
3581 : :
3582 : : end do ! ispinor2
3583 : : end do ! ispinor
3584 : : end do ! isppol
3585 [ # # ]: 0 : ABI_FREE(mat_out_c)
3586 : : endif
3587 : :
3588 : : end if !lpawu
3589 : : end do !atom
3590 : :
3591 : : !=====================================
3592 : : ! Deallocate gathermatlu
3593 : : !=====================================
3594 : :
3595 [ # # ]: 0 : do iatom=1,natom
3596 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
3597 [ # # ]: 0 : ABI_FREE(gathermatlu(iatom)%value)
3598 [ # # ]: 0 : ABI_FREE(muorb(iatom)%value)
3599 : : end if
3600 : : end do
3601 [ # # # # : 0 : ABI_FREE(gathermatlu)
# # # # ]
3602 [ # # # # : 0 : ABI_FREE(muorb)
# # # # ]
3603 : :
3604 : 0 : end subroutine magmomforb_matlu
3605 : :
3606 : : !!***
3607 : :
3608 : :
3609 : : !!***
3610 : : !!****f* m_matlu/magmomfspin_matlu
3611 : : !! NAME
3612 : : !! magmomfspin_matlu
3613 : : !!
3614 : : !! FUNCTION
3615 : : !! return the product of occupation matrix of dimension [(2*ll+1)]**4 in the Ylm basis
3616 : : !! with the matrix of spin angular momentum element for the x,y and z direction.
3617 : : !! Option gives the direction of the magnetic moment x==1, y==2 and z==3.
3618 : : !!
3619 : : !!
3620 : : !! COPYRIGHT
3621 : : !! Copyright (C) 2005-2026 ABINIT group (FGendron)
3622 : : !! This file is distributed under the terms of the
3623 : : !! GNU General Public License, see ~abinit/COPYING
3624 : : !! or http://www.gnu.org/copyleft/gpl.txt .
3625 : : !!
3626 : : !! INPUTS
3627 : : !! matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity in Ylm basis
3628 : : !! natom :: number of atoms
3629 : : !! option = 1 :: x axis
3630 : : !! = 2 :: y axis
3631 : : !! = 3 :: z axis
3632 : : !! optptr > 2 :: print spin angular matrix elements and resulting product
3633 : : !!
3634 : : !! OUTPUT
3635 : : !! matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: product
3636 : : !!
3637 : : !! SIDE EFFECTS
3638 : : !!
3639 : : !! NOTES
3640 : : !!
3641 : : !! SOURCE
3642 : 0 : subroutine magmomfspin_matlu(matlu,mu,natom,option,optprt)
3643 : : use defs_wvltypes
3644 : :
3645 : : !Arguments ------------------------------------
3646 : : !scalars
3647 : : integer, intent(in) :: natom,option,optprt
3648 : : complex(dp), allocatable, intent(inout) :: mu(:)
3649 : : !arrays
3650 : : type(matlu_type), intent(inout) :: matlu(natom)
3651 : : !Local variables-------------------------------
3652 : : !scalars
3653 : : integer :: iatom,im,ispinor,isppol,ispinor2
3654 : : integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,im1,im2,tndim,ndim
3655 : : character(len=500) :: message
3656 : : real(dp) :: xj
3657 : : !arrays
3658 : 0 : complex(dp),allocatable :: mat_out_c(:,:)
3659 : 0 : integer, allocatable :: ind_msml(:,:)
3660 : 0 : complex(dp), allocatable :: temp_mat(:,:)
3661 : 0 : type(coeff2c_type), allocatable :: gathermatlu(:)
3662 : 0 : type(coeff2c_type), allocatable :: muspin(:)
3663 : : !************************************************************************
3664 : :
3665 : : !=====================================
3666 : : ! Allocate matrices
3667 : : !=====================================
3668 : :
3669 [ # # # # : 0 : ABI_MALLOC(gathermatlu,(natom))
# # # # #
# # # #
# ]
3670 [ # # # # : 0 : ABI_MALLOC(muspin,(natom))
# # # # #
# ]
3671 [ # # ]: 0 : do iatom=1,natom
3672 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
3673 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
3674 [ # # # # : 0 : ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
# # # # #
# # # # #
# # # # ]
3675 [ # # # # ]: 0 : gathermatlu(iatom)%value=czero
3676 [ # # # # : 0 : ABI_MALLOC(muspin(iatom)%value,(tndim,tndim))
# # # # #
# # # #
# ]
3677 [ # # # # ]: 0 : muspin(iatom)%value=czero
3678 : : end if
3679 : : end do
3680 : :
3681 [ # # ]: 0 : do iatom=1,natom
3682 : 0 : lpawu=matlu(iatom)%lpawu
3683 [ # # ]: 0 : if(lpawu.ne.-1) then
3684 : 0 : ll=lpawu
3685 : 0 : lcor=lpawu
3686 [ # # # # : 0 : ABI_MALLOC(ind_msml,(2,-ll:ll))
# # # # #
# # # ]
3687 [ # # # # ]: 0 : ind_msml=czero
3688 : : !=====================================
3689 : : !build spin angular momentum matrix along x axis
3690 : : !=====================================
3691 [ # # ]: 0 : if(option==1) then
3692 : : jc1=0
3693 [ # # ]: 0 : do ms1=1,2
3694 [ # # ]: 0 : do ml1=-ll,ll
3695 : 0 : jc1=jc1+1
3696 : 0 : ind_msml(ms1,ml1)=jc1
3697 : : end do
3698 : : end do
3699 : :
3700 : : jc1=0
3701 [ # # ]: 0 : do ms1 =-1,1
3702 : 0 : xj=float(ms1)+half
3703 [ # # ]: 0 : do ml1 = -ll,ll
3704 : 0 : jc1=jc1+1
3705 [ # # ]: 0 : if(xj < 0.0 ) then
3706 : 0 : muspin(iatom)%value(ind_msml(2,ml1),jc1) = 0.5
3707 [ # # ]: 0 : else if(xj > 0.0) then
3708 : 0 : muspin(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = 0.5
3709 : : end if
3710 : : end do
3711 : : end do
3712 : :
3713 : : !=====================================
3714 : : !build spin angular momentum matrix along y axis
3715 : : !up spin is first
3716 : : !=====================================
3717 [ # # ]: 0 : else if(option==2) then
3718 : : jc1=0
3719 [ # # ]: 0 : do ms1=1,2
3720 [ # # ]: 0 : do ml1=-ll,ll
3721 : 0 : jc1=jc1+1
3722 : 0 : ind_msml(ms1,ml1)=jc1
3723 : : end do
3724 : : end do
3725 : :
3726 : : jc1=0
3727 [ # # ]: 0 : do ms1 =-1,1
3728 : 0 : xj=float(ms1)+half
3729 [ # # ]: 0 : do ml1 = -ll,ll
3730 : 0 : jc1=jc1+1
3731 [ # # ]: 0 : if(xj < 0.0 ) then
3732 : 0 : muspin(iatom)%value(ind_msml(2,ml1),jc1) = cmplx(zero,0.5,kind=dp)
3733 [ # # ]: 0 : else if(xj > 0.0) then
3734 : 0 : muspin(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = cmplx(zero,-0.5,kind=dp)
3735 : : end if
3736 : : end do
3737 : : end do
3738 : :
3739 : : !=====================================
3740 : : !build spin angular momentum matrix along z axis
3741 : : !up spin is first
3742 : : !=====================================
3743 [ # # ]: 0 : else if(option==3) then
3744 : : jc1=0
3745 [ # # ]: 0 : do ms1=-1,1
3746 : 0 : xj=float(ms1)+half
3747 [ # # ]: 0 : do ml1=-ll,ll
3748 : 0 : jc1=jc1+1
3749 [ # # ]: 0 : if(jc1 < tndim+1) then
3750 [ # # ]: 0 : if(xj < 0.0 ) then
3751 : 0 : muspin(iatom)%value(jc1,jc1) = -xj
3752 [ # # ]: 0 : else if(xj > 0.0) then
3753 : 0 : muspin(iatom)%value(jc1,jc1) = -xj
3754 : : end if
3755 : : endif
3756 : : end do
3757 : : end do
3758 : : end if
3759 : :
3760 [ # # ]: 0 : if(optprt>2) then
3761 : 0 : write(message,'(a,i4)') "Spin angular momentum matrix elements in |m_l,m_s> basis for axis", option
3762 : 0 : call wrtout(std_out,message,"COLL")
3763 [ # # ]: 0 : do im=1,2*(ll*2+1)
3764 [ # # # # ]: 0 : write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muspin(iatom)%value(im,jm),jm=1,2*(ll*2+1))
3765 : 0 : call wrtout(std_out,message,"COLL")
3766 : : end do
3767 : : end if
3768 : :
3769 [ # # ]: 0 : ABI_FREE(ind_msml)
3770 : : end if !lpawu
3771 : : end do !atom
3772 : :
3773 : : !=====================================
3774 : : ! Reshape input Ylm matlu in one 14x14 matrix
3775 : : !=====================================
3776 : :
3777 : 0 : call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
3778 : :
3779 : : !!printing for debug
3780 : : !! write(std_out,*) "gathermatlu in magmomfspin"
3781 : : !! do im1=1,tndim
3782 : : !! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3783 : : !! (gathermatlu(1)%value(im1,im2),im2=1,tndim)
3784 : : !! call wrtout(std_out,message,'coll')
3785 : : !! end do
3786 : :
3787 : : !=====================================
3788 : : ! Matrix product of Occ and muspin
3789 : : !=====================================
3790 [ # # ]: 0 : do iatom=1,natom
3791 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
3792 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
3793 [ # # # # : 0 : ABI_MALLOC(temp_mat,(tndim,tndim))
# # # # #
# # # # #
# # ]
3794 : :
3795 : : call zgemm('n','n',tndim,tndim,tndim,cone,&
3796 : : & gathermatlu(iatom)%value,tndim,&
3797 : : & muspin(iatom)%value,tndim,czero,&
3798 : 0 : & temp_mat,tndim)
3799 : :
3800 [ # # # # : 0 : gathermatlu(iatom)%value=temp_mat
# # # # #
# # # ]
3801 [ # # ]: 0 : ABI_FREE(temp_mat)
3802 : :
3803 : : !!printing for debug
3804 : : !!write(std_out,*) "gathermatlu in magmomfspin after product"
3805 : : !!do im1=1,tndim
3806 : : !! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3807 : : !! (gathermatlu(1)%value(im1,im2),im2=1,tndim)
3808 : : !! call wrtout(std_out,message,'coll')
3809 : : !!end do
3810 : :
3811 : : !=====================================
3812 : : ! Trace of matrix product
3813 : : !=====================================
3814 : :
3815 [ # # ]: 0 : do im1=1,tndim
3816 [ # # ]: 0 : do im2=1,tndim
3817 [ # # ]: 0 : if(im1==im2) then
3818 : 0 : mu(iatom) = mu(iatom) + gathermatlu(iatom)%value(im1,im2)
3819 : : end if
3820 : : end do
3821 : : end do
3822 : :
3823 : : !=====================================
3824 : : ! Reshape product matrix into matlu format
3825 : : !=====================================
3826 : :
3827 : : !call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
3828 : :
3829 : :
3830 : : !=====================================
3831 : : ! Print matlu
3832 : : !=====================================
3833 [ # # ]: 0 : if(optprt>2) then
3834 [ # # # # : 0 : ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1))
# # # # #
# # # # #
# # # # ]
3835 : 0 : ndim = 2*ll+1
3836 [ # # ]: 0 : do isppol=1,matlu(1)%nsppol
3837 [ # # ]: 0 : do ispinor=1,matlu(1)%nspinor
3838 [ # # ]: 0 : do ispinor2=1,matlu(1)%nspinor
3839 [ # # # # ]: 0 : mat_out_c(:,:) = matlu(iatom)%mat(1+(ispinor-1)*ndim:ndim+(ispinor-1)*ndim,1+(ispinor2-1)*ndim:ndim+(ispinor2-1)*ndim,isppol)
3840 : :
3841 : 0 : write(message,'(2a, i2, a, i2, a, i2)') ch10,"Spin angular momentum matrix, isppol=", isppol, ", ispinor=", ispinor,&
3842 : 0 : & ", ispinor2=", ispinor2
3843 : 0 : call wrtout(std_out,message,'COLL')
3844 [ # # ]: 0 : do im1=1,ll*2+1
3845 : : write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
3846 : 0 : & (mat_out_c(im1,im2),im2=1,ll*2+1)
3847 : 0 : call wrtout(std_out,message,'COLL')
3848 : : end do
3849 : :
3850 : : end do ! im
3851 : : end do ! ispinor
3852 : : end do ! isppol
3853 [ # # ]: 0 : ABI_FREE(mat_out_c)
3854 : : endif
3855 : :
3856 : : end if !lpawu
3857 : : end do !atom
3858 : :
3859 : : !=====================================
3860 : : ! Deallocate gathermatlu
3861 : : !=====================================
3862 : :
3863 [ # # ]: 0 : do iatom=1,natom
3864 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
3865 [ # # ]: 0 : ABI_FREE(gathermatlu(iatom)%value)
3866 [ # # ]: 0 : ABI_FREE(muspin(iatom)%value)
3867 : : end if
3868 : : end do
3869 [ # # # # : 0 : ABI_FREE(gathermatlu)
# # # # ]
3870 [ # # # # : 0 : ABI_FREE(muspin)
# # # # ]
3871 : :
3872 : 0 : end subroutine magmomfspin_matlu
3873 : :
3874 : : !!***
3875 : :
3876 : :
3877 : : !!***
3878 : : !!****f* m_matlu/magmomfzeeman_matlu
3879 : : !! NAME
3880 : : !! magmomfspin_matlu
3881 : : !!
3882 : : !! FUNCTION
3883 : : !! return the product of occupation matrix of dimension [(2*ll+1)]**4 in the Ylm basis
3884 : : !! with the matrix of Zeeman angular momentum element (L_u + 2*S_u) for the x,y and z direction.
3885 : : !! Option gives the direction of the magnetic moment x==1, y==2 and z==3.
3886 : : !!
3887 : : !!
3888 : : !! COPYRIGHT
3889 : : !! Copyright (C) 2005-2026 ABINIT group (FGendron)
3890 : : !! This file is distributed under the terms of the
3891 : : !! GNU General Public License, see ~abinit/COPYING
3892 : : !! or http://www.gnu.org/copyleft/gpl.txt .
3893 : : !!
3894 : : !! INPUTS
3895 : : !! matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity in Ylm basis
3896 : : !! natom :: number of atoms
3897 : : !! option = 1 :: x axis
3898 : : !! = 2 :: y axis
3899 : : !! = 3 :: z axis
3900 : : !! optptr > 2 :: print Zeeman angular matrix elements and resulting product
3901 : : !!
3902 : : !! OUTPUT
3903 : : !! matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: product
3904 : : !!
3905 : : !! SIDE EFFECTS
3906 : : !!
3907 : : !! NOTES
3908 : : !!
3909 : : !! SOURCE
3910 : 0 : subroutine magmomfzeeman_matlu(matlu,mu,natom,option,optprt)
3911 : : use defs_wvltypes
3912 : :
3913 : : !Arguments ------------------------------------
3914 : : !scalars
3915 : : integer, intent(in) :: natom,option,optprt
3916 : : complex(dp), allocatable, intent(inout) :: mu(:)
3917 : : !arrays
3918 : : type(matlu_type), intent(inout) :: matlu(natom)
3919 : : !Local variables-------------------------------
3920 : : !scalars
3921 : : integer :: iatom,im,ispinor,isppol,ispinor2
3922 : : integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,im1,im2,tndim,ndim
3923 : : character(len=500) :: message
3924 : : real(dp) :: xj
3925 : : !arrays
3926 : 0 : complex(dp),allocatable :: mat_out_c(:,:)
3927 : 0 : integer, allocatable :: ind_msml(:,:)
3928 : 0 : complex(dp), allocatable :: temp_mat(:,:)
3929 : 0 : type(coeff2c_type), allocatable :: gathermatlu(:)
3930 : 0 : type(coeff2c_type), allocatable :: muzeeman(:)
3931 : : !************************************************************************
3932 : :
3933 : : !=====================================
3934 : : ! Allocate matrices
3935 : : !=====================================
3936 : :
3937 [ # # # # : 0 : ABI_MALLOC(gathermatlu,(natom))
# # # # #
# # # #
# ]
3938 [ # # # # : 0 : ABI_MALLOC(muzeeman,(natom))
# # # # #
# ]
3939 [ # # ]: 0 : do iatom=1,natom
3940 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
3941 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
3942 [ # # # # : 0 : ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
# # # # #
# # # # #
# # # # ]
3943 [ # # # # ]: 0 : gathermatlu(iatom)%value=czero
3944 [ # # # # : 0 : ABI_MALLOC(muzeeman(iatom)%value,(tndim,tndim))
# # # # #
# # # #
# ]
3945 [ # # # # ]: 0 : muzeeman(iatom)%value=czero
3946 : : end if
3947 : : end do
3948 : :
3949 [ # # ]: 0 : do iatom=1,natom
3950 : 0 : lpawu=matlu(iatom)%lpawu
3951 [ # # ]: 0 : if(lpawu.ne.-1) then
3952 : 0 : ll=lpawu
3953 : 0 : lcor=lpawu
3954 [ # # # # : 0 : ABI_MALLOC(ind_msml,(2,-ll:ll))
# # # # #
# # # ]
3955 [ # # # # ]: 0 : ind_msml=czero
3956 : : !=====================================
3957 : : !build Zeeman angular momentum matrix along x axis
3958 : : !=====================================
3959 [ # # ]: 0 : if(option==1) then
3960 : : jc1=0
3961 [ # # ]: 0 : do ms1=1,2
3962 [ # # ]: 0 : do ml1=-ll,ll
3963 : 0 : jc1=jc1+1
3964 : 0 : ind_msml(ms1,ml1)=jc1
3965 : : end do
3966 : : end do
3967 : :
3968 : : jc1=0
3969 [ # # ]: 0 : do ms1 =-1,1
3970 : 0 : xj=float(ms1)+half
3971 [ # # ]: 0 : do ml1 = -ll,ll
3972 : 0 : jc1=jc1+1
3973 [ # # ]: 0 : if(xj < 0.0 ) then
3974 : 0 : muzeeman(iatom)%value(ind_msml(2,ml1),jc1) = 2*0.5
3975 [ # # ]: 0 : else if(xj > 0.0) then
3976 : 0 : muzeeman(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = 2*0.5
3977 : : end if
3978 [ # # ]: 0 : if(jc1 == 1) then
3979 : 0 : muzeeman(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
3980 : : endif
3981 [ # # # # ]: 0 : if(jc1 > 1 .and.jc1 < 2*(2*ll+1)) then
3982 : 0 : muzeeman(iatom)%value(jc1,jc1+1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5
3983 : 0 : muzeeman(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
3984 [ # # ]: 0 : else if(jc1 == 2*(2*ll+1)) then
3985 : 0 : muzeeman(iatom)%value(jc1,jc1-1) = sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5
3986 : : end if
3987 : : end do
3988 : : end do
3989 : :
3990 : : !=====================================
3991 : : !build Zeeman angular momentum matrix along y axis
3992 : : !up spin is first
3993 : : !=====================================
3994 [ # # ]: 0 : else if(option==2) then
3995 : : jc1=0
3996 [ # # ]: 0 : do ms1=1,2
3997 [ # # ]: 0 : do ml1=-ll,ll
3998 : 0 : jc1=jc1+1
3999 : 0 : ind_msml(ms1,ml1)=jc1
4000 : : end do
4001 : : end do
4002 : :
4003 : : jc1=0
4004 [ # # ]: 0 : do ms1 =-1,1
4005 : 0 : xj=float(ms1)+half
4006 [ # # ]: 0 : do ml1 = -ll,ll
4007 : 0 : jc1=jc1+1
4008 [ # # ]: 0 : if(xj < 0.0 ) then
4009 : 0 : muzeeman(iatom)%value(ind_msml(2,ml1),jc1) = 2*cmplx(zero,0.5,kind=dp)
4010 [ # # ]: 0 : else if(xj > 0.0) then
4011 : 0 : muzeeman(iatom)%value(ind_msml(1,ml1),ind_msml(2,ml1)) = 2*cmplx(zero,-0.5,kind=dp)
4012 : : end if
4013 [ # # ]: 0 : if(jc1 == 1) then
4014 : 0 : muzeeman(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
4015 : : endif
4016 [ # # # # ]: 0 : if(jc1 > 1 .and. jc1 < 2*(2*ll+1)) then
4017 : 0 : muzeeman(iatom)%value(jc1,jc1+1) = cmplx(zero,sqrt(float((lcor*(lcor+1)) - ml1*(ml1 + 1)))*0.5,kind=dp)
4018 : 0 : muzeeman(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp)
4019 [ # # ]: 0 : else if(jc1== 2*(2*ll+1)) then
4020 : 0 : muzeeman(iatom)%value(jc1,jc1-1) = cmplx(zero,-sqrt(float((lcor*(lcor+1)) - ml1*(ml1 - 1)))*0.5,kind=dp)
4021 : : end if
4022 : : end do
4023 : : end do
4024 : :
4025 : :
4026 : : !=====================================
4027 : : !build Zeeman angular momentum matrix along z axis
4028 : : !up spin is first
4029 : : !=====================================
4030 [ # # ]: 0 : else if(option==3) then
4031 : : jc1=0
4032 [ # # ]: 0 : do ms1=-1,1
4033 : 0 : xj=float(ms1)+half
4034 [ # # ]: 0 : do ml1=-ll,ll
4035 : 0 : jc1=jc1+1
4036 [ # # ]: 0 : if(jc1 < tndim+1) then
4037 [ # # ]: 0 : if(xj < 0.0 ) then
4038 : 0 : muzeeman(iatom)%value(jc1,jc1) = ml1-2*xj
4039 [ # # ]: 0 : else if(xj > 0.0) then
4040 : 0 : muzeeman(iatom)%value(jc1,jc1) = ml1-2*xj
4041 : : end if
4042 : : endif
4043 : : end do
4044 : : end do
4045 : : end if
4046 : :
4047 : :
4048 [ # # ]: 0 : if(optprt>2) then
4049 : 0 : write(message,'(a,i4)') "Zeeman angular momentum matrix elements in |m_l,m_s> basis for axis", option
4050 : 0 : call wrtout(std_out,message,"COLL")
4051 [ # # ]: 0 : do im=1,2*(ll*2+1)
4052 [ # # # # ]: 0 : write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muzeeman(iatom)%value(im,jm),jm=1,2*(ll*2+1))
4053 : 0 : call wrtout(std_out,message,"COLL")
4054 : : end do
4055 : : end if
4056 : :
4057 [ # # ]: 0 : ABI_FREE(ind_msml)
4058 : : end if !lpawu
4059 : : end do !atom
4060 : :
4061 : : !=====================================
4062 : : ! Reshape input Ylm matlu in one 14x14 matrix
4063 : : !=====================================
4064 : :
4065 : : !ABI_MALLOC(gathermatlu,(natom))
4066 : : !do iatom=1,natom
4067 : : ! if(matlu(iatom)%lpawu.ne.-1) then
4068 : : ! tndim=2*(2*matlu(iatom)%lpawu+1)
4069 : : ! ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
4070 : : ! gathermatlu(iatom)%value=czero
4071 : : ! end if
4072 : : !end do
4073 : :
4074 : 0 : call gather_matlu(matlu,gathermatlu,natom,option=1,prtopt=1)
4075 : :
4076 : : !!printing for debug
4077 : : !! write(std_out,*) "gathermatlu in magmomfspin"
4078 : : !! do im1=1,tndim
4079 : : !! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
4080 : : !! (gathermatlu(1)%value(im1,im2),im2=1,tndim)
4081 : : !! call wrtout(std_out,message,'coll')
4082 : : !! end do
4083 : :
4084 : : !=====================================
4085 : : ! Matrix product of Occ and muzeeman
4086 : : !=====================================
4087 [ # # ]: 0 : do iatom=1,natom
4088 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
4089 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
4090 [ # # # # : 0 : ABI_MALLOC(temp_mat,(tndim,tndim))
# # # # #
# # # # #
# # ]
4091 : :
4092 : : call abi_zgemm('n','n',tndim,tndim,tndim,cone,&
4093 : : & gathermatlu(iatom)%value,tndim,&
4094 : : & muzeeman(iatom)%value,tndim,czero,&
4095 : 0 : & temp_mat,tndim)
4096 : :
4097 [ # # # # : 0 : gathermatlu(iatom)%value=temp_mat
# # # # #
# # # ]
4098 [ # # ]: 0 : ABI_FREE(temp_mat)
4099 : :
4100 : : !!printing for debug
4101 : : !!write(std_out,*) "gathermatlu in magmomfspin after product"
4102 : : !!do im1=1,tndim
4103 : : !! write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
4104 : : !! (gathermatlu(1)%value(im1,im2),im2=1,tndim)
4105 : : !! call wrtout(std_out,message,'coll')
4106 : : !!end do
4107 : :
4108 : : !=====================================
4109 : : ! Trace of matrix product
4110 : : !=====================================
4111 : :
4112 [ # # ]: 0 : do im1=1,tndim
4113 [ # # ]: 0 : do im2=1,tndim
4114 [ # # ]: 0 : if(im1==im2) then
4115 : 0 : mu(iatom) = mu(iatom) + gathermatlu(iatom)%value(im1,im2)
4116 : : end if
4117 : : end do
4118 : : end do
4119 : :
4120 : : !=====================================
4121 : : ! Reshape product matrix into matlu format
4122 : : !=====================================
4123 : :
4124 : : !call gather_matlu(matlu,gathermatlu,natom,option=-1,prtopt=1)
4125 : :
4126 : :
4127 : : !=====================================
4128 : : ! Print matlu
4129 : : !=====================================
4130 [ # # ]: 0 : if(optprt>2) then
4131 [ # # # # : 0 : ABI_MALLOC(mat_out_c,(2*ll+1,2*ll+1))
# # # # #
# # # # #
# # # # ]
4132 : 0 : ndim = 2*ll+1
4133 [ # # ]: 0 : do isppol=1,matlu(1)%nsppol
4134 [ # # ]: 0 : do ispinor=1,matlu(1)%nspinor
4135 [ # # ]: 0 : do ispinor2=1,matlu(1)%nspinor
4136 [ # # # # ]: 0 : mat_out_c(:,:) = matlu(iatom)%mat(1+(ispinor-1)*ndim:ndim+(ispinor-1)*ndim,1+(ispinor2-1)*ndim:ndim+(ispinor2-1)*ndim,isppol)
4137 : :
4138 : 0 : write(message,'(2a, i2, a, i2, a, i2)') ch10,"Zeeman angular momentum matrix, isppol=", isppol, ", ispinor=",&
4139 : 0 : & ispinor,", ispinor2=", ispinor2
4140 : 0 : call wrtout(std_out,message,'COLL')
4141 [ # # ]: 0 : do im1=1,ll*2+1
4142 : : write(message,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
4143 : 0 : & (mat_out_c(im1,im2),im2=1,ll*2+1)
4144 : 0 : call wrtout(std_out,message,'COLL')
4145 : : end do
4146 : :
4147 : : end do ! im
4148 : : end do ! ispinor
4149 : : end do ! isppol
4150 [ # # ]: 0 : ABI_FREE(mat_out_c)
4151 : : endif ! optprt
4152 : :
4153 : : end if !lpawu
4154 : : end do !atom
4155 : :
4156 : : !=====================================
4157 : : ! Deallocate gathermatlu
4158 : : !=====================================
4159 : :
4160 [ # # ]: 0 : do iatom=1,natom
4161 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
4162 [ # # ]: 0 : ABI_FREE(gathermatlu(iatom)%value)
4163 [ # # ]: 0 : ABI_FREE(muzeeman(iatom)%value)
4164 : : end if
4165 : : end do
4166 [ # # # # : 0 : ABI_FREE(gathermatlu)
# # # # ]
4167 [ # # # # : 0 : ABI_FREE(muzeeman)
# # # # ]
4168 : 0 : end subroutine magmomfzeeman_matlu
4169 : :
4170 : : !!***
4171 : :
4172 : : !!***
4173 : : !!****f* m_matlu/chi_matlu
4174 : : !! NAME
4175 : : !! chi_matlu
4176 : : !!
4177 : : !! FUNCTION
4178 : : !! return the matrix of dimension [(2*ll+1)]**4 in the Ylm basis
4179 : : !! with the matrix elements of the orbital (option=1), spin (option=2) and
4180 : : !! total (option=3) angular momentum for z direction. It is used by the
4181 : : !! QMC for the correlation function of the magnetic moment
4182 : : !!
4183 : : !!
4184 : : !! COPYRIGHT
4185 : : !! Copyright (C) 2005-2026 ABINIT group (FGendron)
4186 : : !! This file is distributed under the terms of the
4187 : : !! GNU General Public License, see ~abinit/COPYING
4188 : : !! or http://www.gnu.org/copyleft/gpl.txt .
4189 : : !!
4190 : : !! INPUTS
4191 : : !! matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) ::
4192 : : !! natom :: number of atoms
4193 : : !! option = 1 :: Orbital angular momentum along z axis
4194 : : !! = 2 :: 2*Spin angluar momentum alonf z axis
4195 : : !! = 3 :: total angular momentum along z axis
4196 : : !! optptr > 2 :: print angular matrix elements
4197 : : !!
4198 : : !! OUTPUT
4199 : : !! matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: quantity in Ylm basis
4200 : : !!
4201 : : !! SIDE EFFECTS
4202 : : !!
4203 : : !! NOTES
4204 : : !!
4205 : : !! SOURCE
4206 : 0 : subroutine chi_matlu(matlu,natom,option,optprt)
4207 : : use defs_wvltypes
4208 : :
4209 : : !Arguments ------------------------------------
4210 : : !scalars
4211 : : integer, intent(in) :: natom,option,optprt
4212 : : !arrays
4213 : : type(matlu_type), intent(inout) :: matlu(natom)
4214 : : !Local variables-------------------------------
4215 : : !scalars
4216 : : integer :: iatom,im
4217 : : integer :: lpawu,ll,jm,ml1,jc1,ms1,lcor,tndim
4218 : : character(len=500) :: message
4219 : : real(dp) :: xj
4220 : : !arrays
4221 : 0 : integer, allocatable :: ind_msml(:,:)
4222 : : ! type(coeff2c_type), allocatable :: gathermatlu(:)
4223 : 0 : type(coeff2c_type), allocatable :: muchi(:)
4224 : : !************************************************************************
4225 : :
4226 : : !=====================================
4227 : : ! Allocate matrices
4228 : : !=====================================
4229 : :
4230 [ # # # # : 0 : ABI_MALLOC(muchi,(natom))
# # # # #
# # # #
# ]
4231 [ # # ]: 0 : do iatom=1,natom
4232 : : !if(matlu(iatom)%lpawu.ne.-1) then
4233 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
4234 [ # # # # : 0 : ABI_MALLOC(muchi(iatom)%value,(tndim,tndim))
# # # # #
# # # # #
# # # # ]
4235 [ # # # # ]: 0 : muchi(iatom)%value=czero
4236 : : !end if
4237 : : end do
4238 : :
4239 [ # # ]: 0 : do iatom=1,natom
4240 : 0 : lpawu=matlu(iatom)%lpawu
4241 [ # # ]: 0 : if(lpawu.ne.-1) then
4242 : 0 : ll=lpawu
4243 : 0 : lcor=lpawu
4244 [ # # # # : 0 : ABI_MALLOC(ind_msml,(2,-ll:ll))
# # # # #
# # # ]
4245 [ # # # # ]: 0 : ind_msml=czero
4246 : : !=====================================
4247 : : !Orbital angular momentum matrix along z axis
4248 : : !=====================================
4249 [ # # ]: 0 : if(option==1) then
4250 : : jc1=0
4251 [ # # ]: 0 : do ms1=1,2
4252 [ # # ]: 0 : do ml1=-ll,ll
4253 : 0 : jc1=jc1+1
4254 : 0 : ind_msml(ms1,ml1)=jc1
4255 : : end do
4256 : : end do
4257 : :
4258 : : jc1=0
4259 [ # # ]: 0 : do ms1=-1,1
4260 [ # # ]: 0 : do ml1=-ll,ll
4261 : 0 : jc1=jc1+1
4262 [ # # ]: 0 : if(jc1 <tndim+1) then
4263 : 0 : muchi(iatom)%value(jc1,jc1) = ml1
4264 : : endif
4265 : : end do
4266 : : end do
4267 : :
4268 : : !=====================================
4269 : : !Spin angular momentum matrix along z axis
4270 : : !up spin is first
4271 : : !=====================================
4272 [ # # ]: 0 : else if(option==2) then
4273 : : jc1=0
4274 [ # # ]: 0 : do ms1=1,2
4275 [ # # ]: 0 : do ml1=-ll,ll
4276 : 0 : jc1=jc1+1
4277 : 0 : ind_msml(ms1,ml1)=jc1
4278 : : end do
4279 : : end do
4280 : :
4281 : : jc1=0
4282 [ # # ]: 0 : do ms1=-1,1
4283 : 0 : xj=float(ms1)+half
4284 [ # # ]: 0 : do ml1=-ll,ll
4285 : 0 : jc1=jc1+1
4286 [ # # ]: 0 : if(jc1<tndim+1) then
4287 [ # # ]: 0 : if(xj < 0.0 ) then
4288 : 0 : muchi(iatom)%value(jc1,jc1) = -2*xj
4289 [ # # ]: 0 : else if(xj > 0.0) then
4290 : 0 : muchi(iatom)%value(jc1,jc1) = -2*xj
4291 : :
4292 : : end if
4293 : : endif
4294 : : end do
4295 : : end do
4296 : :
4297 : : !=====================================
4298 : : !Total angular momentum matrix along z axis
4299 : : !up spin is first
4300 : : !=====================================
4301 [ # # ]: 0 : else if(option==3) then
4302 : : jc1=0
4303 [ # # ]: 0 : do ms1=-1,1
4304 : 0 : xj=float(ms1)+half
4305 [ # # ]: 0 : do ml1=-ll,ll
4306 : 0 : jc1=jc1+1
4307 [ # # ]: 0 : if(jc1<tndim+1) then
4308 [ # # ]: 0 : if(xj < 0.0 ) then
4309 : 0 : muchi(iatom)%value(jc1,jc1) = ml1-2*xj
4310 [ # # ]: 0 : else if(xj > 0.0) then
4311 : 0 : muchi(iatom)%value(jc1,jc1) = ml1-2*xj
4312 : : end if
4313 : : endif
4314 : : end do
4315 : : end do
4316 : : end if
4317 : :
4318 [ # # ]: 0 : if(optprt>2) then
4319 [ # # ]: 0 : if(option==1) then
4320 : 0 : write(message,'(a)') "Orbital angular momentum matrix elements in |m_l,m_s> basis"
4321 [ # # ]: 0 : else if(option==2) then
4322 : 0 : write(message,'(a)') "Spin angular momentum matrix elements in |m_l,m_s> basis"
4323 [ # # ]: 0 : else if(option==3) then
4324 : 0 : write(message,'(a)') "Zeeman angular momentum matrix elements in |m_l,m_s> basis"
4325 : : end if
4326 : 0 : call wrtout(std_out,message,"COLL")
4327 [ # # ]: 0 : do im=1,2*(ll*2+1)
4328 [ # # # # ]: 0 : write(message,'(6(1x,9(1x,f4.1,",",f4.1)))') (muchi(iatom)%value(im,jm),jm=1,2*(ll*2+1))
4329 : 0 : call wrtout(std_out,message,"COLL")
4330 : : end do
4331 : : end if
4332 : :
4333 [ # # ]: 0 : ABI_FREE(ind_msml)
4334 : :
4335 : : !=====================================
4336 : : ! Reshape matrix into matlu format
4337 : : !=====================================
4338 : :
4339 : : ! call gather_matlu(matlu,muchi(iatom),natom=1,option=-1,prtopt=1)
4340 : :
4341 : :
4342 : : end if !lpawu
4343 : : end do !atom
4344 : :
4345 : : !=====================================
4346 : : ! Reshape matrix into matlu format
4347 : : !=====================================
4348 : :
4349 : 0 : call gather_matlu(matlu,muchi,natom,option=-1,prtopt=1)
4350 : :
4351 : : !=====================================
4352 : : ! Deallocate gathermatlu
4353 : : !=====================================
4354 : :
4355 [ # # ]: 0 : do iatom=1,natom
4356 [ # # ]: 0 : if(matlu(iatom)%lpawu.ne.-1) then
4357 [ # # ]: 0 : ABI_FREE(muchi(iatom)%value)
4358 : : end if
4359 : : end do
4360 [ # # # # : 0 : ABI_FREE(muchi)
# # # # ]
4361 : :
4362 : 0 : end subroutine chi_matlu
4363 : :
4364 : : !!***
4365 : :
4366 : : !!****f* m_matlu/trace_prod_matlu
4367 : : !! NAME
4368 : : !! trace_prod_matlu
4369 : : !!
4370 : : !! FUNCTION
4371 : : !! Computes Tr(matlu1*matlu2) for each atom. It is NOT assumed
4372 : : !! that either matlu1 or matlu2 is symmetric, so this routine is
4373 : : !! suboptimal if this is the case.
4374 : : !!
4375 : : !! COPYRIGHT
4376 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
4377 : : !! This file is distributed under the terms of the
4378 : : !! GNU General Public License, see ~abinit/COPYING
4379 : : !! or http://www.gnu.org/copyleft/gpl.txt .
4380 : : !!
4381 : : !! INPUTS
4382 : : !! matlu1(natom),matlu2(natom) :: input quantity
4383 : : !! natom=number of atoms in cell.
4384 : : !! opt_add = if present, add the result to trace
4385 : : !! trace_tot = if present, computes the trace over all atoms
4386 : : !! iatom = if present, only computes the contribution of iatom
4387 : : !! if <=0, add all contributions as usual
4388 : : !!
4389 : : !! OUTPUT
4390 : : !! trace(natom) :: Tr(matlu1*matlu2) for each atom
4391 : : !!
4392 : : !! SIDE EFFECTS
4393 : : !!
4394 : : !! NOTES
4395 : : !!
4396 : : !! SOURCE
4397 : :
4398 : 8138 : subroutine trace_prod_matlu(matlu1,matlu2,natom,trace,trace_tot,iatom)
4399 : :
4400 : : !Arguments ------------------------------------
4401 : : integer, intent(in) :: natom
4402 : : type(matlu_type), intent(in) :: matlu1(natom),matlu2(natom)
4403 : : complex(dp), intent(inout) :: trace(natom)
4404 : : complex(dp), optional, intent(out) :: trace_tot
4405 : : integer, optional, intent(in) :: iatom
4406 : : !Local variables-------------------------------
4407 : : integer :: ia1,ia2,iatom_,isppol,lpawu,nspinor,nsppol
4408 : : !************************************************************************
4409 : :
4410 : 8138 : nspinor = matlu1(1)%nspinor
4411 : 8138 : nsppol = matlu1(1)%nsppol
4412 : :
4413 [ + + ]: 24172 : trace(:) = czero
4414 : :
4415 : 8138 : ia1 = 1 ; ia2 = natom
4416 [ + + ]: 8138 : if (present(iatom)) then
4417 [ + - ]: 8032 : if (iatom > 0) then
4418 : 8138 : ia1 = iatom ; ia2 = iatom
4419 : : end if
4420 : : end if ! present(iatom)
4421 : :
4422 [ + + ]: 24172 : do iatom_=ia1,ia2
4423 : 16034 : lpawu = matlu1(iatom_)%lpawu
4424 [ + + ]: 16034 : if (lpawu == -1) cycle
4425 [ + + ]: 23662 : do isppol=1,nsppol
4426 [ + + + + ]: 495878 : trace(iatom_) = trace(iatom_) + sum(matlu1(iatom_)%mat(:,:,isppol)*transpose(matlu2(iatom_)%mat(:,:,isppol)))
4427 : : end do ! isppol
4428 [ + + ]: 16684 : if (nsppol == 1 .and. nspinor == 1) trace(iatom_) = trace(iatom_) * two
4429 : : end do ! iatom
4430 : :
4431 [ - + - - ]: 8138 : if (present(trace_tot)) trace_tot = sum(trace(:))
4432 : :
4433 : 8138 : end subroutine trace_prod_matlu
4434 : : !!***
4435 : :
4436 : : !!****f* m_matlu/xmpi_matlu
4437 : : !! NAME
4438 : : !! xmpi_sum_matlu
4439 : : !!
4440 : : !! FUNCTION
4441 : : !! Put matlu into a buffer and perform the required MPI operation
4442 : : !! (xmpi_sum or xmpi_bcast) on the input communicator.
4443 : : !!
4444 : : !! INPUTS
4445 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
4446 : : !! natom = number of atoms
4447 : : !! comm = communicator on which to perform the required communication
4448 : : !! master = master node (only used in the case of xmpi_bcast, default is 0)
4449 : : !! option = 1 (default) : xmpi_sum
4450 : : !! = 2 : xmpi_bcast
4451 : : !!
4452 : : !! OUTPUT
4453 : : !!
4454 : : !! SOURCE
4455 : :
4456 : 2134 : subroutine xmpi_matlu(matlu,natom,comm,master,option)
4457 : :
4458 : : !Arguments ------------------------------------
4459 : : integer, intent(in) :: comm,natom
4460 : : type(matlu_type), intent(inout) :: matlu(natom)
4461 : : integer, optional, intent(in) :: master,option
4462 : : !Local variables-------------------------------
4463 : : integer :: iatom,ibuf,im1,ierr,isppol,lpawu
4464 : : integer :: master_node,ndim,nspinor,nsppol,opt,siz_buf
4465 : 2134 : complex(dp), allocatable :: buffer(:)
4466 : : !************************************************************************
4467 : :
4468 : 2134 : nspinor = matlu(1)%nspinor
4469 : 2134 : nsppol = matlu(1)%nsppol
4470 : 2134 : opt = 1
4471 [ + + ]: 2134 : if (present(option)) opt = option
4472 : 2134 : master_node = 0
4473 [ + + ]: 2134 : if (present(master)) master_node = master
4474 : :
4475 : 2134 : siz_buf = 0
4476 [ + + ]: 9549 : do iatom=1,natom
4477 : 7415 : lpawu = matlu(iatom)%lpawu
4478 [ + + ]: 7415 : if (lpawu == -1) cycle
4479 : 9549 : siz_buf = siz_buf + (2*lpawu+1)**2
4480 : : end do ! iatom
4481 : :
4482 : 2134 : siz_buf = siz_buf * (nspinor**2) * nsppol
4483 : :
4484 [ - + + - : 6402 : ABI_MALLOC(buffer,(siz_buf))
+ - + - -
+ - + ]
4485 : :
4486 : 2134 : ibuf = 0
4487 [ + + ]: 9549 : do iatom=1,natom
4488 : 7415 : lpawu = matlu(iatom)%lpawu
4489 [ + + ]: 7415 : if (lpawu == -1) cycle
4490 : 2482 : ndim = nspinor * (2*lpawu+1)
4491 [ + + ]: 8345 : do isppol=1,nsppol
4492 [ + + ]: 28466 : do im1=1,ndim
4493 [ + + ]: 106114 : buffer(ibuf+1:ibuf+ndim) = matlu(iatom)%mat(:,im1,isppol)
4494 : 21051 : ibuf = ibuf + ndim
4495 : : end do ! im1
4496 : : end do ! isppol
4497 : : end do ! iatom
4498 : :
4499 [ + + ]: 2134 : if (opt == 1) then
4500 : 2096 : call xmpi_sum(buffer(:),comm,ierr)
4501 [ + - ]: 38 : else if (opt == 2) then
4502 : 38 : call xmpi_bcast(buffer(:),master_node,comm,ierr)
4503 : : end if ! opt
4504 : :
4505 : : ibuf = 0
4506 [ + + ]: 9549 : do iatom=1,natom
4507 : 7415 : lpawu = matlu(iatom)%lpawu
4508 [ + + ]: 7415 : if (lpawu == -1) cycle
4509 : 2482 : ndim = nspinor * (2*lpawu+1)
4510 [ + + ]: 8345 : do isppol=1,nsppol
4511 [ + + ]: 28466 : do im1=1,ndim
4512 [ + + ]: 106114 : matlu(iatom)%mat(:,im1,isppol) = buffer(ibuf+1:ibuf+ndim)
4513 : 21051 : ibuf = ibuf + ndim
4514 : : end do ! im1
4515 : : end do ! isppol
4516 : : end do ! iatom
4517 : :
4518 [ - + ]: 2134 : ABI_FREE(buffer)
4519 : :
4520 : 2134 : end subroutine xmpi_matlu
4521 : : !!***
4522 : :
4523 : : !!****f* m_matlu/symmetrize_matlu
4524 : : !! NAME
4525 : : !! symmetrize_matlu
4526 : : !!
4527 : : !! FUNCTION
4528 : : !! Symmetrizes matlu (A = (A+A^T)/2
4529 : : !!
4530 : : !! INPUTS
4531 : : !! matlu <type(matlu_type)>= density matrix in the local orbital basis and related variables
4532 : : !! natom = number of atoms
4533 : : !!
4534 : : !! OUTPUT
4535 : : !! err = difference between A and symmetrized A
4536 : : !!
4537 : : !! SOURCE
4538 : :
4539 : 0 : subroutine symmetrize_matlu(matlu,natom,err)
4540 : :
4541 : : !Arguments ------------------------------------
4542 : : integer, intent(in) :: natom
4543 : : type(matlu_type), intent(inout) :: matlu(natom)
4544 : : real(dp), optional, intent(out) :: err
4545 : : !Local variables-------------------------------
4546 : : integer :: iatom,isppol,lpawu,nspinor,nsppol,tndim
4547 : : real(dp) :: err_
4548 : 0 : complex(dp), allocatable :: mat_tmp(:,:)
4549 : : !************************************************************************
4550 : :
4551 : 0 : nspinor = matlu(1)%nspinor
4552 : 0 : nsppol = matlu(1)%nsppol
4553 : :
4554 [ # # ]: 0 : if (present(err)) err = zero
4555 : :
4556 [ # # ]: 0 : do iatom=1,natom
4557 : 0 : lpawu = matlu(iatom)%lpawu
4558 [ # # ]: 0 : if (lpawu == -1) cycle
4559 : 0 : tndim = nspinor * (2*lpawu+1)
4560 [ # # # # : 0 : ABI_MALLOC(mat_tmp,(tndim,tndim))
# # # # #
# # # # #
# # ]
4561 [ # # ]: 0 : do isppol=1,nsppol
4562 : : mat_tmp(:,:) = half * (matlu(iatom)%mat(:,:,isppol)+ &
4563 [ # # # # ]: 0 : & transpose(matlu(iatom)%mat(:,:,isppol)))
4564 [ # # ]: 0 : if (present(err)) then
4565 [ # # # # ]: 0 : err_ = sum(abs(mat_tmp(:,:)-matlu(iatom)%mat(:,:,isppol)))
4566 [ # # ]: 0 : if (err_ > err) err = err_
4567 : : end if ! present(err)
4568 [ # # # # ]: 0 : matlu(iatom)%mat(:,:,isppol) = mat_tmp(:,:)
4569 : : end do ! isppol
4570 [ # # ]: 0 : ABI_FREE(mat_tmp)
4571 : : end do ! iatom
4572 : :
4573 : 0 : end subroutine symmetrize_matlu
4574 : : !!***
4575 : :
4576 : : !!****f* m_matlu/ylm2jmj_matlu
4577 : : !! NAME
4578 : : !! ylm2jmj_matlu
4579 : : !!
4580 : : !! FUNCTION
4581 : : !! Transform mat from Ylm to JmJ basis or vice versa
4582 : : !!
4583 : : !! COPYRIGHT
4584 : : !! Copyright (C) 2005-2026 ABINIT group (BAmadon)
4585 : : !! This file is distributed under the terms of the
4586 : : !! GNU General Public License, see ~abinit/COPYING
4587 : : !! or http://www.gnu.org/copyleft/gpl.txt .
4588 : : !!
4589 : : !! INPUTS
4590 : : !! matlu(natom) :: input quantity
4591 : : !! natom :: number of atoms
4592 : : !! option=1 go from Ylm to JmJ basis
4593 : : !! option=2 go from JmJ to Ylm basis
4594 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
4595 : : !! SIDE EFFECTS
4596 : : !!
4597 : : !! NOTES
4598 : : !!
4599 : : !! SOURCE
4600 : :
4601 : 0 : subroutine ylm2jmj_matlu(matlu,natom,option,paw_dmft)
4602 : :
4603 : : !Arguments ------------------------------------
4604 : : integer, intent(in) :: natom,option
4605 : : type(matlu_type), intent(inout) :: matlu(natom)
4606 : : type(paw_dmft_type), intent(in) :: paw_dmft
4607 : : !Local variables-------------------------------
4608 : : integer :: iatom,lpawu,nspinor,tndim,tndim_max
4609 : 0 : complex(dp), allocatable :: mat_tmp(:,:)
4610 : : character(len=1) :: c1,c2
4611 : : !************************************************************************
4612 : :
4613 : 0 : nspinor = paw_dmft%nspinor
4614 [ # # ]: 0 : if (nspinor == 1) ABI_BUG("nspinor should be equal to 2")
4615 : :
4616 : 0 : tndim_max = nspinor * (2*paw_dmft%maxlpawu+1)
4617 : :
4618 [ # # ]: 0 : if (option == 1) then
4619 : 0 : c1 = "c" ; c2 = "n"
4620 : : else
4621 : 0 : c1 = "n" ; c2 = "c"
4622 : : end if
4623 : :
4624 [ # # # # : 0 : ABI_MALLOC(mat_tmp,(tndim_max,tndim_max))
# # # # #
# # # # #
# # ]
4625 : :
4626 [ # # ]: 0 : do iatom=1,natom
4627 : 0 : lpawu = paw_dmft%lpawu(iatom)
4628 [ # # ]: 0 : if (lpawu == -1) cycle
4629 [ # # ]: 0 : if (lpawu == 0) ABI_BUG("l should not be equal to 0")
4630 : 0 : tndim = nspinor * (2*lpawu+1)
4631 : :
4632 : : call abi_xgemm("n",c2,tndim,tndim,tndim,cone,matlu(iatom)%mat(:,:,1),tndim, &
4633 : 0 : & paw_dmft%jmj2ylm(:,1:tndim,lpawu+1),tndim_max,czero,mat_tmp(:,1:tndim),tndim_max)
4634 : :
4635 : : call abi_xgemm(c1,"n",tndim,tndim,tndim,cone,paw_dmft%jmj2ylm(:,1:tndim,lpawu+1), &
4636 : 0 : & tndim_max,mat_tmp(:,1:tndim),tndim_max,czero,matlu(iatom)%mat(:,:,1),tndim)
4637 : :
4638 : : end do ! iatom
4639 : :
4640 [ # # ]: 0 : ABI_FREE(mat_tmp)
4641 : :
4642 : 0 : end subroutine ylm2jmj_matlu
4643 : : !!***
4644 : :
4645 : : !!***
4646 : : !!****f* m_matlu/magnfield_matlu
4647 : : !! NAME
4648 : : !! magnfield_matlu
4649 : : !!
4650 : : !! FUNCTION
4651 : : !! return the matrix of magnetic moment mz times Bz
4652 : : !!
4653 : : !!
4654 : : !! COPYRIGHT
4655 : : !! Copyright (C) 2005-2026 ABINIT group (FGendron)
4656 : : !! This file is distributed under the terms of the
4657 : : !! GNU General Public License, see ~abinit/COPYING
4658 : : !! or http://www.gnu.org/copyleft/gpl.txt .
4659 : : !!
4660 : : !! INPUTS
4661 : : !! matlu1(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: input quantity in Ylm basis
4662 : : !! natom :: number of atoms
4663 : : !! bfield :: value of magnetic field in Tesla
4664 : : !! option = 1 :: scalar spin angular momentum along z axis
4665 : : !! option = 2 :: SOC total angular momentum (L+2S) along z axis
4666 : : !!
4667 : : !! OUTPUT
4668 : : !! matlu(natom)%(nsppol,nspinor,nspinor,ndim,ndim) :: product
4669 : : !!
4670 : : !! SIDE EFFECTS
4671 : : !!
4672 : : !! NOTES
4673 : : !!
4674 : : !! SOURCE
4675 : 8 : subroutine magnfield_matlu(matlu,natom,bfield,option)
4676 : : use defs_basis
4677 : : use defs_wvltypes
4678 : : implicit none
4679 : :
4680 : : !Arguments ------------------------------------
4681 : : !scalars
4682 : : integer, intent(in) :: natom,option
4683 : : real(dp) :: bfield
4684 : : !arrays
4685 : : type(matlu_type), intent(inout) :: matlu(natom)
4686 : : !Local variables-------------------------------
4687 : : !scalars
4688 : : integer :: iatom,im,ndim,isppol
4689 : : integer :: ll,ml1,jc1,ms1,tndim
4690 : : real(dp) :: xj
4691 : : !arrays
4692 : 8 : type(coeff2c_type), allocatable :: magnmatb(:)
4693 : : !************************************************************************
4694 : :
4695 : : !================================
4696 : : ! Allocate matrices
4697 : : !================================
4698 : :
4699 [ - + + - : 32 : ABI_MALLOC(magnmatb,(natom))
+ - + - -
+ - + +
+ ]
4700 [ + + ]: 16 : do iatom=1,natom
4701 [ + - ]: 16 : if(matlu(iatom)%lpawu .ne. -1) then
4702 : 8 : tndim=2*(2*matlu(iatom)%lpawu+1)
4703 [ - + + - : 32 : ABI_MALLOC(magnmatb(iatom)%value,(tndim,tndim))
- + + - +
- + - - +
- + - + ]
4704 [ + + + + ]: 888 : magnmatb(iatom)%value=czero
4705 : : endif
4706 : : enddo
4707 : :
4708 [ + - ]: 8 : if(option .eq. 1) then
4709 : :
4710 : : !================================
4711 : : ! Scalar magnetism (Spin only case)
4712 : : ! H = mu_B*g_e*S_Z*B_z
4713 : : !================================
4714 : :
4715 [ + + ]: 16 : do iatom=1,natom
4716 [ + - ]: 16 : if(matlu(iatom)%lpawu .ne. -1) then
4717 : 8 : ndim=2*matlu(iatom)%lpawu+1
4718 [ + + ]: 24 : do isppol=1,matlu(iatom)%nsppol
4719 [ + + ]: 104 : do im=1,ndim
4720 [ + + ]: 96 : if (isppol .eq. 1) then
4721 : 40 : matlu(iatom)%mat(im,im,isppol) = half*bfield
4722 : : else
4723 : 40 : matlu(iatom)%mat(im,im,isppol) = -half*bfield
4724 : : endif
4725 : : enddo ! im
4726 : : enddo ! isppol
4727 : : endif ! lpawu
4728 : : enddo ! natom
4729 : :
4730 : :
4731 [ # # ]: 0 : elseif(option .eq. 2) then
4732 : :
4733 : : !================================
4734 : : ! Spin-orbit magnetism
4735 : : ! H = mu_B*(L_z+g_e*S_Z)*B_z
4736 : : !================================
4737 : :
4738 [ # # ]: 0 : do iatom=1,natom
4739 [ # # ]: 0 : if(matlu(iatom)%lpawu .ne. -1) then
4740 : 0 : tndim=2*(2*matlu(iatom)%lpawu+1)
4741 : 0 : ll=matlu(iatom)%lpawu
4742 : :
4743 : 0 : jc1=0
4744 [ # # ]: 0 : do ms1=-1,1
4745 : 0 : xj=float(ms1)+half
4746 [ # # ]: 0 : do ml1=-ll,ll
4747 : 0 : jc1=jc1+1
4748 [ # # ]: 0 : if(jc1 < tndim+1) then
4749 [ # # ]: 0 : if (xj < 0.0) then
4750 : 0 : magnmatb(iatom)%value(jc1,jc1) = half*(ml1-2*xj)*bfield
4751 [ # # ]: 0 : elseif(xj > 0.0) then
4752 : 0 : magnmatb(iatom)%value(jc1,jc1) = half*(ml1-2*xj)*bfield
4753 : : endif
4754 : : endif
4755 : : enddo !ml1
4756 : : enddo ! ms1
4757 : : endif !lpawu
4758 : : enddo !natom
4759 : : endif !option
4760 : :
4761 : : !=======================
4762 : : ! reshape matrix
4763 : : !=======================
4764 : :
4765 [ - + ]: 8 : if(option .eq. 2) then
4766 : 0 : call gather_matlu(matlu,magnmatb(natom),natom,option=-1,prtopt=1)
4767 : : endif
4768 : :
4769 : : !================================
4770 : : ! Deallocate matrices
4771 : : !================================
4772 : :
4773 [ + + ]: 16 : do iatom=1,natom
4774 [ + - ]: 16 : if(matlu(iatom)%lpawu .ne. -1) then
4775 [ - + ]: 8 : ABI_FREE(magnmatb(iatom)%value)
4776 : : endif
4777 : : enddo
4778 : :
4779 [ + - + + : 16 : ABI_FREE(magnmatb)
- + - + ]
4780 : :
4781 : 8 : end subroutine magnfield_matlu
4782 : : !!***
4783 : :
4784 : :
4785 : : !!****f* m_matlu/magmomjmj_matlu
4786 : : !! NAME
4787 : : !! magmomjmj_matlu
4788 : : !!
4789 : : !! FUNCTION
4790 : : !! return the matrix of magnetic moments in the Jmj basis
4791 : : !!
4792 : : !!
4793 : : !! COPYRIGHT
4794 : : !! Copyright (C) 2005-2026 ABINIT group (FGendron)
4795 : : !! This file is distributed under the terms of the
4796 : : !! GNU General Public License, see ~abinit/COPYING
4797 : : !! or http://www.gnu.org/copyleft/gpl.txt .
4798 : : !!
4799 : : !! INPUTS
4800 : : !!
4801 : : !! OUTPUT
4802 : : !!
4803 : : !! SIDE EFFECTS
4804 : : !!
4805 : : !! NOTES
4806 : : !!
4807 : : !! SOURCE
4808 : 0 : subroutine magmomjmj_matlu(matlu,natom)
4809 : : use defs_basis
4810 : : use defs_wvltypes
4811 : : implicit none
4812 : :
4813 : : !Arguments ------------------------------------
4814 : : !scalars
4815 : : integer, intent(in) :: natom
4816 : : !arrays
4817 : : type(matlu_type), intent(inout) :: matlu(natom)
4818 : : !Local variables-------------------------------
4819 : : !scalars
4820 : : integer :: iatom,lpawu,ll,ml1,ms1,jm,jc1,tndim,jj
4821 : : real(dp) :: xj,xmj
4822 : : !arrays
4823 : 0 : integer, allocatable :: ind_msml(:,:)
4824 : 0 : type(coeff2c_type), allocatable :: gathermatlu(:)
4825 : 0 : complex(dpc),allocatable :: mlms2jmj(:,:)
4826 : : !************************************************************************
4827 : :
4828 : : !=====================================
4829 : : ! Allocate Matrices
4830 : : !=====================================
4831 : :
4832 [ # # # # : 0 : ABI_MALLOC(gathermatlu,(natom))
# # # # #
# # # #
# ]
4833 : :
4834 [ # # ]: 0 : do iatom=1,natom
4835 : 0 : lpawu=matlu(iatom)%lpawu
4836 [ # # ]: 0 : if(lpawu.ne.-1) then
4837 : 0 : ll=lpawu
4838 : 0 : tndim=2*(2*ll+1)
4839 : :
4840 [ # # # # : 0 : ABI_MALLOC(gathermatlu(iatom)%value,(tndim,tndim))
# # # # #
# # # # #
# # # # ]
4841 [ # # # # ]: 0 : gathermatlu(iatom)%value=czero
4842 [ # # # # : 0 : ABI_MALLOC(mlms2jmj,(tndim,tndim))
# # # # #
# # # #
# ]
4843 [ # # # # ]: 0 : mlms2jmj=czero
4844 [ # # # # : 0 : ABI_MALLOC(ind_msml,(2,-ll:ll))
# # # # #
# # # #
# ]
4845 [ # # # # ]: 0 : mlms2jmj=czero
4846 : :
4847 : : !=====================================
4848 : : ! Build J,M_J matrix
4849 : : !=====================================
4850 : :
4851 : : jc1=0
4852 [ # # ]: 0 : do ms1=1,2
4853 [ # # ]: 0 : do ml1=-ll,ll
4854 : 0 : jc1=jc1+1
4855 : 0 : ind_msml(ms1,ml1)=jc1
4856 : : end do
4857 : : end do
4858 : :
4859 : 0 : jc1=0
4860 [ # # ]: 0 : do jj=ll,ll+1
4861 : 0 : xj=float(jj)-half ! xj is in {ll-0.5, ll+0.5}
4862 [ # # ]: 0 : do jm=-jj,jj-1
4863 : 0 : xmj=float(jm)+half ! xmj is in {-xj,xj}
4864 : 0 : jc1=jc1+1 ! Global index for JMJ
4865 [ # # ]: 0 : if(nint(xj+0.5)==ll+1) then ! if xj=ll+0.5
4866 : 0 : mlms2jmj(jc1,jc1)=xmj ! J=L+0.5 and m_J=L+0.5
4867 [ # # ]: 0 : else if(nint(xj-0.5)==ll-1) then
4868 : 0 : mlms2jmj(jc1,jc1)=xmj ! J=L+0.5 and m_J=-L-0.5
4869 : : end if
4870 : : end do
4871 : : end do
4872 : :
4873 : : !print to debug
4874 : : !write(message,'(3a)') ch10,"JMJ Matrix"
4875 : : !call wrtout(std_out,message,"COLL")
4876 : : !do im=1,2*(ll*2+1)
4877 : : ! write(message,'(12(1x,18(1x,f5.2,f5.2)))') (mlms2jmj(im,jm),jm=1,2*(ll*2+1))
4878 : : ! call wrtout(std_out,message,"COLL")
4879 : : !end do
4880 : :
4881 : : !=====================================
4882 : : ! Put back into matlu format
4883 : : !=====================================
4884 : :
4885 [ # # # # : 0 : gathermatlu(iatom)%value=mlms2jmj
# # # # #
# # # ]
4886 : :
4887 : 0 : call gather_matlu(matlu,gathermatlu(iatom),natom=1,option=-1,prtopt=0)
4888 : :
4889 : : !=====================================
4890 : : ! Deallocate Matrices
4891 : : !=====================================
4892 : :
4893 [ # # ]: 0 : ABI_FREE(gathermatlu(iatom)%value)
4894 : : end if !lpawu
4895 : : end do !natom
4896 : :
4897 [ # # ]: 0 : ABI_FREE(mlms2jmj)
4898 [ # # ]: 0 : ABI_FREE(ind_msml)
4899 [ # # # # : 0 : ABI_FREE(gathermatlu)
# # # # ]
4900 : :
4901 : 0 : end subroutine magmomjmj_matlu
4902 : :
4903 : :
4904 [ # # # # : 0 : END MODULE m_matlu
# # # # #
# # # ]
4905 : : !!***
|
