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