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