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