Line data Source code
1 : !!****m* ABINIT/m_hu
2 : !! NAME
3 : !! m_hu
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 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 :
24 : #include "abi_common.h"
25 :
26 : MODULE m_hu
27 :
28 : use defs_basis
29 : use m_abicore
30 : use m_errors
31 :
32 : use m_abi_linalg, only : abi_xgemm
33 : use m_matlu, only : matlu_type
34 : use m_paw_dmft, only : paw_dmft_type
35 : use m_pawtab, only : pawtab_type
36 : use m_crystal, only : crystal_t
37 :
38 : implicit none
39 :
40 : private
41 :
42 : public :: init_vee
43 : public :: destroy_vee
44 : public :: init_hu
45 : public :: copy_hu
46 : public :: destroy_hu
47 : !public :: qmc_hu
48 : public :: print_hu
49 : public :: vee2udens_hu
50 : public :: rotatevee_hu
51 : public :: printvee_hu
52 : public :: vee2udensatom_hu
53 : public :: vee_slm2ylm_hu
54 : public :: vee_ndim2tndim_hu
55 : public :: vee_ndim2tndim_hu_r
56 : public :: udens_slatercondon_hu
57 : public :: udens_inglis_hu
58 :
59 : !!***
60 :
61 : !!****t* m_hu/vee_type
62 : !! NAME
63 : !! vee_type
64 : !!
65 : !! FUNCTION
66 : !! Structured datatype to store the U tensor
67 : !! for each atom.
68 : !!
69 : !! SOURCE
70 :
71 : type, public :: vee_type ! for each atom
72 :
73 : complex(dp), allocatable :: mat(:,:,:,:)
74 :
75 : end type vee_type
76 : !!***
77 :
78 : !----------------------------------------------------------------------
79 :
80 :
81 : !!****t* m_hu/hu_type
82 : !! NAME
83 : !! hu_type
84 : !!
85 : !! FUNCTION
86 : !! This structured datatype contains interaction matrices for the correlated subspace
87 : !!
88 : !! SOURCE
89 :
90 : type, public :: hu_type ! for each typat
91 :
92 : integer :: lpawu
93 :
94 : !logical :: jmjbasis
95 :
96 : real(dp) :: upawu ! => upaw
97 :
98 : real(dp) :: jpawu ! => jpaw
99 :
100 : logical :: jpawu_zero ! true if all jpawu are zero
101 : ! false if one of the jpaw is not zero
102 :
103 : real(dp), allocatable :: fk(:)
104 :
105 : complex(dp), allocatable :: udens(:,:)
106 :
107 : complex(dp), allocatable :: uqmc(:)
108 :
109 : complex(dp), allocatable :: vee(:,:,:,:)
110 :
111 : complex(dp), allocatable :: veeslm2(:,:,:,:)
112 :
113 : end type hu_type
114 :
115 : !----------------------------------------------------------------------
116 :
117 : CONTAINS !========================================================================================
118 : !!***
119 :
120 : !!****f* m_hu/init_vee
121 : !! NAME
122 : !! init_vee
123 : !!
124 : !! FUNCTION
125 : !! Allocate variables used in type vee_type.
126 : !!
127 : !! INPUTS
128 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
129 : !! vee = tensor for the interactions
130 : !!
131 : !! OUTPUTS
132 : !!
133 : !! SOURCE
134 :
135 103 : subroutine init_vee(paw_dmft,vee)
136 :
137 : !Arguments ------------------------------------
138 : type(paw_dmft_type), intent(in) :: paw_dmft
139 : type(vee_type), intent(inout) :: vee(paw_dmft%natom)
140 : !Local variables ------------------------------------
141 : integer :: iatom,lpawu,ndim
142 : !************************************************************************
143 :
144 467 : do iatom=1,paw_dmft%natom
145 364 : lpawu = paw_dmft%lpawu(iatom)
146 364 : if (lpawu == -1) cycle
147 127 : ndim = 2 * (2*lpawu+1)
148 762 : ABI_MALLOC(vee(iatom)%mat,(ndim,ndim,ndim,ndim))
149 1506985 : vee(iatom)%mat(:,:,:,:) = czero
150 : end do ! iatom
151 :
152 103 : end subroutine init_vee
153 : !!***
154 :
155 : !!****f* m_hu/destroy_vee
156 : !! NAME
157 : !! destroy_vee
158 : !!
159 : !! FUNCTION
160 : !! Destroy vee
161 : !!
162 : !! INPUTS
163 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
164 : !! vee = tensor for the interactions
165 : !!
166 : !! OUTPUTS
167 : !!
168 : !! SOURCE
169 :
170 103 : subroutine destroy_vee(paw_dmft,vee)
171 :
172 : !Arguments ------------------------------------
173 : type(paw_dmft_type), intent(in) :: paw_dmft
174 : type(vee_type), intent(inout) :: vee(paw_dmft%natom)
175 : !Local variables ------------------------------------
176 : integer :: iatom,lpawu
177 : !************************************************************************
178 :
179 467 : do iatom=1,paw_dmft%natom
180 364 : lpawu = paw_dmft%lpawu(iatom)
181 364 : if (lpawu == -1) cycle
182 467 : ABI_FREE(vee(iatom)%mat)
183 : end do ! iatom
184 :
185 103 : end subroutine destroy_vee
186 : !!***
187 :
188 : !!****f* m_hu/init_hu
189 : !! NAME
190 : !! init_hu
191 : !!
192 : !! FUNCTION
193 : !! Allocate variables used in type hu_type.
194 : !!
195 : !! INPUTS
196 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
197 : !! pawtab <type(pawtab)>=paw related data
198 : !!
199 : !! OUTPUTS
200 : !! hu <type(hu_type)>= U interaction
201 : !!
202 : !! SOURCE
203 :
204 93 : subroutine init_hu(hu,paw_dmft,pawtab)
205 :
206 : !Arguments ------------------------------------
207 : type(paw_dmft_type), intent(in) :: paw_dmft
208 : type(pawtab_type), intent(in) :: pawtab(paw_dmft%ntypat)
209 : type(hu_type), intent(inout) :: hu(paw_dmft%ntypat)
210 : !Local variables ------------------------------------
211 : integer :: i,ij,ij1,ij2,itypat,lpawu,m
212 : integer :: m1,ms,ms1,ndim,ntypat,tndim
213 : logical :: dmft_optim,t2g,x2my2d
214 : real(dp) :: jpawu,upawu,xtemp
215 : integer, parameter :: mt2g(3) = (/1,2,4/)
216 93 : integer, allocatable :: xij(:,:)
217 : character(len=4) :: tag
218 : character(len=500) :: message
219 : !************************************************************************
220 :
221 93 : ntypat = paw_dmft%ntypat
222 93 : t2g = (paw_dmft%dmft_t2g == 1)
223 93 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
224 :
225 93 : dmft_optim = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
226 :
227 93 : write(message,'(2a)') ch10," == Compute Interactions for DMFT"
228 93 : call wrtout(std_out,message,'COLL')
229 :
230 93 : xtemp = zero
231 :
232 : ! ====================================
233 : ! Compute hu(iatom)%uqmc from vee
234 : ! ====================================
235 93 : hu(1)%jpawu_zero = .true.
236 296 : do itypat=1,ntypat
237 203 : lpawu = pawtab(itypat)%lpawu
238 203 : hu(itypat)%upawu = zero
239 203 : hu(itypat)%jpawu = zero
240 : !hu(itypat)%jmjbasis = .false.
241 203 : if (t2g .and. lpawu == 2) lpawu = 1
242 203 : if (x2my2d .and. lpawu == 2) lpawu = 0
243 203 : hu(itypat)%lpawu = lpawu
244 203 : if (lpawu == -1) cycle
245 93 : ndim = 2*lpawu + 1
246 93 : tndim = 2 * ndim
247 93 : hu(itypat)%upawu = pawtab(itypat)%upawu
248 93 : hu(itypat)%jpawu = pawtab(itypat)%jpawu
249 :
250 93 : if (hu(itypat)%jpawu > tol4) hu(1)%jpawu_zero = .false.
251 : ! ndim1=2*hu(itypat)%lpawu+1
252 :
253 : ! allocate(hu(itypat)%vee(ndim,ndim,ndim,ndim))
254 :
255 558 : ABI_MALLOC(hu(itypat)%vee,(ndim,ndim,ndim,ndim))
256 :
257 : ! t2g case begin
258 93 : if (t2g) then
259 188 : do ms1=1,ndim
260 611 : do m1=1,ndim
261 1833 : do ms=1,ndim
262 5499 : do m=1,ndim
263 : hu(itypat)%vee(m,ms,m1,ms1) = &
264 5076 : & cmplx(pawtab(itypat)%vee(mt2g(m),mt2g(ms),mt2g(m1),mt2g(ms1)),zero,kind=dp)
265 : end do ! m
266 : end do ! ms
267 : end do ! m1
268 : end do ! ms1
269 : ! t2g case end
270 : ! x2my2d case begin
271 46 : else if (x2my2d) then
272 0 : hu(itypat)%vee(1,1,1,1) = cmplx(pawtab(itypat)%upawu,zero,kind=dp)
273 : ! x2my2d case end
274 : else
275 70266 : hu(itypat)%vee(:,:,:,:) = cmplx(pawtab(itypat)%vee(:,:,:,:),zero,kind=dp)
276 : end if ! t2g or xymy2d
277 : ! x2my2d case end
278 :
279 558 : ABI_MALLOC(hu(itypat)%veeslm2,(tndim,tndim,tndim,tndim))
280 93 : call vee_ndim2tndim_hu(lpawu,hu(itypat)%vee(:,:,:,:),hu(itypat)%veeslm2(:,:,:,:))
281 :
282 : ! This is copied from pawpuxinit: it would be better not to duplicate
283 : ! these lines.
284 279 : ABI_MALLOC(hu(itypat)%fk,(0:lpawu)) ! not used in the t2g and x2my2d cases
285 93 : hu(itypat)%fk(0) = hu(itypat)%upawu
286 93 : if (lpawu == 1) then
287 47 : hu(itypat)%fk(1) = hu(itypat)%jpawu * dble(5)
288 46 : else if (lpawu == 2) then
289 29 : hu(itypat)%fk(1) = hu(itypat)%jpawu * dble(14) / (one+pawtab(itypat)%f4of2_sla)
290 29 : hu(itypat)%fk(2) = hu(itypat)%fk(1) * pawtab(itypat)%f4of2_sla
291 17 : else if (lpawu == 3) then
292 : hu(itypat)%fk(1) = hu(itypat)%jpawu * dble(6435) / (dble(286)+&
293 17 : & dble(195)*pawtab(itypat)%f4of2_sla+dble(250)*pawtab(itypat)%f6of2_sla)
294 17 : hu(itypat)%fk(2) = hu(itypat)%fk(1) * pawtab(itypat)%f4of2_sla
295 17 : hu(itypat)%fk(3) = hu(itypat)%fk(1) * pawtab(itypat)%f6of2_sla
296 : end if ! lpawu
297 :
298 93 : write(tag,'(i4)') itypat
299 93 : write(message,'(3a)') ch10,' -------> For Correlated Species ',adjustl(tag)
300 93 : call wrtout(std_out,message,'COLL')
301 :
302 279 : ABI_MALLOC(hu(itypat)%uqmc,(ndim*(tndim-1)))
303 372 : ABI_MALLOC(hu(itypat)%udens,(tndim,tndim))
304 372 : ABI_MALLOC(xij,(tndim,tndim))
305 :
306 8827 : hu(itypat)%udens(:,:) = czero
307 : ij = 0
308 810 : do ms=1,tndim-1
309 717 : xij(ms,ms) = 0
310 717 : m = mod(ms-1,ndim) + 1
311 4367 : do ms1=ms+1,tndim
312 3557 : ij = ij + 1
313 3557 : xij(ms,ms1) = ij
314 3557 : xij(ms1,ms) = ij
315 3557 : m1 = mod(ms1-1,ndim) + 1
316 3557 : if ((ms <= ndim) .and. (ms1 > ndim)) then
317 1981 : hu(itypat)%uqmc(ij) = hu(itypat)%vee(m,m1,m,m1)
318 : else
319 1576 : hu(itypat)%uqmc(ij) = hu(itypat)%vee(m,m1,m,m1) - hu(itypat)%vee(m,m1,m1,m)
320 : end if
321 3557 : hu(itypat)%udens(ms,ms1) = hu(itypat)%uqmc(ij)
322 4274 : hu(itypat)%udens(ms1,ms) = hu(itypat)%udens(ms,ms1)
323 : end do ! ms1
324 : end do ! ms
325 :
326 93 : if (t2g .and. dmft_optim) then
327 : upawu = zero
328 : jpawu = zero
329 0 : do ms1=1,ndim
330 0 : do ms=1,ndim
331 0 : upawu = upawu + dble(hu(itypat)%vee(ms,ms1,ms,ms1))
332 0 : jpawu = jpawu + dble(hu(itypat)%vee(ms,ms1,ms,ms1)-hu(itypat)%vee(ms,ms1,ms1,ms))
333 : end do ! ms
334 : end do ! ms1
335 0 : upawu = upawu / dble(ndim**2)
336 0 : jpawu = upawu - jpawu/dble(2*lpawu*ndim)
337 0 : hu(itypat)%upawu = upawu
338 0 : hu(itypat)%jpawu = jpawu
339 : end if ! t2g and dmft_optim=1
340 :
341 93 : xij(tndim,tndim) = 0
342 93 : write(message,'(a,5x,a)') ch10,"-------- Interactions in the density-density representation, in the cubic basis "
343 93 : call wrtout(std_out,message,'COLL')
344 903 : write(message,'(1x,14(2x,i5))') (m,m=1,tndim)
345 93 : call wrtout(std_out,message,'COLL')
346 : ! xtemp1b=0.d0
347 : ! ====================================
348 : ! Print hu(iatom)%uqmc
349 : ! ====================================
350 93 : ij2 = 0
351 903 : do i=1,tndim
352 810 : if (i < tndim) then
353 717 : ij1 = ij2 + 1
354 717 : ij2 = ij2 + tndim - i
355 : end if
356 : ! write(std_out,*) itypat
357 : ! do m=1,i
358 : ! write(std_out,*) i,m
359 : ! write(std_out,*) xij(i,m)
360 : ! write(std_out,*) ij1,ij2
361 : ! enddo
362 1527 : if (i == 1) write(message,'(i3,14f7.3)') i,xtemp,(dble(hu(itypat)%uqmc(m)),m=ij1,ij2)
363 1434 : if (i /= tndim .and. i /= 1) write(message,'(i3,14f7.3)') i, &
364 6928 : & (dble(hu(itypat)%uqmc(xij(i,m))),m=1,i-1),xtemp,(dble(hu(itypat)%uqmc(m)),m=ij1,ij2)
365 1527 : if (i == tndim) write(message,'(i3,14f7.3)') i,(dble(hu(itypat)%uqmc(xij(i,m))),m=1,i-1),xtemp
366 903 : call wrtout(std_out,message,'COLL')
367 : end do ! i
368 93 : write(message,'(5x,a)') "--------------------------------------------------------"
369 93 : call wrtout(std_out,message,'COLL')
370 296 : ABI_FREE(xij)
371 :
372 : !imag part for debug
373 : ! ij2 = 0
374 : ! do i=1,tndim
375 : ! if (i < tndim) then
376 : ! ij1 = ij2 + 1
377 : ! ij2 = ij2 + tndim - i
378 : ! end if
379 : ! if (i == 1) write(message,'(i3,14f7.3)') i,xtemp,(aimag(hu(itypat)%uqmc(m)),m=ij1,ij2)
380 : ! if (i /= tndim .and. i /= 1) write(message,'(i3,14f7.3)') i, &
381 : ! & (aimag(hu(itypat)%uqmc(xij(i,m))),m=1,i-1),xtemp,(aimag(hu(itypat)%uqmc(m)),m=ij1,ij2)
382 : ! if (i == tndim) write(message,'(i3,14f7.3)') i,(aimag(hu(itypat)%uqmc(xij(i,m))),m=1,i-1),xtemp
383 : ! call wrtout(std_out,message,'COLL')
384 : ! end do ! i
385 : ! write(message,'(5x,a)') "--------------------------------------------------------"
386 : ! call wrtout(std_out,message,'COLL')
387 : ! ABI_FREE(xij)
388 :
389 :
390 : end do ! itypat
391 :
392 93 : end subroutine init_hu
393 : !!***
394 :
395 : !!****f* m_hu/copy_hu
396 : !! NAME
397 : !! copy_hu
398 : !!
399 : !! FUNCTION
400 : !! Copy hu into hu_new
401 : !!
402 : !! INPUTS
403 : !! hu <type(hu_type)>= U interaction
404 : !!
405 : !! OUTPUTS
406 : !! hu_new <type(hu_type)>= U interaction
407 : !!
408 : !! SOURCE
409 :
410 0 : subroutine copy_hu(ntypat,hu,hu_new)
411 :
412 : !Arguments ------------------------------------
413 : !type
414 : integer, intent(in) :: ntypat
415 : type(hu_type), intent(in) :: hu(ntypat)
416 : type(hu_type), intent(inout) :: hu_new(ntypat)
417 : !Local variables ------------------------------------
418 : integer :: itypat,ndim
419 : !************************************************************************
420 :
421 0 : do itypat=1,ntypat
422 0 : hu_new(itypat)%lpawu = hu(itypat)%lpawu
423 : !hu_new(itypat)%jmjbasis = hu(itypat)%jmjbasis
424 0 : hu_new(itypat)%upawu = hu(itypat)%upawu
425 0 : hu_new(itypat)%jpawu = hu(itypat)%jpawu
426 0 : hu_new(itypat)%jpawu_zero = hu(itypat)%jpawu_zero
427 0 : ndim=2*hu_new(itypat)%lpawu+1
428 0 : ABI_MALLOC(hu_new(itypat)%uqmc,(ndim*(2*ndim-1)))
429 0 : ABI_MALLOC(hu_new(itypat)%udens,(2*ndim,2*ndim))
430 0 : ABI_MALLOC(hu_new(itypat)%vee,(ndim,ndim,ndim,ndim))
431 0 : ABI_MALLOC(hu_new(itypat)%fk,(0:hu_new(itypat)%lpawu))
432 0 : hu_new(itypat)%vee = hu(itypat)%vee
433 0 : hu_new(itypat)%udens = hu(itypat)%udens
434 0 : hu_new(itypat)%uqmc = hu(itypat)%uqmc
435 0 : hu_new(itypat)%fk = hu(itypat)%fk
436 : end do ! itypat
437 :
438 0 : end subroutine copy_hu
439 : !!***
440 :
441 : !!****f* m_hu/destroy_hu
442 : !! NAME
443 : !! destroy_hu
444 : !!
445 : !! FUNCTION
446 : !! Deallocate hu
447 : !!
448 : !! INPUTS
449 : !! hu <type(hu_type)> = data for the interaction in DMFT.
450 : !! ntypat = number of species
451 : !!
452 : !! OUTPUT
453 : !!
454 : !! SOURCE
455 :
456 93 : subroutine destroy_hu(hu,ntypat)
457 :
458 : !Arguments ------------------------------------
459 : integer, intent(in) :: ntypat
460 : type(hu_type), intent(inout) :: hu(ntypat)
461 : !Local variables-------------------------------
462 : integer :: itypat
463 : ! *********************************************************************
464 :
465 296 : do itypat=1,ntypat
466 203 : ABI_SFREE(hu(itypat)%uqmc)
467 203 : ABI_SFREE(hu(itypat)%udens)
468 203 : ABI_SFREE(hu(itypat)%fk)
469 203 : ABI_SFREE(hu(itypat)%vee)
470 296 : ABI_SFREE(hu(itypat)%veeslm2)
471 : end do ! itypat
472 :
473 93 : end subroutine destroy_hu
474 : !!***
475 :
476 : !!****f* m_hu/print_hu
477 : !! NAME
478 : !! print_hu
479 : !!
480 : !! FUNCTION
481 : !! print density density interaction (used for DFT+DMFT)
482 : !!
483 : !! INPUTS
484 : !! ntypat = number of species
485 : !! prtopt = option for printing
486 : !! hu <type(hu_type)> = data for the interaction in DMFT.
487 : !!
488 : !! OUTPUT
489 : !!
490 : !! SOURCE
491 :
492 0 : subroutine print_hu(hu,ntypat,prtopt)
493 :
494 : !Arguments ------------------------------------
495 : !type
496 : integer, intent(in):: ntypat
497 : type(hu_type),intent(in) :: hu(ntypat)
498 : integer :: prtopt
499 :
500 : !Local variables-------------------------------
501 : integer :: itypat
502 : integer :: lpawu,ms,ms1,m,ndim
503 : character(len=500) :: message
504 : ! *********************************************************************
505 :
506 0 : do itypat = 1 , ntypat
507 0 : lpawu=hu(itypat)%lpawu
508 0 : if(lpawu/=-1) then
509 0 : ndim=2*lpawu+1
510 0 : write(message,'(2a,i4)') ch10,' -------> For Correlated species'
511 0 : call wrtout(std_out, message,'COLL')
512 0 : if(prtopt==0) then
513 0 : write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in cubic basis "
514 0 : else if(prtopt==1) then
515 0 : write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in diagonal basis"
516 0 : else if(prtopt==2) then
517 0 : write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in Ylm basis"
518 0 : else if(prtopt==3) then
519 0 : write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in JMJ basis"
520 : endif
521 0 : call wrtout(std_out, message,'COLL')
522 0 : write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
523 0 : call wrtout(std_out, message,'COLL')
524 0 : do ms=1,2*ndim
525 : write(message,'(i3,14f7.3)') &
526 0 : & ms, (dble(hu(itypat)%udens(ms,ms1)),ms1=1,2*ndim)
527 0 : call wrtout(std_out, message,'COLL')
528 : enddo
529 0 : write(message,'(5x,a)') "--------------------------------------------------------"
530 0 : call wrtout(std_out, message,'COLL')
531 : endif ! lpawu/=1
532 : enddo ! ntypat
533 :
534 :
535 0 : end subroutine print_hu
536 : !!***
537 :
538 : !!****f* m_hu/vee2udens_hu
539 : !! NAME
540 : !! print_hu
541 : !!
542 : !! FUNCTION
543 : !! interaction udens in recomputed from new vee.
544 : !!
545 : !! INPUTS
546 : !! ntypat = number of species
547 : !! prtopt = option for printing
548 : !!
549 : !! OUTPUT
550 : !!
551 : !! SIDE EFFECT
552 : !! hu <type(hu_type)> = data for the interaction in DMFT.
553 : !!
554 : !! SOURCE
555 :
556 0 : subroutine vee2udens_hu(hu,ntypat,prtopt)
557 :
558 : !Arguments ------------------------------------
559 : !type
560 : integer, intent(in):: ntypat
561 : type(hu_type),intent(inout) :: hu(ntypat)
562 : integer :: prtopt
563 :
564 : !Local variables-------------------------------
565 : integer :: ij,itypat
566 : integer :: lpawu,m1,ms,ms1,m,ndim
567 : character(len=500) :: message
568 : ! *********************************************************************
569 0 : do itypat=1,ntypat
570 0 : lpawu=hu(itypat)%lpawu
571 0 : if(lpawu.ne.-1) then
572 0 : ndim=2*lpawu+1
573 0 : write(message,'(2a,i4)') ch10,' -------> For Correlated Species', itypat
574 0 : call wrtout(std_out, message,'COLL')
575 :
576 0 : hu(itypat)%udens=zero
577 0 : ij=0
578 0 : do ms=1,2*ndim-1
579 : ! xij(ms,ms)=0
580 0 : do ms1=ms+1,2*ndim
581 0 : ij=ij+1
582 : ! xij(ms,ms1)=ij
583 : ! xij(ms1,ms)=ij
584 0 : if(ms<=ndim.and.ms1>ndim) then
585 0 : m1 = ms1 - ndim
586 0 : m = ms
587 0 : hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)
588 0 : hu(itypat)%udens(ms,ms1)= hu(itypat)%vee(m,m1,m,m1)
589 0 : hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
590 0 : else if(ms<=ndim.and.ms1<=ndim) then
591 0 : m1 = ms1
592 0 : m = ms
593 0 : hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)-hu(itypat)%vee(m,m1,m1,m)
594 0 : hu(itypat)%udens(ms,ms1)= hu(itypat)%uqmc(ij)
595 0 : hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
596 : else
597 0 : m1 = ms1 - ndim
598 0 : m = ms - ndim
599 0 : hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)-hu(itypat)%vee(m,m1,m1,m)
600 0 : hu(itypat)%udens(ms,ms1)= hu(itypat)%uqmc(ij)
601 0 : hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
602 : endif
603 : enddo
604 : enddo
605 : ! xij(2*ndim,2*ndim)=0
606 : ! write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation "
607 : ! call wrtout(std_out, message,'COLL')
608 : ! write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
609 : ! call wrtout(std_out, message,'COLL')
610 : endif
611 : enddo ! itypat
612 0 : call print_hu(hu,ntypat,prtopt)
613 :
614 :
615 0 : end subroutine vee2udens_hu
616 : !!***
617 :
618 : !!****f* m_hu/rotatevee_hu
619 : !! NAME
620 : !! rotatevee_hu
621 : !!
622 : !! FUNCTION
623 : !! Rotate U matrix to new basis
624 : !!
625 : !! INPUTS
626 : !! hu <type(hu_type)>= U interaction
627 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
628 : !! pawprtvol >= 3 : print different quantities
629 : !! rot_mat = rotation matrix
630 : !! rot_type = 0 ! keep original Slm basis
631 : !! = 1 ! use the rotation matrix rot_mat from diago of dmat, green, levels..
632 : !! = 2 ! rotation to the Ylm basis
633 : !! = 3 ! rotation to the JmJ Basis
634 : !! = 4 ! same as 1 but rot_mat is applied from the Ylm basis instead of the Slm basis
635 : !!
636 : !! OUTPUT
637 : !! udens_atoms = rotated udens
638 : !! vee_rotated = rotated vee
639 : !!
640 : !! SIDE EFFECT
641 : !! hu <type(hu_type)> = data for the interaction in DMFT.
642 : !!
643 : !! SOURCE
644 :
645 83 : subroutine rotatevee_hu(hu,paw_dmft,pawprtvol,rot_mat,rot_type,udens_atoms,vee_rotated)
646 :
647 : !Arguments ------------------------------------
648 : type(paw_dmft_type), intent(in) :: paw_dmft
649 : integer, intent(in) :: pawprtvol,rot_type
650 : type(hu_type), target, intent(inout) :: hu(paw_dmft%ntypat)
651 : type(matlu_type), intent(in) :: rot_mat(paw_dmft%natom)
652 : type(matlu_type), intent(inout) :: udens_atoms(paw_dmft%natom)
653 : type(vee_type), target, intent(inout) :: vee_rotated(paw_dmft%natom)
654 : !Local variables-------------------------------
655 : integer :: iatom,itypat,lpawu,m,m1,m2,mi,ms,ms1,nat_correl
656 : integer :: natom,ndim,nflavor,nspinor,nsppol,nsppol_,prtonly,tndim,prtopt
657 : logical :: triqs
658 : real(dp) :: f2,jpawu,xsum,xsum2,xsum2new,xsumnew
659 : character(len=4) :: tag_at
660 : character(len=30) :: basis_vee
661 : character(len=500) :: message
662 83 : complex(dp), target, allocatable :: veeylm(:,:,:,:)
663 : complex(dp), pointer :: veeslm(:,:,:,:) => null(),veetemp(:,:,:,:) => null()
664 : complex(dp), pointer :: veetemp2(:,:,:,:) => null(),veetemp3(:,:,:,:) => null()
665 : complex(dp), pointer :: veeylm2(:,:,:,:) => null()
666 : ! *********************************************************************
667 :
668 83 : natom = paw_dmft%natom
669 83 : nspinor = paw_dmft%nspinor
670 83 : nsppol = paw_dmft%nsppol
671 :
672 83 : triqs = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7 .or. paw_dmft%dmft_solv == 10)
673 :
674 83 : if(paw_dmft%dmft_solv .eq. 10) then
675 0 : prtopt = 4
676 : else
677 83 : prtopt = 1
678 : endif
679 :
680 83 : write(message,'(a,3x,a)') ch10,"== Rotate interaction to the CTQMC basis"
681 83 : call wrtout(std_out,message,"COLL")
682 :
683 : !================================================
684 : ! NSPINOR = 2
685 : !================================================
686 :
687 83 : if (nspinor == 2) then
688 :
689 4 : do iatom=1,natom
690 2 : lpawu = paw_dmft%lpawu(iatom)
691 2 : if (lpawu == -1) cycle
692 2 : itypat = paw_dmft%typat(iatom)
693 2 : jpawu = hu(itypat)%jpawu
694 2 : ndim = 2*lpawu + 1
695 2 : tndim = nspinor * ndim
696 : !if(pawprtvol>=3) then
697 : ! write(message,'(2a)') ch10," VEE INPUT AVANT TRANSFORMATION"
698 : ! call wrtout(std_out, message,'COLL')
699 : ! call printvee_hu(ndim,hu(itypat)%vee,4,'Slm')
700 : !endif
701 :
702 2 : write(tag_at,'(i4)') iatom
703 2 : write(message,'(3a)') ch10,' -------> For Correlated atom ',adjustl(tag_at)
704 2 : call wrtout(std_out,message,'COLL')
705 :
706 : ! ==================================
707 : ! First print veeslm
708 : ! ==================================
709 :
710 : ! Print udens in the Slm basis
711 2 : call vee2udensatom_hu(ndim,hu(itypat)%udens(:,:),hu(itypat)%vee(:,:,:,:),"cubic",prtonly=1)
712 :
713 2 : if (.not. paw_dmft%dmft_solv .eq. 10) then
714 2 : basis_vee = "cubic"
715 : ! First print veeslm
716 : !call printvee_hu(ndim,real(veeslm),1,basis_vee)
717 2 : call printvee_hu(tndim,hu(itypat)%veeslm2(:,:,:,:),prtopt,basis_vee)
718 : end if
719 :
720 : ! ==================================
721 : ! Then compute veerotated
722 : ! ==================================
723 :
724 : ! In the basis where levels/density matrix/green function is diagonal
725 : ! ================================================================================
726 :
727 : ! When J=0, vee(:,i,:,j) and vee(i,:,j,:) are homotheties for all i,j
728 : ! when using Slater parametrization ; so no need to rotate them
729 2 : if (rot_type == 1 .and. jpawu > tol10) then
730 : ! ---------------------
731 :
732 : !veerotated=czero
733 : !do m1=1,tndim
734 : ! do m2=1,tndim
735 : ! do m3=1,tndim
736 : ! do m4=1,tndim
737 : ! do mi=1,tndim
738 : ! do mj=1,tndim
739 : ! do mk=1,tndim
740 : ! do ml=1,tndim
741 : ! veerotated(m1,m2,m3,m4)= veerotated(m1,m2,m3,m4) + &
742 : !& conjg(rot_mat(iatom,1)%value(mi,m1))* &
743 : !& conjg(rot_mat(iatom,1)%value(mj,m2))* &
744 : !& rot_mat(iatom,1)%value(mk,m3)* &
745 : !& rot_mat(iatom,1)%value(ml,m4)* &
746 : !& veeslm2(mi,mj,mk,ml)
747 : ! enddo
748 : ! enddo
749 : ! enddo
750 : ! enddo
751 : ! enddo
752 : ! enddo
753 : ! enddo
754 : ! enddo
755 :
756 0 : call rotate_hu(rot_mat(iatom)%mat(:,:,:),1,tndim,hu(itypat)%veeslm2(:,:,:,:),vee_rotated(iatom)%mat(:,:,:,:))
757 0 : basis_vee = "CTQMC basis from cubic"
758 :
759 : ! In the Ylm basis
760 : ! ================================================================================
761 2 : else if ((rot_type == 2 .or. rot_type == 3 .or. rot_type == 4) .and. jpawu > tol10) then
762 : ! ---------------------------
763 :
764 0 : ABI_MALLOC(veeylm,(ndim,ndim,ndim,ndim))
765 0 : if (rot_type == 2) then
766 0 : veeylm2 => vee_rotated(iatom)%mat(:,:,:,:)
767 : else
768 0 : ABI_MALLOC(veeylm2,(tndim,tndim,tndim,tndim))
769 : end if
770 : ! Change basis from slm to ylm basis
771 : !if (triqs) then
772 0 : veeslm => hu(itypat)%vee(:,:,:,:)
773 : !else
774 : ! ABI_MALLOC(veeslm,(ndim,ndim,ndim,ndim))
775 : ! veeslm(:,:,:,:) = cmplx(real(hu(itypat)%vee(:,:,:,:)),zero,kind=sp)
776 : !end if
777 :
778 0 : call vee_slm2ylm_hu(lpawu,veeslm(:,:,:,:),veeylm(:,:,:,:),paw_dmft,1,2)
779 :
780 : !if (.not. triqs) then
781 : ! ABI_FREE(veeslm)
782 : !end if
783 :
784 0 : veeslm => null()
785 :
786 : ! The line below is not really useful
787 : !if (.not. triqs) veeylm(:,:,:,:) = cmplx(dble(veeylm(:,:,:,:)),zero,kind=dp)
788 :
789 0 : basis_vee = 'Ylm'
790 0 : if (.not. paw_dmft%dmft_solv .eq. 10) then
791 : ! Print interaction matrix in the ylm basis
792 0 : call printvee_hu(ndim,veeylm(:,:,:,:),prtopt,basis_vee,hu(itypat)%upawu)
793 :
794 : ! Print interaction matrix in the ylm basis from Slater tables
795 0 : if (pawprtvol >= 4) then
796 0 : call udens_slatercondon_hu(hu(itypat)%fk(:),lpawu)
797 : end if
798 : end if
799 :
800 : ! Build large matrix
801 0 : call vee_ndim2tndim_hu(lpawu,veeylm(:,:,:,:),veeylm2(:,:,:,:))
802 :
803 0 : if(.not. paw_dmft%dmft_solv .eq. 10) then
804 0 : if (rot_type == 3 .or. rot_type == 4) then
805 0 : call printvee_hu(tndim,veeylm2(:,:,:,:),prtopt,basis_vee)
806 : end if
807 : endif
808 :
809 : ! ---------------------------
810 : !
811 : ! In the JmJ basis
812 : ! ================================================================================
813 0 : if (rot_type == 3) then
814 :
815 : ! apply change of basis
816 0 : call vee_ylm2jmj_hu(lpawu,veeylm2(:,:,:,:),vee_rotated(iatom)%mat(:,:,:,:),1,paw_dmft)
817 :
818 : ! print interaction matrix in the JMJ basis from Inglis and Julien tables
819 0 : if (pawprtvol >= 3) then
820 0 : call udens_inglis_hu(hu(itypat)%fk(:),lpawu)
821 : end if
822 :
823 : ! new dimension
824 :
825 0 : basis_vee = "JmJ"
826 :
827 0 : else if (rot_type == 4) then
828 :
829 : !veerotated=czero
830 :
831 : !do m1=1,tndim
832 : ! do m2=1,tndim
833 : ! do m3=1,tndim
834 : ! do m4=1,tndim
835 : ! do mi=1,tndim
836 : ! do mj=1,tndim
837 : ! do mk=1,tndim
838 : ! do ml=1,tndim
839 : ! veerotated(m1,m2,m3,m4)= veerotated(m1,m2,m3,m4) +
840 : ! &
841 : !& conjg(rot_mat(iatom,1)%value(mi,m1))* &
842 : !& conjg(rot_mat(iatom,1)%value(mj,m2))* &
843 : !& rot_mat(iatom,1)%value(mk,m3)* &
844 : !& rot_mat(iatom,1)%value(ml,m4)* &
845 : !& veeylm2(mi,mj,mk,ml)
846 : ! enddo
847 : ! enddo
848 : ! enddo
849 : ! enddo
850 : ! enddo
851 : ! enddo
852 : ! enddo
853 : !enddo
854 0 : if (.not. paw_dmft%dmft_solv .eq. 10) then
855 0 : call udens_inglis_hu(hu(itypat)%fk(:),lpawu)
856 : endif
857 0 : call rotate_hu(rot_mat(iatom)%mat(:,:,:),1,tndim,veeylm2(:,:,:,:),vee_rotated(iatom)%mat(:,:,:,:))
858 0 : basis_vee = 'CTQMC basis from Ylm'
859 :
860 : end if ! rot_type
861 :
862 : end if ! rot_type
863 :
864 2 : if (rot_type == 0 .or. (jpawu <= tol10)) then
865 82742 : vee_rotated(iatom)%mat(:,:,:,:) = hu(itypat)%veeslm2(:,:,:,:)
866 422 : udens_atoms(iatom)%mat(:,:,1) = hu(itypat)%udens(:,:)
867 : end if ! rot_type
868 :
869 2 : ABI_SFREE(veeylm)
870 2 : if (rot_type >= 3 .and. jpawu > tol10) then
871 0 : ABI_FREE(veeylm2)
872 : end if
873 2 : veeylm2 => null()
874 :
875 2 : f2 = zero
876 2 : if (lpawu /= 0) f2 = hu(itypat)%fk(1)
877 :
878 2 : if (.not. paw_dmft%dmft_solv .eq. 10) then
879 2 : call printvee_hu(tndim,vee_rotated(iatom)%mat(:,:,:,:),prtopt,basis_vee,hu(itypat)%upawu,f2)
880 : ! call printvee_hu(dim_vee,real(veeylm),1,hu(itypat)%upawu)
881 : endif
882 :
883 : !uaver=zero
884 2 : if (rot_type /= 0 .and. jpawu > tol10) then
885 : ! Careful, since vee now has spin off-diagonal elements, this is wrong to
886 : ! use vee2udensatom_hu in order to compute udens
887 0 : do ms1=1,tndim
888 0 : do ms=1,tndim
889 0 : udens_atoms(iatom)%mat(ms,ms1,1) = vee_rotated(iatom)%mat(ms,ms1,ms,ms1) - vee_rotated(iatom)%mat(ms,ms1,ms1,ms)
890 : !uaver=uaver+udens_atoms(iatom)%value(ms,ms1)
891 : end do ! ms
892 : end do ! ms1
893 : end if ! rot_type /= 0 and jpawu /= zero
894 :
895 4 : if (paw_dmft%dmft_solv .eq. 10) then
896 0 : call vee2udensatom_hu(ndim,udens_atoms(iatom)%mat(:,:,1),vee_rotated(iatom)%mat(:,:,:,:),"CTQMC",prtonly=2)
897 :
898 : !xsum=0
899 : !xsum2=0
900 : !do m1=1,ndim
901 : ! do m2=1,ndim
902 : ! xsum = xsum + aimag(vee_rotated(iatom)%mat(m1,m2,m1,m2))
903 : ! xsum2 = xsum2 + aimag(udens_atoms(iatom)%mat(m1,m2,1))
904 : ! end do ! m2
905 : !end do ! m1
906 :
907 : !write(*,*)"Sum of Imag part of Rotated Interaction Vee",xsum
908 : !write(message,'(5a)') " == Sum of Imag part of Rotated Interaction Udens",xsum2
909 : !call wrtout(std_out,message,'COLL')
910 :
911 : else
912 2 : call vee2udensatom_hu(ndim,udens_atoms(iatom)%mat(:,:,1),vee_rotated(iatom)%mat(:,:,:,:),basis_vee,prtonly=1)
913 : endif
914 :
915 : end do ! iatom
916 : !ABI_ERROR("Aborting now!")
917 :
918 : !================================================
919 : ! NSPINOR = 1
920 : !================================================
921 :
922 81 : else if (nspinor == 1) then
923 :
924 81 : nat_correl = 0
925 395 : do iatom=1,natom
926 314 : lpawu = paw_dmft%lpawu(iatom)
927 314 : if (lpawu == -1) cycle
928 105 : itypat = paw_dmft%typat(iatom)
929 105 : jpawu = hu(itypat)%jpawu
930 105 : nat_correl = nat_correl + 1
931 105 : if (nat_correl > 1 .and. (hu(itypat)%jpawu > tol4)) then
932 0 : write(message,'(3a)') ch10,' -------> Warning: several atoms: ',' not extensively tested '
933 0 : ABI_WARNING(message)
934 : end if
935 :
936 105 : write(tag_at,'(i4)') iatom
937 105 : write(message,'(3a)') ch10,' -------> For Correlated atom ',adjustl(tag_at)
938 105 : call wrtout(std_out,message,'COLL')
939 :
940 : ! ! ================================================================
941 : ! ! If rotation for spin 2 and rotation for spin 1 are not equal
942 : ! ! then print a warning
943 : ! ! useful only for magnetic case
944 : ! ! ================================================================
945 105 : ndim = 2*lpawu + 1
946 105 : tndim = nspinor * ndim
947 105 : nflavor = 2 * ndim
948 105 : if (nsppol == 2 .and. pawprtvol >= 3 .and. (.not. triqs)) then
949 0 : do m2=1,tndim
950 0 : do m1=1,tndim
951 0 : if (abs(rot_mat(iatom)%mat(m1,m2,1)-rot_mat(iatom)%mat(m1,m2,2)) > tol4) then
952 0 : write(message,'(2a,i4)') ch10,' rot_mat differs for value of isppol but value for isppol=2 not used'
953 0 : call wrtout(std_out,message,'COLL')
954 0 : write(message,'(a,4e16.8)') ch10,rot_mat(iatom)%mat(m1,m2,1),rot_mat(iatom)%mat(m1,m2,2)
955 0 : call wrtout(std_out,message,'COLL',do_flush=.True.)
956 : end if
957 : end do ! m1
958 : end do ! m2
959 : end if ! nsppol=2 and pawprtvol>=3
960 :
961 : ! ! =================================================
962 : ! ! See if rotation is complex or real
963 : ! ! =================================================
964 105 : if (pawprtvol >= 3 .and. (.not. triqs)) then
965 272 : do m1=1,ndim
966 1318 : do mi=1,ndim
967 1264 : if (abs(aimag(rot_mat(iatom)%mat(mi,m1,1))) > tol8) then
968 163 : write(message,'(2a,2i6,2e14.3)') ch10,"rot_mat is complex for", &
969 326 : & mi,m1,rot_mat(iatom)%mat(mi,m1,1)
970 163 : call wrtout(std_out,message,'COLL')
971 : end if
972 : end do ! mi
973 : end do ! m1
974 : end if ! pawprtvol
975 :
976 : ! write vee for information with a classification.
977 105 : if (pawprtvol >= 3) then
978 54 : call printvee_hu(ndim,hu(itypat)%vee(:,:,:,:),2,'cubic')
979 : end if
980 :
981 105 : basis_vee = 'cubic'
982 105 : prtonly = 1
983 105 : if (jpawu > tol10 .and. rot_type /= 0) then
984 :
985 : ! ! Compute rotated vee.
986 : !veetemp=zero
987 : !do m1=1,ndim
988 : ! do m2=1,ndim
989 : ! do m3=1,ndim
990 : ! do m4=1,ndim
991 : ! do mi=1,ndim
992 : ! do mj=1,ndim
993 : ! do mk=1,ndim
994 : ! do ml=1,ndim
995 : ! ! if((mi==mk.and.mj==ml).or.(mi==ml.and.mj==mk)) then
996 : ! veetemp(m1,m2,m3,m4)= veetemp(m1,m2,m3,m4) + &
997 : !& real( &
998 : !& conjg(rot_mat(iatom,1)%value(mi,m1))* &
999 : !& conjg(rot_mat(iatom,1)%value(mj,m2))* &
1000 : !& rot_mat(iatom,1)%value(mk,m3)* &
1001 : !& rot_mat(iatom,1)%value(ml,m4)* &
1002 : !& hu(itypat)%vee(mi,mj,mk,ml)&
1003 : ! )
1004 : ! endif
1005 : ! enddo
1006 : ! enddo
1007 : ! enddo
1008 : ! enddo
1009 : ! enddo
1010 : ! enddo
1011 : ! enddo
1012 : ! enddo
1013 :
1014 40 : nsppol_ = 1
1015 40 : prtonly = 0
1016 :
1017 : ! Use different rotation matrices for each spin with TRIQS
1018 40 : if (triqs .and. rot_type /= 2) nsppol_ = nsppol
1019 :
1020 40 : if (rot_type == 2 .or. rot_type == 4) then
1021 0 : basis_vee = 'Ylm'
1022 0 : ABI_MALLOC(veeylm,(ndim,ndim,ndim,ndim))
1023 0 : call vee_slm2ylm_hu(lpawu,hu(itypat)%vee(:,:,:,:),veeylm(:,:,:,:),paw_dmft,1,2)
1024 0 : if (rot_type == 2 .or. nsppol_ == 2) then
1025 0 : if (rot_type == 2) then
1026 0 : veeylm2 => vee_rotated(iatom)%mat(:,:,:,:)
1027 0 : else if (nsppol_ == 2) then
1028 0 : ABI_MALLOC(veeylm2,(nflavor,nflavor,nflavor,nflavor))
1029 : end if
1030 0 : call vee_ndim2tndim_hu(lpawu,veeylm(:,:,:,:),veeylm2(:,:,:,:))
1031 : end if ! rot_type
1032 0 : veetemp3 => veeylm(:,:,:,:)
1033 : end if ! rot_type=2 or 4
1034 :
1035 40 : if (rot_type == 1 .or. rot_type == 4) then
1036 40 : basis_vee = 'CTQMC basis from cubic'
1037 40 : if (rot_type == 4) basis_vee = 'CTQMC basis from Ylm'
1038 40 : if (nsppol_ == 2) then
1039 0 : if (rot_type == 1) then
1040 0 : veetemp => hu(itypat)%veeslm2(:,:,:,:)
1041 : else
1042 0 : veetemp => veeylm2
1043 : end if ! rot_type
1044 0 : veetemp2 => vee_rotated(iatom)%mat(:,:,:,:)
1045 : else
1046 40 : if (rot_type == 1) then
1047 40 : veetemp => hu(itypat)%vee(:,:,:,:)
1048 : else
1049 0 : veetemp => veeylm(:,:,:,:)
1050 : end if ! rot_type
1051 240 : ABI_MALLOC(veetemp2,(ndim,ndim,ndim,ndim))
1052 : end if ! nsppol_
1053 40 : call rotate_hu(rot_mat(iatom)%mat(:,:,:),nsppol_,ndim,veetemp(:,:,:,:),veetemp2(:,:,:,:))
1054 4840 : if (.not. triqs) veetemp2(:,:,:,:) = cmplx(dble(veetemp2(:,:,:,:)),zero,kind=dp) ! neglect imaginary part in Abinit
1055 40 : if (nsppol_ == 1) then
1056 40 : call vee_ndim2tndim_hu(lpawu,veetemp2(:,:,:,:),vee_rotated(iatom)%mat(:,:,:,:))
1057 40 : veetemp3 => veetemp2(:,:,:,:)
1058 : end if
1059 : end if ! rot_type=1 or 4
1060 :
1061 40 : if (nsppol_ == 2) then
1062 0 : prtonly = 1
1063 : ! It is wrong to use vee2udensatom_hu to build udens_atoms here
1064 0 : do m=1,nflavor
1065 0 : do m1=1,nflavor
1066 0 : udens_atoms(iatom)%mat(m,m1,1) = vee_rotated(iatom)%mat(m,m1,m,m1) - vee_rotated(iatom)%mat(m,m1,m1,m)
1067 : end do ! m1
1068 : end do ! m
1069 : end if ! nsppol_=2
1070 : else
1071 : prtonly = 1
1072 8815 : udens_atoms(iatom)%mat(:,:,1) = hu(itypat)%udens(:,:)
1073 1206375 : vee_rotated(iatom)%mat(:,:,:,:) = hu(itypat)%veeslm2(:,:,:,:)
1074 : end if ! jpawu=zero
1075 :
1076 105 : xsum = zero
1077 105 : xsum2 = zero
1078 105 : xsumnew = zero
1079 105 : xsum2new = zero
1080 582 : do m1=1,ndim
1081 2951 : do m2=1,ndim
1082 2369 : xsum = xsum + dble(hu(itypat)%vee(m1,m2,m1,m2))
1083 2369 : xsum2 = xsum2 + dble(hu(itypat)%vee(m1,m2,m2,m1))
1084 2369 : xsumnew = xsumnew + dble(vee_rotated(iatom)%mat(m1,m2,m1,m2))
1085 2846 : xsum2new = xsum2new + dble(vee_rotated(iatom)%mat(m1,m2,m2,m1))
1086 : end do ! m2
1087 : end do ! m1
1088 105 : if (abs(xsum-xsumnew) > tol5 .or. abs(xsum2-xsum2new) > tol5) then
1089 0 : write(message,'(2a)') ch10," BUG: New interaction after rotation do not respect sum rules"
1090 0 : call wrtout(std_out,message,'COLL')
1091 0 : write(message,'(2a,2f14.3)') ch10,' Comparison of \sum_{m1,m3} vee(m1,m3,m1,m3) before and after rotation is',&
1092 0 : & xsum,xsumnew
1093 0 : call wrtout(std_out,message,'COLL')
1094 0 : write(message,'(2a,2f14.3)') ch10,' Comparison of \sum_{m1,m3} vee(m1,m3,m3,m1) before and after rotation is',&
1095 0 : & xsum2,xsum2new
1096 0 : call wrtout(std_out,message,'COLL')
1097 : end if ! abs(xsum-xsumnew)>tol5
1098 105 : if (pawprtvol >= 3) then
1099 54 : write(message,'(2a)') ch10," VEE ROTATED"
1100 54 : call wrtout(std_out,message,'COLL')
1101 54 : call printvee_hu(2*ndim,vee_rotated(iatom)%mat(:,:,:,:),2,'CTQMC')
1102 54 : write(message,'(a)') ch10
1103 54 : call wrtout(std_out,message,'COLL')
1104 : end if ! pawprtvol>=3
1105 :
1106 105 : call vee2udensatom_hu(ndim,udens_atoms(iatom)%mat(:,:,1),veetemp3(:,:,:,:),basis_vee,prtonly=prtonly)
1107 :
1108 105 : veetemp => null()
1109 105 : veetemp3 => null()
1110 :
1111 105 : ABI_SFREE(veeylm)
1112 :
1113 105 : if (jpawu > tol10 .and. rot_type /= 0) then
1114 40 : if (rot_type == 4 .and. nsppol_ == 2) then
1115 0 : ABI_FREE(veeylm2)
1116 : end if
1117 40 : if (rot_type /= 2 .and. nsppol_ == 1) then
1118 40 : ABI_FREE(veetemp2)
1119 : end if
1120 : end if ! jpawu>tol10 and rot_type/=0
1121 :
1122 105 : veetemp2 => null()
1123 395 : veeylm2 => null()
1124 :
1125 : ! udens_atoms(iatom)%value=zero
1126 : ! ij=0
1127 : ! do ms=1,2*ndim-1
1128 : ! do ms1=ms+1,2*ndim
1129 : ! ij=ij+1
1130 : ! if(ms<=ndim.and.ms1>ndim) then
1131 : ! m1 = ms1 - ndim
1132 : ! m = ms
1133 : ! hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)
1134 : ! udens_atoms(iatom)%value(ms,ms1)= veetemp(m,m1,m,m1)
1135 : ! udens_atoms(iatom)%value(ms1,ms)= udens_atoms(iatom)%value(ms,ms1)
1136 : ! else if(ms<=ndim.and.ms1<=ndim) then
1137 : ! m1 = ms1
1138 : ! m = ms
1139 : ! hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1140 : ! udens_atoms(iatom)%value(ms,ms1)= hu(itypat)%uqmc(ij)
1141 : ! udens_atoms(iatom)%value(ms1,ms)= udens_atoms(iatom)%value(ms,ms1)
1142 : ! else
1143 : ! m1 = ms1 - ndim
1144 : ! m = ms - ndim
1145 : ! hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1146 : ! udens_atoms(iatom)%value(ms,ms1)= hu(itypat)%uqmc(ij)
1147 : ! udens_atoms(iatom)%value(ms1,ms)= udens_atoms(iatom)%value(ms,ms1)
1148 : ! endif
1149 : ! enddo
1150 : ! enddo
1151 : ! write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation "
1152 : ! call wrtout(std_out, message,'COLL')
1153 : ! write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
1154 : ! call wrtout(std_out, message,'COLL')
1155 : ! do ms=1,2*ndim
1156 : ! write(message,'(i3,14f7.3)') &
1157 : ! & ms, (udens_atoms(iatom)%value(ms,ms1),ms1=1,2*ndim)
1158 : ! call wrtout(std_out, message,'COLL')
1159 : ! enddo
1160 : ! write(message,'(5x,a)') "--------------------------------------------------------"
1161 : ! call wrtout(std_out, message,'COLL')
1162 : !ABI_FREE(veetemp)
1163 : ! endif ! lpawu/=1
1164 : ! call print_hu(hu,cryst_struc%ntypat,1)
1165 :
1166 : end do ! iatom
1167 : ! call print_hu(hu,cryst_struc%ntypat,1)
1168 : ! call vee2udens_hu(hu,cryst_struc%ntypat,2)
1169 : end if ! nspinor
1170 :
1171 83 : end subroutine rotatevee_hu
1172 : !!***
1173 :
1174 : !!****f* m_hu/rotate_hu
1175 : !! NAME
1176 : !! rotate_hu
1177 : !!
1178 : !! FUNCTION
1179 : !! Rotate an interaction tensor
1180 : !!
1181 : !! INPUTS
1182 : !! rot_mat = rotation matrix
1183 : !! tndim = dimension of tensor
1184 : !! vee = input tensor
1185 : !!
1186 : !! OUTPUT
1187 : !! vee_rotated = rotated tensor
1188 : !!
1189 : !! SOURCE
1190 :
1191 40 : subroutine rotate_hu(rot_mat,nsppol,tndim,vee,vee_rotated)
1192 :
1193 : !Arguments ------------------------------------
1194 : integer, intent(in) :: nsppol,tndim
1195 : complex(dp), intent(in) :: rot_mat(tndim,tndim,nsppol)
1196 : complex(dp), intent(in) :: vee(tndim*nsppol,tndim*nsppol,tndim*nsppol,tndim*nsppol)
1197 : complex(dp), intent(inout) :: vee_rotated(tndim*nsppol,tndim*nsppol,tndim*nsppol,tndim*nsppol)
1198 : !Local variables-------------------------------
1199 : integer :: is1,is2,loop,m1,m2,ms1,ms2
1200 40 : complex(dp), allocatable :: mat_tmp(:,:),vee_tmp(:,:)
1201 : ! *********************************************************************
1202 :
1203 160 : ABI_MALLOC(mat_tmp,(tndim,tndim))
1204 120 : ABI_MALLOC(vee_tmp,(tndim,tndim))
1205 :
1206 120 : do loop=1,2
1207 200 : do is2=1,nsppol
1208 400 : do m2=1,tndim
1209 240 : ms2 = m2 + (is2-1)*tndim
1210 1040 : do m1=1,tndim
1211 720 : ms1 = m1 + (is2-1)*tndim
1212 1680 : do is1=1,nsppol
1213 :
1214 : ! Make copy here to prevent creation of temporary when calling zgemm
1215 720 : if (loop == 1) then
1216 4680 : vee_tmp(:,:) = vee(1+(is1-1)*tndim:is1*tndim,ms1,1+(is1-1)*tndim:is1*tndim,ms2)
1217 : else
1218 4680 : vee_tmp(:,:) = vee_rotated(ms1,1+(is1-1)*tndim:is1*tndim,ms2,1+(is1-1)*tndim:is1*tndim)
1219 : end if ! loop
1220 :
1221 : call abi_xgemm("c","n",tndim,tndim,tndim,cone,rot_mat(:,:,is1),tndim, &
1222 720 : & vee_tmp(:,:),tndim,czero,mat_tmp(:,:),tndim)
1223 : call abi_xgemm("n","n",tndim,tndim,tndim,cone,mat_tmp(:,:),tndim, &
1224 720 : & rot_mat(:,:,is1),tndim,czero,vee_tmp(:,:),tndim)
1225 :
1226 1440 : if (loop == 1) then
1227 4680 : vee_rotated(1+(is1-1)*tndim:is1*tndim,ms1,1+(is1-1)*tndim:is1*tndim,ms2) = vee_tmp(:,:)
1228 : else
1229 4680 : vee_rotated(ms1,1+(is1-1)*tndim:is1*tndim,ms2,1+(is1-1)*tndim:is1*tndim) = vee_tmp(:,:)
1230 : end if ! loop
1231 :
1232 : end do ! is1
1233 : end do ! m1
1234 : end do ! ms2
1235 : end do ! is2
1236 : end do ! loop
1237 :
1238 40 : ABI_FREE(mat_tmp)
1239 40 : ABI_FREE(vee_tmp)
1240 :
1241 40 : end subroutine rotate_hu
1242 : !!***
1243 :
1244 : !!****f* m_hu/printvee_hu
1245 : !! NAME
1246 : !! printvee_hu
1247 : !!
1248 : !! FUNCTION
1249 : !! Print vee
1250 : !!
1251 : !! INPUTS
1252 : !! vee = tensor for Coulomb interactions
1253 : !!
1254 : !! OUTPUT
1255 : !!
1256 : !! SOURCE
1257 :
1258 112 : subroutine printvee_hu(ndim,vee,prtopt,basis,upawu,f2)
1259 :
1260 : !Arguments ------------------------------------
1261 : !type
1262 : integer, intent(in) :: ndim,prtopt
1263 : complex(dp), intent(in) :: vee(ndim,ndim,ndim,ndim)
1264 : real(dp), optional, intent(in) :: f2,upawu
1265 : character(len=*), intent(in) :: basis
1266 : !Local variables-------------------------------
1267 : integer :: abcomp,m1,m2,mi,mj,mk,ml
1268 112 : real(dp), allocatable :: a2pp(:,:),b0(:,:),b2pp(:,:)
1269 : character(len=2000) :: message
1270 : ! *********************************************************************
1271 :
1272 112 : write(message,'(5a)') ch10,&
1273 224 : & ' Coulomb interaction in the ',trim(basis),' basis'
1274 112 : call wrtout(std_out,message,'COLL')
1275 :
1276 112 : if (prtopt .eq. 2) then
1277 :
1278 108 : write(message,'(2a)') ch10," <mi,mi|vee|mi mi> : U1"
1279 108 : call wrtout(std_out,message,'COLL')
1280 762 : do mi=1,ndim
1281 654 : write(message,'(4i4,3x,e10.3)') mi,mi,mi,mi,dble(vee(mi,mi,mi,mi))
1282 762 : call wrtout(std_out,message,'COLL')
1283 : end do ! mi
1284 :
1285 108 : write(message,'(2a)') ch10," <mi,mj|vee|mi mj> : U2"
1286 108 : call wrtout(std_out,message,'COLL')
1287 762 : do mi=1,ndim
1288 3050 : do mj=mi+1,ndim
1289 2288 : write(message,'(4i4,3x,e10.3)') mi,mj,mi,mj,dble(vee(mi,mj,mi,mj))
1290 2942 : call wrtout(std_out,message,'COLL')
1291 : end do ! mj
1292 : end do ! mi
1293 :
1294 108 : write(message,'(2a)') ch10," <mi,mj|vee|mj mi> : J"
1295 108 : call wrtout(std_out,message,'COLL')
1296 762 : do mi=1,ndim
1297 3050 : do mj=mi+1,ndim
1298 2288 : write(message,'(4i4,3x,e10.3)') mi,mj,mj,mi,dble(vee(mi,mj,mj,mi))
1299 2942 : call wrtout(std_out,message,'COLL')
1300 : end do ! mj
1301 : end do ! mi
1302 :
1303 108 : write(message,'(2a)') ch10," <mi,mi|vee|mj mj> : J"
1304 108 : call wrtout(std_out,message,'COLL')
1305 762 : do mi=1,ndim
1306 3050 : do mj=mi+1,ndim
1307 2288 : write(message,'(4i4,3x,e10.3)') mi,mi,mj,mj,dble(vee(mi,mi,mj,mj))
1308 2942 : call wrtout(std_out,message,'COLL')
1309 : end do ! mj
1310 : end do ! mi
1311 :
1312 108 : write(message,'(2a)') ch10," vee is non zero also for"
1313 108 : call wrtout(std_out,message,'COLL')
1314 762 : do mi=1,ndim
1315 5992 : do mj=1,ndim
1316 58822 : do mk=1,ndim
1317 684686 : do ml=1,ndim
1318 626518 : if ((.not. (mi == mk .and. mj == ml)) .and. (.not. (mi == ml .and. mj == mk)) .and.&
1319 52938 : & .not. (mi == mj .and. mk == ml)) then
1320 612136 : if (dble(vee(mi,mj,mk,ml)) > tol8) then
1321 960 : write(message,'(4i4,3x,e10.3)') mi,mj,mk,ml,dble(vee(mi,mj,mk,ml))
1322 960 : call wrtout(std_out,message,'COLL')
1323 : end if
1324 : end if
1325 : end do ! ml
1326 : end do ! mk
1327 : end do ! mj
1328 : end do ! mi
1329 108 : write(message,'(a)') ch10
1330 108 : call wrtout(std_out,message,'COLL')
1331 :
1332 : end if ! prtopt>=2
1333 :
1334 112 : if (prtopt .eq. 1) then
1335 :
1336 4 : write(message,'(2x,a,3x,14f10.4)') "Um1m2=Vee(m1,m2,m1,m2)"
1337 4 : call wrtout(std_out,message,'COLL')
1338 60 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1339 4 : call wrtout(std_out,message,'COLL')
1340 60 : do m1=1,ndim
1341 840 : write(message,'(2x,i4,3x,14f10.4)') m1,(dble(vee(m1,m2,m1,m2)),m2=1,ndim)
1342 60 : call wrtout(std_out,message,'COLL')
1343 : end do ! m1
1344 4 : write(message,'(a)') ch10
1345 4 : call wrtout(std_out,message,'COLL')
1346 :
1347 4 : write(message,'(2x,a,3x,14f10.4)') "Jm1m2=Vee(m1,m2,m2,m1)"
1348 4 : call wrtout(std_out,message,'COLL')
1349 60 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1350 4 : call wrtout(std_out,message,'COLL')
1351 60 : do m1=1,ndim
1352 840 : write(message,'(2x,i4,3x,14f10.4)') m1,(dble(vee(m1,m2,m2,m1)),m2=1,ndim)
1353 60 : call wrtout(std_out,message,'COLL')
1354 : end do ! m1
1355 4 : write(message,'(a)') ch10
1356 4 : call wrtout(std_out,message,'COLL')
1357 :
1358 4 : write(message,'(2x,a,3x,14f10.4)') "Udens(m1,m2)"
1359 4 : call wrtout(std_out,message,'COLL')
1360 60 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1361 4 : call wrtout(std_out,message,'COLL')
1362 60 : do m1=1,ndim
1363 840 : write(message,'(2x,i4,3x,14f10.4)') m1,(dble(vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)),m2=1,ndim)
1364 60 : call wrtout(std_out,message,'COLL')
1365 : end do ! m1
1366 4 : write(message,'(a)') ch10
1367 4 : call wrtout(std_out,message,'COLL')
1368 :
1369 : !if (prtopt >= 3) then
1370 : ! if(ndim==7) then
1371 : ! averu=zero
1372 : ! do m1=1,7
1373 : ! do m2=1,7
1374 : ! averu=vee(m1,m2,m1,m2)+averu
1375 : ! enddo
1376 : ! enddo
1377 : ! averu=averu/49.d0
1378 : ! averurestricted=zero
1379 : ! do m1=1,7
1380 : ! do m2=1,7
1381 : ! if(m1/=m2) averurestricted=vee(m1,m2,m1,m2)+averurestricted
1382 : ! enddo
1383 : ! enddo
1384 : ! averurestricted=averurestricted/42.d0
1385 : ! averj=zero
1386 : ! do m1=1,7
1387 : ! do m2=1,7
1388 : ! if(m1/=m2) averj=vee(m1,m2,m2,m1)+averj
1389 : ! enddo
1390 : ! enddo
1391 : ! averj=averj/42
1392 : ! averall=zero
1393 : ! do m1=1,7
1394 : ! do m2=1,7
1395 : ! if(m1/=m2) averall=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+averall
1396 : ! enddo
1397 : ! enddo
1398 : ! averall=averall/42
1399 : !write(6,*) "averages Ylm U, U restricted,J, U-J",averu,averurestricted,averj,averall
1400 :
1401 : ! endif
1402 :
1403 : ! if(ndim==14.and.trim(basis)=='jmj') then
1404 : ! aver52=zero
1405 : ! do m1=1,6
1406 : ! do m2=1,6
1407 : ! aver52=vee(m1,m2,m1,m2)+aver52
1408 : ! enddo
1409 : ! enddo
1410 : ! aver52=aver52/36.d0
1411 : ! aver72=zero
1412 : ! do m1=7,14
1413 : ! do m2=7,14
1414 : ! aver72=vee(m1,m2,m1,m2)+aver72
1415 : ! enddo
1416 : ! enddo
1417 : ! aver72=aver72/64.d0
1418 : ! avernondiag=zero
1419 : ! do m1=1,6
1420 : ! do m2=7,14
1421 : ! avernondiag=vee(m1,m2,m1,m2)+avernondiag
1422 : ! enddo
1423 : ! enddo
1424 : ! avernondiag=avernondiag/48.d0
1425 : ! averall=zero
1426 : ! do m1=1,14
1427 : ! do m2=1,14
1428 : ! averall=vee(m1,m2,m1,m2)+averall
1429 : ! enddo
1430 : ! enddo
1431 : ! averall=averall/196.d0
1432 : !write(6,*) "U averages",aver52,aver72,avernondiag,averall
1433 :
1434 :
1435 :
1436 : ! aver52=zero
1437 : ! do m1=1,6
1438 : ! do m2=1,6
1439 : ! if(m1/=m2) aver52=vee(m1,m2,m2,m1)+aver52
1440 : ! enddo
1441 : ! enddo
1442 : ! aver52=aver52/30.d0
1443 : ! aver72=zero
1444 : ! do m1=7,14
1445 : ! do m2=7,14
1446 : ! if(m1/=m2) aver72=vee(m1,m2,m2,m1)+aver72
1447 : ! enddo
1448 : ! enddo
1449 : ! aver72=aver72/56.d0
1450 : ! avernondiag=zero
1451 : ! do m1=1,6
1452 : ! do m2=7,14
1453 : ! avernondiag=vee(m1,m2,m2,m1)+avernondiag
1454 : ! enddo
1455 : ! enddo
1456 : ! avernondiag=avernondiag/48.d0
1457 : ! averall=zero
1458 : ! do m1=1,14
1459 : ! do m2=1,14
1460 : ! if(m1/=m2) averall=vee(m1,m2,m2,m1)+averall
1461 : ! enddo
1462 : ! enddo
1463 : ! averall=averall/182.d0
1464 : !write(6,*) "J averages",aver52,aver72,avernondiag,averall
1465 :
1466 :
1467 :
1468 :
1469 : ! aver52=zero
1470 : ! do m1=1,6
1471 : ! do m2=1,6
1472 : ! if(m1/=m2) aver52=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+aver52
1473 : ! enddo
1474 : ! enddo
1475 : ! aver52=aver52/30.d0
1476 : ! aver72=zero
1477 : ! do m1=7,14
1478 : ! do m2=7,14
1479 : ! if(m1/=m2) aver72=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+aver72
1480 : ! enddo
1481 : ! enddo
1482 : ! aver72=aver72/56.d0
1483 : ! avernondiag=zero
1484 : ! do m1=1,6
1485 : ! do m2=7,14
1486 : ! avernondiag=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+avernondiag
1487 : ! enddo
1488 : ! enddo
1489 : ! avernondiag=avernondiag/48.d0
1490 : ! averall=zero
1491 : ! do m1=1,14
1492 : ! do m2=1,14
1493 : ! if(m1/=m2) averall=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+averall
1494 : ! enddo
1495 : ! enddo
1496 : ! averall=averall/182.d0
1497 : !write(6,*) "U-J averages",aver52,aver72,avernondiag,averall
1498 :
1499 : ! endif
1500 : !endif
1501 :
1502 4 : if (present(upawu)) then
1503 :
1504 8 : ABI_MALLOC(a2pp,(ndim,ndim))
1505 6 : ABI_MALLOC(b2pp,(ndim,ndim))
1506 6 : ABI_MALLOC(b0,(ndim,ndim))
1507 :
1508 : ! write(message,'(2x,a,3x,14f10.4)') "For check with respect to Slater's paper"
1509 : ! call wrtout(std_out, message,'COLL')
1510 : ! ABI_MALLOC(f0,(ndim,ndim,ndim,ndim))
1511 : ! write(message,'(2x,a,3x,14f10.4)') "Vee(m1,m2,m1,m2)-F0*ao(m1,m2)"
1512 : ! call wrtout(std_out, message,'COLL')
1513 : ! write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1514 : ! call wrtout(std_out, message,'COLL')
1515 : !
1516 : ! do m1=1,ndim
1517 : ! write(message,'(2x,i4,3x,14f10.4)') m1,(vee(m1,m2,m1,m2)-upawu,m2=1,ndim)
1518 : ! call wrtout(std_out, message,'COLL')
1519 : ! enddo
1520 : !
1521 : ! f0=zero
1522 : ! do m1=1,ndim
1523 : ! f0(m1,m1,m1,m1)=upawu
1524 : ! enddo
1525 : ! write(message,'(a)') ch10
1526 : ! call wrtout(std_out, message,'COLL')
1527 : ! write(message,'(2x,a,3x,14f10.4)') "Vee(m1,m2,m2,m1)-F0*b0(m1,m2)"
1528 : ! call wrtout(std_out, message,'COLL')
1529 : ! write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1530 : ! call wrtout(std_out, message,'COLL')
1531 : ! do m1=1,ndim
1532 : ! write(message,'(2x,i4,3x,14f10.4)') m1,(vee(m1,m2,m2,m1)-f0(m1,m2,m2,m1),m2=1,ndim)
1533 : ! call wrtout(std_out, message,'COLL')
1534 : ! enddo
1535 : ! ABI_FREE(f0)
1536 :
1537 :
1538 422 : b0(:,:) = zero
1539 30 : do m1=1,ndim
1540 30 : b0(m1,m1) = upawu
1541 : end do
1542 2 : abcomp = 0
1543 2 : if (ndim == 3 .and. present(f2) .and. (trim(basis)=='slm')) then
1544 0 : a2pp(:,:) = RESHAPE((/1,-2,1,-2,4,-2,-1,-2,-1/),(/3,3/))
1545 0 : a2pp(:,:) = a2pp(:,:)/dble(25)*f2 + upawu
1546 0 : b2pp(:,:) = RESHAPE((/1,3,6,3,4,3,6,3,1/),(/3,3/))
1547 0 : b2pp(:,:) = b2pp(:,:)/dble(25)*f2 + b0(:,:)
1548 : abcomp = 1
1549 2 : else if (ndim == 6 .and. present(f2) .and. (trim(basis)=='jmj')) then
1550 0 : ABI_MALLOC(a2pp,(6,6))
1551 0 : ABI_MALLOC(b2pp,(6,6))
1552 : a2pp(:,:) = RESHAPE((/0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1,-1,1,0,&
1553 0 : & 0,-1,1,1,-1,0,0,-1,1,1,-1,0,0,1,-1,-1,1/),(/6,6/))
1554 0 : a2pp(:,:) = a2pp(:,:)/dble(25)*f2 + upawu
1555 : b2pp(:,:) = RESHAPE((/0,0,1,2,3,4,0,0,4,3,2,1,1,4,1,2,2,0,2,3,&
1556 0 : & 2,1,0,2,3,2,2,0,1,2,4,1,0,2,2,1/),(/6,6/))
1557 0 : b2pp(:,:) = b2pp(:,:)/dble(25)*f2 + b0(:,:)
1558 : abcomp = 1
1559 : end if
1560 2 : if (mod(ndim,3) == 0 .and. present(f2) .and. abcomp == 1) then
1561 0 : write(message,'(2x,a)') "Exact result for Umm is"
1562 0 : call wrtout(std_out,message,'COLL')
1563 0 : do m1=1,ndim
1564 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(a2pp(m1,m2),m2=1,ndim)
1565 0 : call wrtout(std_out,message,'COLL')
1566 : end do ! m1
1567 0 : write(message,'(a)') ch10
1568 0 : call wrtout(std_out,message,'COLL')
1569 0 : write(message,'(2x,a,3x,14f10.4)') "Exact result for Jmm is"
1570 0 : call wrtout(std_out,message,'COLL')
1571 0 : do m1=1,ndim
1572 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(b2pp(m1,m2),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 if
1578 2 : ABI_FREE(a2pp)
1579 2 : ABI_FREE(b2pp)
1580 2 : ABI_FREE(b0)
1581 : end if ! present(upawu)
1582 :
1583 : end if ! prtopt>=1
1584 :
1585 : if (prtopt .eq. 4 ) then
1586 : ! write(message,'(2x,a,3x,14f10.4)') "Real part of Um1m2=Vee(m1,m2,m1,m2)"
1587 : ! call wrtout(std_out,message,'COLL')
1588 : ! write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1589 : ! call wrtout(std_out,message,'COLL')
1590 : ! do m1=1,ndim
1591 : ! write(message,'(2x,i4,3x,14f10.4)') m1,(dble(vee(m1,m2,m1,m2)),m2=1,ndim)
1592 : ! call wrtout(std_out,message,'COLL')
1593 : ! end do ! m1
1594 : ! write(message,'(a)') ch10
1595 : ! call wrtout(std_out,message,'COLL')
1596 :
1597 : !write(message,'(2x,a,3x,14f10.4)') "Imag part of Um1m2=Vee(m1,m2,m1,m2)"
1598 : !call wrtout(std_out,message,'COLL')
1599 : !write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1600 : !call wrtout(std_out,message,'COLL')
1601 : !do m1=1,ndim
1602 : ! write(message,'(2x,i4,3x,14f10.6)') m1,(aimag(vee(m1,m2,m1,m2)),m2=1,ndim)
1603 : ! call wrtout(std_out,message,'COLL')
1604 : !end do ! m1
1605 : !write(message,'(a)') ch10
1606 : !call wrtout(std_out,message,'COLL')
1607 :
1608 : ! write(message,'(2x,a,3x,14f10.4)') "Real part of Jm1m2=Vee(m1,m2,m2,m1)"
1609 : ! call wrtout(std_out,message,'COLL')
1610 : ! write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1611 : ! call wrtout(std_out,message,'COLL')
1612 : ! do m1=1,ndim
1613 : ! write(message,'(2x,i4,3x,14f10.4)') m1,(dble(vee(m1,m2,m2,m1)),m2=1,ndim)
1614 : ! call wrtout(std_out,message,'COLL')
1615 : ! end do ! m1
1616 : ! write(message,'(a)') ch10
1617 : ! call wrtout(std_out,message,'COLL')
1618 :
1619 : !write(message,'(2x,a,3x,14f10.4)') "Imag part of Jm1m2=Vee(m1,m2,m2,m1)"
1620 : !call wrtout(std_out,message,'COLL')
1621 : !write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1622 : !call wrtout(std_out,message,'COLL')
1623 : !do m1=1,ndim
1624 : ! write(message,'(2x,i4,3x,14f10.6)') m1,(aimag(vee(m1,m2,m2,m1)),m2=1,ndim)
1625 : ! call wrtout(std_out,message,'COLL')
1626 : !end do ! m1
1627 : !write(message,'(a)') ch10
1628 : !call wrtout(std_out,message,'COLL')
1629 :
1630 : endif !prtopt == 4
1631 112 : end subroutine printvee_hu
1632 : !!***
1633 :
1634 : !!****f* m_hu/vee2udensatom_hu
1635 : !! NAME
1636 : !! vee2udensatom_hu
1637 : !!
1638 : !! FUNCTION
1639 : !! Compute and print density density interaction from full tensor (used for DFT+DMFT)
1640 : !!
1641 : !! INPUTS
1642 : !! ndim = number of orbitals (without counting the spin)
1643 : !! udens_atoms = density-density interactions
1644 : !! veetemp = full interaction tensor
1645 : !! basis = basis of the interaction tensor
1646 : !! prtonly = 0 (default) : compute and print udens_atoms
1647 : !! = 1 : only print udens_atoms
1648 : !! = 2 : print also complex part
1649 : !! OUTPUT
1650 : !!
1651 : !! SOURCE
1652 :
1653 109 : subroutine vee2udensatom_hu(ndim,udens_atoms,veetemp,basis,prtonly)
1654 :
1655 : !Arguments ------------------------------------
1656 : integer, intent(in) :: ndim
1657 : complex(dp), intent(inout) :: udens_atoms(2*ndim,2*ndim)
1658 : !real(dp), intent(in) :: veetemp(nspinor*ndim,nspinor*ndim,nspinor*ndim,nspinor*ndim)
1659 : complex(dp), intent(in) :: veetemp(ndim,ndim,ndim,ndim)
1660 : character(len=*), intent(in) :: basis
1661 : integer, intent(in), optional :: prtonly
1662 : !Local variables-------------------------------
1663 : integer :: m,m1,ms,ms1,prt_only,tndim
1664 : character(len=1000) :: message
1665 : ! *********************************************************************
1666 :
1667 109 : tndim = 2 * ndim
1668 109 : prt_only = 0
1669 109 : if (present(prtonly)) prt_only = prtonly
1670 109 : if (prt_only == 0) then
1671 1720 : udens_atoms(:,:) = czero
1672 : !ij = 0
1673 240 : do ms=1,tndim-1
1674 200 : m = mod(ms-1,ndim) + 1
1675 840 : do ms1=ms+1,tndim
1676 600 : m1 = mod(ms1-1,ndim) + 1
1677 : !ij = ij + 1
1678 600 : if (ms <= ndim .and. ms1 > ndim) then
1679 : ! hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)
1680 360 : udens_atoms(ms,ms1) = veetemp(m,m1,m,m1)
1681 : ! write(6,*)"A", ms,ms1,udens_atoms(ms,ms1)
1682 : else
1683 : ! hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1684 240 : udens_atoms(ms,ms1) = veetemp(m,m1,m,m1) - veetemp(m,m1,m1,m)
1685 : ! write(6,*)"B", ms,ms1,udens_atoms(ms,ms1)
1686 : end if
1687 800 : udens_atoms(ms1,ms) = udens_atoms(ms,ms1)
1688 : end do ! ms1
1689 : end do ! ms
1690 :
1691 : ! else if(nspinor==2) then
1692 : !
1693 : ! do ms=1,2*ndim
1694 : ! do ms1=1,2*ndim
1695 : ! udens_atoms(ms,ms1)=veetemp(ms,ms1,ms,ms1)-veetemp(ms,ms1,ms1,ms)
1696 : ! enddo
1697 : ! enddo
1698 :
1699 : end if ! prt_only=0
1700 :
1701 109 : message=''
1702 109 : write(message,'(4a)') ch10," -------- Interactions in the ",trim(basis)," basis "
1703 109 : call wrtout(std_out,message,'COLL')
1704 109 : if (prtonly .eq. 2) then
1705 0 : write(message,'(14(i14))') (m,m=1,tndim)
1706 0 : call wrtout(std_out,message,'COLL')
1707 0 : do ms=1,tndim
1708 0 : write(message,'(i3,14(2x,28f7.3))') ms,((udens_atoms(ms,ms1)),ms1=1,tndim)
1709 0 : call wrtout(std_out,message,'COLL')
1710 : enddo
1711 : else
1712 1119 : write(message,'(1x,14(2x,i5))') (m,m=1,tndim)
1713 109 : call wrtout(std_out,message,'COLL')
1714 1119 : do ms=1,tndim
1715 11270 : write(message,'(i3,14f7.3)') ms,(dble(udens_atoms(ms,ms1)),ms1=1,tndim)
1716 1119 : call wrtout(std_out,message,'COLL')
1717 : enddo !ms
1718 : endif
1719 109 : write(message,'(3a)') "--------------------------------------------------------"
1720 109 : call wrtout(std_out,message,'COLL')
1721 :
1722 109 : end subroutine vee2udensatom_hu
1723 : !!***
1724 :
1725 : !!****f* m_hu/reddd
1726 : !! NAME
1727 : !! reddd
1728 : !!
1729 : !! FUNCTION
1730 : !!
1731 : !! INPUTS
1732 : !!
1733 : !! OUTPUT
1734 : !!
1735 : !! SOURCE
1736 :
1737 : !function reddd(mi,ndim)
1738 :
1739 : ! use defs_basis
1740 :
1741 : !Arguments ------------------------------------
1742 : !scalars
1743 : ! integer,intent(in) :: mi,ndim
1744 : ! integer :: reddd
1745 : ! *************************************************************************
1746 :
1747 : ! if(mi<ndim+1) reddd=mi
1748 : ! if(mi>=ndim+1) reddd=mi-ndim
1749 :
1750 : !end function reddd
1751 : !!***
1752 :
1753 : !!****f* m_hu/vee_slm2ylm_hu
1754 : !! NAME
1755 : !! vee_slm2ylm_hu
1756 : !!
1757 : !! FUNCTION
1758 : !! For a given angular momentum lcor, change a matrix of interaction of dimension (2*lcor+1)
1759 : !! from the Slm to the Ylm basis if option==1 or from Ylm to Slm if !option==2
1760 : !!
1761 : !! COPYRIGHT
1762 : !! Copyright (C) 1998-2026 ABINIT group (BA)
1763 : !! This file is distributed under the terms of the
1764 : !! GNU General Public License, see ~abinit/COPYING
1765 : !! or http://www.gnu.org/copyleft/gpl.txt .
1766 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
1767 : !!
1768 : !! INPUTS
1769 : !! lcor= angular momentum, size of the matrix is 2*lcor+1
1770 : !! mat_inp_c= Input matrix
1771 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1772 : !! option= 1 Change matrix from Slm to Ylm basis
1773 : !! 2 Change matrix from Ylm to Slm basis
1774 : !! prtvol=printing volume
1775 : !!
1776 : !! OUTPUT
1777 : !! mat_out_c= Output matrix in Ylm or Slm basis according to option
1778 : !!
1779 : !! NOTES
1780 : !!
1781 : !! SOURCE
1782 :
1783 0 : subroutine vee_slm2ylm_hu(lcor,mat_inp_c,mat_out_c,paw_dmft,option,prtvol)
1784 :
1785 : !Arguments ---------------------------------------------
1786 : integer,intent(in) :: lcor,option,prtvol
1787 : complex(dp), intent(in) :: mat_inp_c(2*lcor+1,2*lcor+1,2*lcor+1,2*lcor+1)
1788 : complex(dp), intent(out) :: mat_out_c(2*lcor+1,2*lcor+1,2*lcor+1,2*lcor+1)
1789 : type(paw_dmft_type) , target, intent(in) :: paw_dmft
1790 : !Local variables ---------------------------------------
1791 : integer :: ndim
1792 0 : complex(dp), allocatable :: slm2ylm(:,:)
1793 : character(len=500) :: message
1794 : ! *********************************************************************
1795 :
1796 0 : if (option /= 1 .and. option /= 2) then
1797 0 : message = ' option=/1 or 2 !'
1798 0 : ABI_BUG(message)
1799 : end if
1800 :
1801 0 : if (abs(prtvol) > 2) then
1802 0 : write(message,'(3a)') ch10," vee_slm2ylm_hu"
1803 0 : call wrtout(std_out,message,'COLL')
1804 : end if
1805 :
1806 0 : if (abs(prtvol) > 2) then
1807 0 : if (option == 1) then
1808 0 : write(message,'(3a)') ch10,"matrix in cubic basis is changed into Ylm basis"
1809 0 : else if (option == 2) then
1810 0 : write(message,'(3a)') ch10,"matrix in Ylm basis is changed into cubic basis"
1811 : end if
1812 0 : call wrtout(std_out,message,'COLL')
1813 : end if ! prtvol>2
1814 :
1815 0 : ndim = 2*lcor + 1
1816 :
1817 0 : ABI_MALLOC(slm2ylm,(ndim,ndim))
1818 :
1819 0 : if (option == 1) then
1820 0 : slm2ylm(:,:) = conjg(transpose(paw_dmft%slm2ylm(1:ndim,1:ndim,lcor+1)))
1821 0 : else if (option == 2) then
1822 : ! Make copy here to prevent creation of temporary in case ndim /= ndim_max
1823 0 : slm2ylm(:,:) = paw_dmft%slm2ylm(1:ndim,1:ndim,lcor+1)
1824 : end if
1825 0 : call rotate_hu(slm2ylm(:,:),1,ndim,mat_inp_c(:,:,:,:),mat_out_c(:,:,:,:))
1826 :
1827 0 : ABI_FREE(slm2ylm)
1828 :
1829 : !ll=lcor
1830 : !ABI_MALLOC(slm2ylm,(2*ll+1,2*ll+1))
1831 : !slm2ylm=czero
1832 : !mat_out_c=czero
1833 :
1834 : ! ===== Definitions of slm2ylm
1835 : !do im=1,2*ll+1
1836 : ! mm=im-ll-1;jm=-mm+ll+1 ! mmj=-mm
1837 : !! im is in {1,....2*ll+1}
1838 : !! mm is in {-ll,....+ll}
1839 : !! jm is in {2*ll+1,....,1}
1840 : ! onem=dble((-1)**mm)
1841 : ! if (mm> 0) then ! im in {ll+1,2ll+1} and jm in {ll+1,1}
1842 : ! slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
1843 : ! slm2ylm(jm,im)= cmplx(invsqrt2, zero,kind=dp)
1844 : ! end if
1845 : ! if (mm==0) then
1846 : ! slm2ylm(im,im)=cone
1847 : ! end if
1848 : ! if (mm< 0) then
1849 : ! slm2ylm(im,im)= cmplx(zero, invsqrt2,kind=dp)
1850 : ! slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
1851 : ! end if
1852 : ! end do
1853 : ! do im=1,2*ll+1
1854 : ! write(message,'(7(2f14.5))') (slm2ylm(im,jm),jm=1,2*ll+1)
1855 : ! call wrtout(std_out,message,'COLL')
1856 : ! end do
1857 :
1858 : ! ===== Definitions of slm2ylm
1859 : !!!! pawtab(itypat)%vee(m11,m31,m21,m41)= <m11 m31| vee| m21 m41 >
1860 : !!!! pawtab(itypat)%vee(m11,m21,m31,m41)= <m11 m21| vee| m31 m41 >
1861 :
1862 : !do jm=1,2*ll+1
1863 : ! do im=1,2*ll+1
1864 : ! do hm=1,2*ll+1
1865 : ! do gm=1,2*ll+1
1866 : ! tmp2=czero
1867 : ! do gg=1,2*ll+1
1868 : ! do hh=1,2*ll+1
1869 : ! do ii=1,2*ll+1
1870 : ! do jj=1,2*ll+1
1871 : ! if(option==1) then
1872 : ! tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*(slm2ylm(im,ii))*CONJG(slm2ylm(jm,jj))&
1873 : !& *(slm2ylm(gm,gg))*CONJG(slm2ylm(hm,hh))
1874 : ! if(gm==1.and.hm==1.and.im==1.and.jm==1) then
1875 : ! write(6,'(4i4,2f10.5,2f10.5)') gg,hh,ii,jj,tmp2,mat_inp_c(gg,hh,ii,jj)
1876 : ! write(6,*) "i1"
1877 : ! endif
1878 : ! else if(option==2) then
1879 : ! tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*CONJG(slm2ylm(ii,im))*(slm2ylm(jj,jm))&
1880 : !& *CONJG(slm2ylm(gg,gm))*(slm2ylm(hh,hm))
1881 : ! end if
1882 : ! end do
1883 : ! end do
1884 : ! end do
1885 : ! end do
1886 : ! mat_out_c(gm,hm,im,jm)=tmp2
1887 : ! mat_out_c(gm,im,hm,jm)=tmp2
1888 : ! end do
1889 : ! end do
1890 : ! end do
1891 : !end do
1892 :
1893 : !ABI_FREE(slm2ylm)
1894 :
1895 0 : end subroutine vee_slm2ylm_hu
1896 : !!***
1897 :
1898 : !!****f* m_hu/vee_ndim2tndim_hu_r
1899 : !! NAME
1900 : !! vee_ndim2tndim_hu_r
1901 : !!
1902 : !! FUNCTION
1903 : !! Change a matrix of interaction of dimension [(2*lcor+1)]**2
1904 : !! into a full spin and orbital interaction matrix of dimension [2*(2l+1)]**4
1905 : !!
1906 : !! COPYRIGHT
1907 : !! Copyright (C) 1998-2026 ABINIT group (BA)
1908 : !! This file is distributed under the terms of the
1909 : !! GNU General Public License, see ~abinit/COPYING
1910 : !! or http://www.gnu.org/copyleft/gpl.txt .
1911 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
1912 : !!
1913 : !! INPUTS
1914 : !! lcor= angular momentum, size of the matrix is 2(2*lcor+1)
1915 : !! mat_inp_c= real input matrix
1916 : !! prtvol=printing volume
1917 : !! option= 1 : Vout_(s1m1,s2m2,s3m3,s4m4)=Vinp_(m1,m2,m3,m4)*delta_s1s3*delta_s2s4
1918 : !!
1919 : !!
1920 : !! OUTPUT
1921 : !! mat_out_c= real output matrix
1922 : !!
1923 : !! NOTES
1924 : !!
1925 : !! SOURCE
1926 :
1927 0 : subroutine vee_ndim2tndim_hu_r(lcor,mat_inp_c,mat_out_c,option)
1928 :
1929 : !Arguments ---------------------------------------------
1930 : !scalars
1931 : integer,intent(in) :: lcor,option
1932 : !arrays
1933 : real(dp), intent(in) :: mat_inp_c(:,:,:,:) !real(dp), intent(inout) :: mat_inp_c(:,:,:,:) !(2*lcor+1,2*lcor+1,2*lcor+1,2*lcor+1) real(dp), intent(in), pointer :: mat_inp_c(:,:,:,:)
1934 : real(dp), intent(out) :: mat_out_c(:,:,:,:) !(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
1935 :
1936 : !Local variables ---------------------------------------
1937 : !scalars
1938 : integer :: m1,m2,m3,m4,is1,is2,is3,is4,ndim,s1,s2,s3,s4
1939 :
1940 : ! *********************************************************************
1941 0 : ndim=2*lcor+1
1942 0 : mat_out_c=czero
1943 :
1944 0 : if(option==1) then
1945 0 : do m1=1,ndim
1946 0 : do m2=1,ndim
1947 0 : do m3=1,ndim
1948 0 : do m4=1,ndim
1949 0 : do is1=1,2
1950 0 : do is2=1,2
1951 :
1952 0 : is3=is1 ; is4=is2
1953 :
1954 0 : s1=(is1-1)*ndim ; s2=(is2-1)*ndim ; s3=(is3-1)*ndim ; s4=(is4-1)*ndim
1955 :
1956 0 : mat_out_c(m1+s1,m2+s2,m3+s3,m4+s4)= mat_inp_c(m1,m2,m3,m4)
1957 :
1958 : enddo
1959 : enddo
1960 : enddo
1961 : enddo
1962 : enddo
1963 : enddo
1964 : endif
1965 :
1966 :
1967 0 : end subroutine vee_ndim2tndim_hu_r
1968 : !!***
1969 :
1970 : !!****f* m_hu/vee_ndim2tndim_hu
1971 : !! NAME
1972 : !! vee_ndim2tndim_hu
1973 : !!
1974 : !! FUNCTION
1975 : !! Change a matrix of interaction of dimension [(2*lcor+1)]**4
1976 : !! into a full spin and orbital interaction matrix of dimension [2*(2l+1)]**4
1977 : !!
1978 : !! COPYRIGHT
1979 : !! Copyright (C) 1998-2026 ABINIT group (BA)
1980 : !! This file is distributed under the terms of the
1981 : !! GNU General Public License, see ~abinit/COPYING
1982 : !! or http://www.gnu.org/copyleft/gpl.txt .
1983 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
1984 : !!
1985 : !! INPUTS
1986 : !! lcor= angular momentum, size of the matrix is 2(2*lcor+1)
1987 : !! mat_inp_c= input matrix
1988 : !!
1989 : !! OUTPUT
1990 : !! mat_out_c= Complex output matrix
1991 : !!
1992 : !! NOTES
1993 : !!
1994 : !! SOURCE
1995 :
1996 133 : subroutine vee_ndim2tndim_hu(lcor,mat_inp_c,mat_out_c)
1997 :
1998 : !Arguments ---------------------------------------------
1999 : !scalars
2000 : integer, intent(in) :: lcor
2001 : !arrays
2002 : complex(dp), intent(in) :: mat_inp_c(2*lcor+1,2*lcor+1,2*lcor+1,2*lcor+1)
2003 : complex(dp), intent(inout) :: mat_out_c(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
2004 : !Local variables ---------------------------------------
2005 : !scalars
2006 : integer :: is1,is2,m1,m2,m3,m4,ndim,s1,s2
2007 : ! *********************************************************************
2008 :
2009 133 : ndim = 2*lcor + 1
2010 1160811 : mat_out_c(:,:,:,:) = czero
2011 :
2012 399 : do is2=1,2
2013 266 : s2 = (is2-1) * ndim
2014 1449 : do m4=1,ndim
2015 3416 : do is1=1,2
2016 2100 : s1 = (is1-1) * ndim
2017 12514 : do m3=1,ndim
2018 58684 : do m2=1,ndim
2019 320540 : do m1=1,ndim
2020 311176 : mat_out_c(m1+s1,m2+s2,m3+s1,m4+s2) = mat_inp_c(m1,m2,m3,m4)
2021 : end do ! m1
2022 : end do ! m2
2023 : end do ! m3
2024 : end do ! is1
2025 : end do ! m4
2026 : end do ! is2
2027 :
2028 133 : end subroutine vee_ndim2tndim_hu
2029 : !!***
2030 :
2031 : !!****f* ABINIT/vee_ylm2jmj_hu
2032 : !! NAME
2033 : !! vee_ylm2jmj_hu
2034 : !!
2035 : !! FUNCTION
2036 : !! For a given angular momentum lcor, change a matrix of dimension [2(2*lcor+1)]**4
2037 : !! from the Ylm basis to the J,M_J basis if option==1
2038 : !!
2039 : !! COPYRIGHT
2040 : !! Copyright (C) 1998-2026 ABINIT group (BA)
2041 : !! This file is distributed under the terms of the
2042 : !! GNU General Public License, see ~abinit/COPYING
2043 : !! or http://www.gnu.org/copyleft/gpl.txt .
2044 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
2045 : !!
2046 : !! INPUTS
2047 : !! lcor= angular momentum
2048 : !! mat_inp_c = input tensor
2049 : !! mat_out_c = output tensor
2050 : !! option= 1 matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis
2051 : !! 2 matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis
2052 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
2053 : !!
2054 : !! SIDE EFFECTS
2055 : !!
2056 : !! NOTES
2057 : !! useful only in ndij==4
2058 : !!
2059 : !! SOURCE
2060 :
2061 0 : subroutine vee_ylm2jmj_hu(lcor,mat_inp_c,mat_out_c,option,paw_dmft)
2062 :
2063 : !Arguments ---------------------------------------------
2064 : integer, intent(in) :: lcor,option
2065 : complex(dp), intent(in) :: mat_inp_c(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
2066 : complex(dp), intent(inout) :: mat_out_c(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
2067 : type(paw_dmft_type), intent(in) :: paw_dmft
2068 : !Local variables ---------------------------------------
2069 : integer :: im,jm,tndim
2070 0 : complex(dp), allocatable :: jmj2ylm(:,:)
2071 : character(len=500) :: message
2072 : !*********************************************************************
2073 :
2074 0 : if (option /= 1 .and. option /= 2) then
2075 0 : message = ' option=/1 and =/2 !'
2076 0 : ABI_BUG(message)
2077 : end if
2078 :
2079 0 : if (lcor == 0) ABI_BUG("l should not be equal to 0")
2080 :
2081 0 : if (option == 1) then
2082 0 : write(message,'(3a)') ch10,"matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis"
2083 0 : else if (option == 2) then
2084 0 : write(message,'(3a)') ch10,"matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis"
2085 : end if
2086 0 : call wrtout(std_out,message,"COLL")
2087 :
2088 0 : tndim = 2 * (2*lcor+1)
2089 :
2090 : ! Make copy to prevent creation of temporary in the case ndim /= ndim_max
2091 0 : ABI_MALLOC(jmj2ylm,(tndim,tndim))
2092 :
2093 0 : if (option == 1) then
2094 0 : jmj2ylm(:,:) = paw_dmft%jmj2ylm(1:tndim,1:tndim,lcor+1)
2095 0 : else if (option == 2) then
2096 0 : jmj2ylm(:,:) = conjg(transpose(paw_dmft%jmj2ylm(1:tndim,1:tndim,lcor+1)))
2097 : end if
2098 :
2099 0 : write(message,'(3a)') ch10,"Matrix to go from |J,M_J> to |M_L,M_S>"
2100 0 : call wrtout(std_out,message,"COLL")
2101 0 : do im=1,tndim
2102 0 : write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (jmj2ylm(im,jm),jm=1,tndim)
2103 0 : call wrtout(std_out,message,"COLL")
2104 : end do
2105 :
2106 0 : call rotate_hu(jmj2ylm(:,:),1,tndim,mat_inp_c(:,:,:,:),mat_out_c(:,:,:,:))
2107 :
2108 0 : ABI_FREE(jmj2ylm)
2109 :
2110 : !--------------- Built indices + allocations
2111 : !ll=lcor
2112 : !ABI_MALLOC(mlms2jmj,(2*(2*ll+1),2*(2*ll+1)))
2113 : !mlms2jmj=czero
2114 : !ABI_MALLOC(ind_msml,(2,-ll:ll))
2115 : !mlms2jmj=czero
2116 : !jc1=0
2117 : !do ms1=1,2
2118 : ! do ml1=-ll,ll
2119 : ! jc1=jc1+1
2120 : ! ind_msml(ms1,ml1)=jc1
2121 : ! end do
2122 : !end do
2123 :
2124 : !--------------- built mlms2jmj
2125 : !do jj=ll,ll+1 ! the physical value of j are ll-0.5,ll+0.5
2126 : !xj(jj)=jj-0.5
2127 : !if(ll==0)then
2128 : ! message=' ll should not be equal to zero !'
2129 : ! ABI_BUG(message)
2130 : !end if
2131 : !jc1=0
2132 : !invsqrt2lp1=one/sqrt(float(2*lcor+1))
2133 : !do jj=ll,ll+1
2134 : ! xj=float(jj)-half ! xj is in {ll-0.5, ll+0.5}
2135 : ! do jm=-jj,jj-1
2136 : ! xmj=float(jm)+half ! xmj is in {-xj,xj}
2137 : ! jc1=jc1+1 ! Global index for JMJ
2138 : ! if(nint(xj+0.5)==ll+1) then ! if xj=ll+0.5
2139 : ! if(nint(xmj+0.5)==ll+1) then
2140 : ! mlms2jmj(ind_msml(2,ll),jc1)=1.0 ! J=L+0.5 and m_J=L+0.5
2141 : ! else if(nint(xmj-0.5)==-ll-1) then
2142 : ! mlms2jmj(ind_msml(1,-ll),jc1)=1.0 ! J=L+0.5 and m_J=-L-0.5
2143 : ! else
2144 : ! mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
2145 : ! mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
2146 : ! end if
2147 : ! end if
2148 : ! if(nint(xj+0.5)==ll) then ! if xj=ll-0.5
2149 : ! mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
2150 : ! mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
2151 : ! end if
2152 : ! end do
2153 : !end do
2154 : !write(message,'(3a)') ch10,"Matrix to go from |M_L,M_S> to |J,M_J>"
2155 : !call wrtout(std_out,message,"COLL")
2156 : !do im=1,2*(ll*2+1)
2157 : ! write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mlms2jmj(im,jm),jm=1,2*(ll*2+1))
2158 : ! call wrtout(std_out,message,"COLL")
2159 : !end do
2160 :
2161 : !--------------- compute change of basis
2162 : ! do jm=1,2*(2*ll+1)
2163 : ! do im=1,2*(2*ll+1)
2164 : ! do hm=1,2*(2*ll+1)
2165 : ! do gm=1,2*(2*ll+1)
2166 : ! tmp2=czero
2167 : ! do gg=1,2*(2*ll+1)
2168 : ! do hh=1,2*(2*ll+1)
2169 : ! do ii=1,2*(2*ll+1)
2170 : ! do jj=1,2*(2*ll+1)
2171 : ! if(option==1) then
2172 : ! tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*CONJG(mlms2jmj(ii,im))*(mlms2jmj(jj,jm))&
2173 : !& *CONJG(mlms2jmj(gg,gm))*(mlms2jmj(hh,hm))
2174 : ! else if(option==2) then
2175 : ! tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*CONJG(mlms2jmj(ii,im))*(mlms2jmj(jj,jm))& ! inv=t*
2176 : !& *CONJG(mlms2jmj(gg,gm))*(mlms2jmj(hh,hm)) ! inv=t*
2177 : ! tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*(mlms2jmj(im,ii))*CONJG(mlms2jmj(jm,jj))& ! inv=t*
2178 : !& *(mlms2jmj(gm,gg))*CONJG(mlms2jmj(hm,hh)) ! inv=t*
2179 : ! end if
2180 : ! end do
2181 : ! end do
2182 : ! end do
2183 : ! end do
2184 : ! mat_out_c(gm,im,hm,jm)=tmp2
2185 : ! end do
2186 : ! end do
2187 : ! end do
2188 : ! end do
2189 : ! ABI_FREE(mlms2jmj)
2190 : ! ABI_FREE(ind_msml)
2191 :
2192 0 : end subroutine vee_ylm2jmj_hu
2193 : !!***
2194 :
2195 : !!****f* ABINIT/udens_slatercondon_hu
2196 : !! NAME
2197 : !! udens_slatercondon_hu
2198 : !!
2199 : !! FUNCTION
2200 : !! For a given angular momentum l and Slater integrals, give the
2201 : !! density density interactions U(m,m') and J(m,m') from Slater and
2202 : !! Condon tables
2203 : !!
2204 : !! COPYRIGHT
2205 : !! Copyright (C) 1998-2026 ABINIT group (BA)
2206 : !! This file is distributed under the terms of the
2207 : !! GNU General Public License, see ~abinit/COPYING
2208 : !! or http://www.gnu.org/copyleft/gpl.txt .
2209 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
2210 : !!
2211 : !! INPUTS
2212 : !! lcor= angular momentum
2213 : !! fk(lcor+1)= Slater integrals
2214 : !!
2215 : !! SIDE EFFECTS
2216 : !!
2217 : !! NOTES
2218 : !!
2219 : !! SOURCE
2220 :
2221 0 : subroutine udens_slatercondon_hu(fk,lcor)
2222 :
2223 : !Arguments ---------------------------------------------
2224 : !scalars
2225 : integer, intent(in) :: lcor
2226 : real(dp), intent(in) :: fk(0:lcor)
2227 : !Local variables ---------------------------------------
2228 : !scalars
2229 : character(len=500) :: message
2230 : integer :: m1,m2
2231 : !arrays
2232 0 : real(dp), allocatable :: aklmlmp(:,:,:),bklmlmp(:,:,:),jdens(:,:),udens(:,:)
2233 : !*********************************************************************
2234 :
2235 0 : ABI_MALLOC(aklmlmp,(0:lcor,-lcor:lcor,-lcor:lcor)) ! k,m,m'
2236 0 : ABI_MALLOC(bklmlmp,(0:lcor,-lcor:lcor,-lcor:lcor)) ! k,m,m'
2237 0 : ABI_MALLOC(udens,(-lcor:lcor,-lcor:lcor)) ! m,m'
2238 0 : ABI_MALLOC(jdens,(-lcor:lcor,-lcor:lcor)) ! m,m'
2239 : ! k=2*(lcor)
2240 0 : aklmlmp(:,:,:) = zero
2241 0 : bklmlmp(:,:,:) = zero
2242 0 : udens(:,:) = zero
2243 0 : jdens(:,:) = zero
2244 0 : if (lcor == 0) then
2245 0 : aklmlmp(0,0,0)=1
2246 : !
2247 0 : bklmlmp(0,0,0)=0
2248 0 : else if (lcor == 1) then
2249 0 : aklmlmp(0, :, :)=1
2250 0 : aklmlmp(1, 1, 1)= one/25._dp
2251 0 : aklmlmp(1,-1,-1)= one/25._dp
2252 0 : aklmlmp(1,-1, 1)= one/25._dp
2253 0 : aklmlmp(1, 1,-1)= one/25._dp
2254 0 : aklmlmp(1, 1, 0)=-two/25._dp
2255 0 : aklmlmp(1,-1, 0)=-two/25._dp
2256 0 : aklmlmp(1, 0,-1)=-two/25._dp
2257 0 : aklmlmp(1, 0, 1)=-two/25._dp
2258 0 : aklmlmp(1, 0, 0)= four/25._dp
2259 : !
2260 0 : bklmlmp(0, 1, 1)= one
2261 0 : bklmlmp(0,-1,-1)= one
2262 0 : bklmlmp(0, 0, 0)= one
2263 0 : bklmlmp(1, 1, 1)= one/25._dp
2264 0 : bklmlmp(1,-1,-1)= one/25._dp
2265 0 : bklmlmp(1,-1,+1)= six/25._dp
2266 0 : bklmlmp(1,+1,-1)= six/25._dp
2267 0 : bklmlmp(1,+1, 0)= three/25._dp
2268 0 : bklmlmp(1,-1, 0)= three/25._dp
2269 0 : bklmlmp(1, 0,-1)= three/25._dp
2270 0 : bklmlmp(1, 0, 1)= three/25._dp
2271 0 : bklmlmp(1, 0, 0)= four/25._dp
2272 0 : else if (lcor == 2) then
2273 0 : aklmlmp(0, :, :)=1
2274 0 : aklmlmp(1, 2, 2)= four / 49._dp
2275 0 : aklmlmp(1, 2, 1)= -two / 49._dp
2276 0 : aklmlmp(1, 2, 0)= -four / 49._dp
2277 0 : aklmlmp(1, 2,-1)= -two / 49._dp
2278 0 : aklmlmp(1, 2,-2)= four / 49._dp
2279 0 : aklmlmp(1, 1, 2)= -two / 49._dp
2280 0 : aklmlmp(1, 1, 1)= one / 49._dp
2281 0 : aklmlmp(1, 1, 0)= two / 49._dp
2282 0 : aklmlmp(1, 1,-1)= one / 49._dp
2283 0 : aklmlmp(1, 1,-2)= -two / 49._dp
2284 0 : aklmlmp(1, 0, 2)= -four / 49._dp
2285 0 : aklmlmp(1, 0, 1)= two / 49._dp
2286 0 : aklmlmp(1, 0, 0)= four / 49._dp
2287 0 : aklmlmp(1, 0,-1)= two / 49._dp
2288 0 : aklmlmp(1, 0,-2)= -four / 49._dp
2289 0 : aklmlmp(1,-1, 2)= -two / 49._dp
2290 0 : aklmlmp(1,-1, 1)= one / 49._dp
2291 0 : aklmlmp(1,-1, 0)= two / 49._dp
2292 0 : aklmlmp(1,-1,-1)= one / 49._dp
2293 0 : aklmlmp(1,-1,-2)= -two / 49._dp
2294 0 : aklmlmp(1,-2, 2)= four / 49._dp
2295 0 : aklmlmp(1,-2, 1)= -two / 49._dp
2296 0 : aklmlmp(1,-2, 0)= -four / 49._dp
2297 0 : aklmlmp(1,-2,-1)= -two / 49._dp
2298 0 : aklmlmp(1,-2,-2)= four / 49._dp
2299 :
2300 0 : aklmlmp(2, 2, 2)= one / 441._dp
2301 0 : aklmlmp(2, 2, 1)= -four / 441._dp
2302 0 : aklmlmp(2, 2, 0)= six / 441._dp
2303 0 : aklmlmp(2, 2,-1)= -four / 441._dp
2304 0 : aklmlmp(2, 2,-2)= one / 441._dp
2305 0 : aklmlmp(2, 1, 2)= -four / 441._dp
2306 0 : aklmlmp(2, 1, 1)= 16._dp / 441._dp
2307 0 : aklmlmp(2, 1, 0)= -24._dp / 441._dp
2308 0 : aklmlmp(2, 1,-1)= 16._dp / 441._dp
2309 0 : aklmlmp(2, 1,-2)= -four / 441._dp
2310 0 : aklmlmp(2, 0, 2)= six / 441._dp
2311 0 : aklmlmp(2, 0, 1)= -24._dp / 441._dp
2312 0 : aklmlmp(2, 0, 0)= 36._dp / 441._dp
2313 0 : aklmlmp(2, 0,-1)= -24._dp / 441._dp
2314 0 : aklmlmp(2, 0,-2)= six / 441._dp
2315 0 : aklmlmp(2,-1, 2)= -four / 441._dp
2316 0 : aklmlmp(2,-1, 1)= 16._dp / 441._dp
2317 0 : aklmlmp(2,-1, 0)= -24._dp / 441._dp
2318 0 : aklmlmp(2,-1,-1)= 16._dp / 441._dp
2319 0 : aklmlmp(2,-1,-2)= -four / 441._dp
2320 0 : aklmlmp(2,-2, 2)= one / 441._dp
2321 0 : aklmlmp(2,-2, 1)= -four / 441._dp
2322 0 : aklmlmp(2,-2, 0)= six / 441._dp
2323 0 : aklmlmp(2,-2,-1)= -four / 441._dp
2324 0 : aklmlmp(2,-2,-2)= one / 441._dp
2325 : !do m1=lcor,-lcor,-1
2326 : ! do m2=lcor,-lcor,-1
2327 : ! write(6,*) m1,m2,aklmlmp(2,m1,m2)
2328 : ! enddo
2329 : !enddo
2330 :
2331 : !write(message,'(2x,a,3x,14f10.4)') " Slater aklmlmp(2,:,:)"
2332 : !call wrtout(std_out, message,'COLL')
2333 : !write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2334 : !call wrtout(std_out, message,'COLL')
2335 : !do m1=-lcor,lcor,1
2336 : ! write(message,'(2x,i4,3x,14f10.4)') m1,(aklmlmp(2,m1,m2),m2=-lcor,lcor,1)
2337 : ! call wrtout(std_out, message,'COLL')
2338 : !enddo
2339 :
2340 0 : bklmlmp(0, 2, 2)=1
2341 0 : bklmlmp(0, 1, 1)=1
2342 0 : bklmlmp(0, 0, 0)=1
2343 0 : bklmlmp(0,-1,-1)=1
2344 0 : bklmlmp(0,-2,-2)=1
2345 0 : bklmlmp(1, 2, 2)= four / 49._dp
2346 0 : bklmlmp(1, 2, 1)= six / 49._dp
2347 0 : bklmlmp(1, 2, 0)= four / 49._dp
2348 0 : bklmlmp(1, 2,-1)= zero
2349 0 : bklmlmp(1, 2,-2)= zero
2350 0 : bklmlmp(1, 1, 2)= six / 49._dp
2351 0 : bklmlmp(1, 1, 1)= one / 49._dp
2352 0 : bklmlmp(1, 1, 0)= one / 49._dp
2353 0 : bklmlmp(1, 1,-1)= six / 49._dp
2354 0 : bklmlmp(1, 1,-2)= zero
2355 0 : bklmlmp(1, 0, 2)= four / 49._dp
2356 0 : bklmlmp(1, 0, 1)= one / 49._dp
2357 0 : bklmlmp(1, 0, 0)= four / 49._dp
2358 0 : bklmlmp(1, 0,-1)= one / 49._dp
2359 0 : bklmlmp(1, 0,-2)= four / 49._dp
2360 0 : bklmlmp(1,-1, 2)= zero
2361 0 : bklmlmp(1,-1, 1)= six / 49._dp
2362 0 : bklmlmp(1,-1, 0)= one / 49._dp
2363 0 : bklmlmp(1,-1,-1)= one / 49._dp
2364 0 : bklmlmp(1,-1,-2)= six / 49._dp
2365 0 : bklmlmp(1,-2, 2)= zero
2366 0 : bklmlmp(1,-2, 1)= zero
2367 0 : bklmlmp(1,-2, 0)= four / 49._dp
2368 0 : bklmlmp(1,-2,-1)= six / 49._dp
2369 0 : bklmlmp(1,-2,-2)= four / 49._dp
2370 :
2371 0 : bklmlmp(2, 2, 2)= one / 441._dp
2372 0 : bklmlmp(2, 2, 1)= five / 441._dp
2373 0 : bklmlmp(2, 2, 0)= 15._dp / 441._dp
2374 0 : bklmlmp(2, 2,-1)= 35._dp / 441._dp
2375 0 : bklmlmp(2, 2,-2)= 70._dp / 441._dp
2376 0 : bklmlmp(2, 1, 2)= five / 441._dp
2377 0 : bklmlmp(2, 1, 1)= 16._dp / 441._dp
2378 0 : bklmlmp(2, 1, 0)= 30._dp / 441._dp
2379 0 : bklmlmp(2, 1,-1)= 40._dp / 441._dp
2380 0 : bklmlmp(2, 1,-2)= 35._dp / 441._dp
2381 0 : bklmlmp(2, 0, 2)= 15._dp / 441._dp
2382 0 : bklmlmp(2, 0, 1)= 30._dp / 441._dp
2383 0 : bklmlmp(2, 0, 0)= 36._dp / 441._dp
2384 0 : bklmlmp(2, 0,-1)= 30._dp / 441._dp
2385 0 : bklmlmp(2, 0,-2)= 15._dp / 441._dp
2386 0 : bklmlmp(2,-1, 2)= 35._dp / 441._dp
2387 0 : bklmlmp(2,-1, 1)= 40._dp / 441._dp
2388 0 : bklmlmp(2,-1, 0)= 30._dp / 441._dp
2389 0 : bklmlmp(2,-1,-1)= 16._dp / 441._dp
2390 0 : bklmlmp(2,-1,-2)= five / 441._dp
2391 0 : bklmlmp(2,-2, 2)= 70._dp / 441._dp
2392 0 : bklmlmp(2,-2, 1)= 35._dp / 441._dp
2393 0 : bklmlmp(2,-2, 0)= 15._dp / 441._dp
2394 0 : bklmlmp(2,-2,-1)= five / 441._dp
2395 0 : bklmlmp(2,-2,-2)= one / 441._dp
2396 :
2397 : !write(message,'(2x,a,3x,14f10.4)') " Slater bklmlmp(2,:,:)"
2398 : !call wrtout(std_out, message,'COLL')
2399 : !write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2400 : !call wrtout(std_out, message,'COLL')
2401 : !do m1=-lcor,lcor,1
2402 : ! write(message,'(2x,i4,3x,14f10.4)') m1,(bklmlmp(2,m1,m2),m2=-lcor,lcor,1)
2403 : ! call wrtout(std_out, message,'COLL')
2404 : !enddo
2405 : ! if f4of2_sla: these data agree with the explicit calculation
2406 : !else if(lcor==3) then
2407 : end if ! lcor
2408 :
2409 0 : do m2=-lcor,lcor,1
2410 0 : do m1=-lcor,lcor,1
2411 0 : udens(m1,m2) = sum(fk(:)*aklmlmp(:,m1,m2))
2412 0 : jdens(m1,m2) = sum(fk(:)*bklmlmp(:,m1,m2))
2413 : !write(6,*) kk,m1,m2
2414 : !write(6,*) "--",fk(kk),aklmlmp(kk,m1,m2)
2415 : !write(6,*) "--",fk(kk),bklmlmp(kk,m1,m2)
2416 : !udens(m1,m2)=udens(m1,m2)+fk(kk)*aklmlmp(kk,m1,m2)
2417 : !jdens(m1,m2)=jdens(m1,m2)+fk(kk)*bklmlmp(kk,m1,m2)
2418 : end do ! m1
2419 : end do ! m2
2420 0 : write(message,'(2x,a,3x,14f10.4)') " Direct Interaction Matrix from Slater tables (in the Ylm basis) "
2421 0 : call wrtout(std_out,message,'COLL')
2422 0 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2423 0 : call wrtout(std_out,message,'COLL')
2424 0 : do m1=-lcor,lcor,1
2425 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2),m2=-lcor,lcor,1)
2426 0 : call wrtout(std_out,message,'COLL')
2427 : end do ! m1
2428 :
2429 0 : write(message,'(a,2x,a,3x,14f10.4)') ch10," Exchange Interaction Matrix from Slater tables (in the Ylm basis) "
2430 0 : call wrtout(std_out,message,'COLL')
2431 0 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2432 0 : call wrtout(std_out,message,'COLL')
2433 0 : do m1=-lcor,lcor,1
2434 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(jdens(m1,m2),m2=-lcor,lcor,1)
2435 0 : call wrtout(std_out,message,'COLL')
2436 : end do ! m1
2437 :
2438 0 : write(message,'(a,2x,a,3x,14f10.4)') ch10," Density density Interaction Matrix from Slater tables (in the Ylm basis) "
2439 0 : call wrtout(std_out,message,'COLL')
2440 0 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2441 0 : call wrtout(std_out,message,'COLL')
2442 0 : do m1=-lcor,lcor,1
2443 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2)-jdens(m1,m2),m2=-lcor,lcor,1)
2444 0 : call wrtout(std_out,message,'COLL')
2445 : end do ! m1
2446 :
2447 0 : ABI_FREE(jdens)
2448 0 : ABI_FREE(udens)
2449 0 : ABI_FREE(aklmlmp)
2450 0 : ABI_FREE(bklmlmp)
2451 :
2452 0 : end subroutine udens_slatercondon_hu
2453 : !!***
2454 :
2455 : !!****f* ABINIT/udens_inglis_hu
2456 : !! NAME
2457 : !! udens_inglis_hu
2458 : !!
2459 : !! FUNCTION
2460 : !! For a given angular momentum l and Slater integrals, give the
2461 : !! density density interactions U(m,m') and J(m,m') from Inglis tables
2462 : !! in JMJ Basis
2463 : !!
2464 : !! COPYRIGHT
2465 : !! Copyright (C) 1998-2026 ABINIT group (BA)
2466 : !! This file is distributed under the terms of the
2467 : !! GNU General Public License, see ~abinit/COPYING
2468 : !! or http://www.gnu.org/copyleft/gpl.txt .
2469 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
2470 : !!
2471 : !! INPUTS
2472 : !! lcor= angular momentum
2473 : !! fk(lcor+1)= Slater integrals
2474 : !!
2475 : !! SIDE EFFECTS
2476 : !!
2477 : !! NOTES
2478 : !!
2479 : !! SOURCE
2480 :
2481 0 : subroutine udens_inglis_hu(fk,lcor)
2482 :
2483 : !Arguments ---------------------------------------------
2484 : !scalars
2485 : integer, intent(in) :: lcor
2486 : real(dp), intent(in) :: fk(0:lcor)
2487 : !Local variables ---------------------------------------
2488 : !scalars
2489 : character(len=500) :: message
2490 : integer :: m1,m2,tndim
2491 : !arrays
2492 0 : real(dp), allocatable :: a2pp(:,:),app(:,:,:),b2pp(:,:),bpp(:,:,:),jdens(:,:),udens(:,:)
2493 : !*********************************************************************
2494 :
2495 0 : tndim = 2 * (2*lcor+1)
2496 0 : ABI_MALLOC(app,(0:lcor,tndim,tndim))
2497 0 : ABI_MALLOC(bpp,(0:lcor,tndim,tndim))
2498 0 : ABI_MALLOC(a2pp,(tndim,tndim))
2499 0 : ABI_MALLOC(b2pp,(tndim,tndim))
2500 :
2501 0 : ABI_MALLOC(udens,(tndim,tndim))
2502 0 : ABI_MALLOC(jdens,(tndim,tndim))
2503 :
2504 0 : udens(:,:) = zero
2505 0 : jdens(:,:) = zero
2506 0 : a2pp(:,:) = zero
2507 0 : b2pp(:,:) = zero
2508 0 : app(:,:,:) = zero
2509 0 : bpp(:,:,:) = zero
2510 0 : if (lcor == 1) then
2511 0 : app(0,:,:)=one
2512 : a2pp(:,:)=RESHAPE((/0._dp,0._dp, 0._dp, 0._dp, 0._dp, 0._dp,&
2513 : & 0._dp,0._dp, 0._dp, 0._dp, 0._dp, 0._dp,&
2514 : & 0._dp,0._dp, 1._dp,-1._dp,-1._dp, 1._dp,&
2515 : & 0._dp,0._dp,-1._dp, 1._dp, 1._dp,-1._dp,&
2516 : & 0._dp,0._dp,-1._dp, 1._dp, 1._dp,-1._dp,&
2517 0 : & 0._dp,0._dp, 1._dp,-1._dp,-1._dp, 1._dp/),(/6,6/))
2518 0 : app(1,:,:)=a2pp(:,:)
2519 0 : app(1,:,:)=app(1,:,:)/25._dp
2520 :
2521 : b2pp(:,:)=RESHAPE((/1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2522 : & 0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,&
2523 : & 0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,&
2524 : & 0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,&
2525 : & 0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,&
2526 0 : & 0._dp,0._dp,0._dp,0._dp,0._dp,1._dp /),(/6,6/))
2527 0 : bpp(0,:,:)=b2pp(:,:)
2528 : b2pp(:,:)=RESHAPE((/0._dp,0._dp,1._dp,2._dp,3._dp,4._dp,&
2529 : & 0._dp,0._dp,4._dp,3._dp,2._dp,1._dp,&
2530 : & 1._dp,4._dp,1._dp,2._dp,2._dp,0._dp,&
2531 : & 2._dp,3._dp,2._dp,1._dp,0._dp,2._dp,&
2532 : & 3._dp,2._dp,2._dp,0._dp,1._dp,2._dp,&
2533 0 : & 4._dp,1._dp,0._dp,2._dp,2._dp,1._dp /),(/6,6/))
2534 0 : bpp(1,:,:)=b2pp(:,:)
2535 0 : bpp(1,:,:)=bpp(1,:,:)/25._dp
2536 : end if ! lcor=1
2537 0 : if (lcor == 2) then
2538 0 : app(0,:,:)=one
2539 : a2pp(:,:)=RESHAPE((/ 49._dp,-49._dp,-49._dp, 49._dp, 70._dp,-14._dp,-56._dp,-56._dp,-14._dp, 70._dp,&
2540 : & -49._dp, 49._dp, 49._dp,-49._dp,-70._dp, 14._dp, 56._dp, 56._dp, 14._dp,-70._dp,&
2541 : & -49._dp, 49._dp, 49._dp,-49._dp,-70._dp, 14._dp, 56._dp, 56._dp, 14._dp,-70._dp,&
2542 : & 49._dp,-49._dp,-49._dp, 49._dp, 70._dp,-14._dp,-56._dp,-56._dp,-14._dp, 70._dp,&
2543 : & 70._dp,-70._dp,-70._dp, 70._dp,100._dp,-20._dp,-80._dp,-80._dp,-20._dp,100._dp,&
2544 : & -14._dp, 14._dp, 14._dp,-14._dp,-20._dp, 4._dp, 16._dp, 16._dp, 4._dp,-20._dp,&
2545 : & -56._dp, 56._dp, 56._dp,-56._dp,-80._dp, 16._dp, 64._dp, 64._dp, 16._dp,-80._dp,&
2546 : & -56._dp, 56._dp, 56._dp,-56._dp,-80._dp, 16._dp, 64._dp, 64._dp, 16._dp,-80._dp,&
2547 : & -14._dp, 14._dp, 14._dp,-14._dp,-20._dp, 4._dp, 16._dp, 16._dp, 4._dp,-20._dp,&
2548 0 : & 70._dp,-70._dp,-70._dp, 70._dp,100._dp,-20._dp,-80._dp,-80._dp,-20._dp,100._dp/),(/10,10/))
2549 0 : app(1,:,:)=a2pp(:,:)/1225._dp
2550 0 : app(2,:,:)=zero
2551 :
2552 : b2pp(:,:)=RESHAPE((/1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2553 : & 0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2554 : & 0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2555 : & 0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2556 : & 0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2557 : & 0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,&
2558 : & 0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,&
2559 : & 0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,&
2560 : & 0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,&
2561 0 : & 0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp/),(/10,10/))
2562 0 : bpp(0,:,:)=b2pp(:,:)
2563 : b2pp(:,:)=RESHAPE((/ 49._dp, 98._dp, 98._dp, 0._dp, 30._dp, 36._dp, 27._dp, 12._dp, 0._dp, 0._dp,&
2564 : & 98._dp, 49._dp, 0._dp, 98._dp, 40._dp, 2._dp, 6._dp, 25._dp, 32._dp, 0._dp,&
2565 : & 98._dp, 0._dp, 49._dp, 98._dp, 0._dp, 32._dp, 25._dp, 6._dp, 2._dp, 40._dp,&
2566 : & 0._dp, 98._dp, 98._dp, 49._dp, 0._dp, 0._dp, 12._dp, 27._dp, 36._dp, 30._dp,&
2567 : & 30._dp, 40._dp, 0._dp, 0._dp,100._dp,120._dp, 60._dp, 0._dp, 0._dp, 0._dp,&
2568 : & 36._dp, 2._dp, 32._dp, 0._dp,120._dp, 4._dp, 48._dp,108._dp, 0._dp, 0._dp,&
2569 : & 27._dp, 6._dp, 25._dp, 12._dp, 60._dp, 48._dp, 64._dp, 0._dp,108._dp, 0._dp,&
2570 : & 12._dp, 25._dp, 6._dp, 27._dp, 0._dp,108._dp, 0._dp, 64._dp, 48._dp, 60._dp,&
2571 : & 0._dp, 32._dp, 2._dp, 36._dp, 0._dp, 0._dp,108._dp, 48._dp, 4._dp,120._dp,&
2572 0 : & 0._dp, 0._dp, 40._dp, 30._dp, 0._dp, 0._dp, 0._dp, 60._dp,120._dp,100._dp/),(/10,10/))
2573 0 : bpp(1,:,:)=b2pp(:,:)/1225._dp
2574 0 : write(std_out,*) "warning: this test is only valid if f4of2_sla=0"
2575 : end if ! lcor=2
2576 0 : if (lcor == 3) then
2577 : ! app(0,:,:)=one
2578 : ! a2pp(:,:)=RESHAPE((/ 100.0/1225.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2579 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2580 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2581 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2582 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2583 : !& 0 ,0, 0, 0, 0, 100.0/1225.0, 0, 0, 0, 0, 0, 0, 0, 0,&
2584 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2585 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2586 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2587 : !& 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2588 : !& 0 ,0, 1,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0,&
2589 : !& 0 ,0,-1, 1, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0,&
2590 : !& 0 ,0,-1, 1, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0,&
2591 : !& 0 ,0, 1,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0,&/),(/14,14/))
2592 : b2pp(:,:)=RESHAPE((/ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2593 : & 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2594 : & 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2595 : & 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2596 : & 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2597 : & 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2598 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2599 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2600 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2601 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,&
2602 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,&
2603 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,&
2604 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,&
2605 0 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/),(/14,14/))
2606 0 : bpp(0,:,:)=b2pp(:,:)
2607 : b2pp(:,:)=RESHAPE((/100.0, 120.0, 60.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2608 : & 120.0, 4.0, 48.0, 108.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2609 : & 60.0, 48.0, 64.0, 0.0, 108.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2610 : & 0.0, 108.0, 0.0, 64.0, 48.0, 60.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2611 : & 0.0, 0.0, 108.0, 48.0, 4.0, 120.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2612 : & 0.0, 0.0, 0.0, 60.0, 120.0, 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2613 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2614 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2615 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2616 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2617 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2618 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2619 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,&
2620 0 : & 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/),(/14,14/))
2621 0 : bpp(1,:,:)=b2pp(:,:)/1225._dp
2622 0 : write(std_out,*) "warning: only 5/2 5/2 elements are given and only for exchange and F2 !"
2623 : ! app(1,:,:)=a2pp(:,:)
2624 : ! app(1,:,:)=app(1,:,:)/25
2625 : !
2626 : ! bpp(0,:,:)=one
2627 : ! b2pp(:,:)=RESHAPE((/0,0,1,2,3,4,&
2628 : !& 0,0,4,3,2,1,&
2629 : !& 1,4,1,2,2,0,&
2630 : !& 2,3,2,1,0,2,&
2631 : !& 3,2,2,0,1,2,&
2632 : !& 4,1,0,2,2,1 /),(/6,6/))
2633 : ! bpp(1,:,:)=b2pp(:,:)
2634 : ! bpp(1,:,:)=bpp(1,:,:)/25
2635 : end if ! lcor=3
2636 :
2637 0 : do m2=1,tndim
2638 0 : do m1=1,tndim
2639 0 : udens(m1,m2) = sum(fk(:)*app(:,m1,m2))
2640 0 : jdens(m1,m2) = sum(fk(:)*bpp(:,m1,m2))
2641 : !write(6,*) kk,m1,m2
2642 : !write(6,*) "--",fk(kk),aklmlmp(kk,m1,m2)
2643 : !write(6,*) "--",fk(kk),bklmlmp(kk,m1,m2)
2644 : !udens(m1,m2)=udens(m1,m2)+fk(kk)*app(kk,m1,m2)
2645 : !jdens(m1,m2)=jdens(m1,m2)+fk(kk)*bpp(kk,m1,m2)
2646 : end do ! m1
2647 : end do ! m2
2648 0 : write(message,'(2x,a,3x,14f10.4)') " Direct Interaction Matrix from Inglis tables (in the JMJ basis) "
2649 0 : call wrtout(std_out,message,'COLL')
2650 0 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,tndim,1)
2651 0 : call wrtout(std_out,message,'COLL')
2652 0 : do m1=1,tndim
2653 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2),m2=1,tndim,1)
2654 0 : call wrtout(std_out,message,'COLL')
2655 : end do ! m1
2656 :
2657 0 : write(message,'(a,2x,a,3x,14f10.4)') ch10," Exchange Interaction Matrix from Inglis tables (in the JMJ basis) "
2658 0 : call wrtout(std_out,message,'COLL')
2659 0 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,tndim,1)
2660 0 : call wrtout(std_out,message,'COLL')
2661 0 : do m1=1,tndim
2662 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(jdens(m1,m2),m2=1,tndim,1)
2663 0 : call wrtout(std_out,message,'COLL')
2664 : end do ! m1
2665 :
2666 0 : write(message,'(a,2x,a,3x,14f10.4)') ch10, " Density Density interactions from Inglis tables (in the JMJ basis) "
2667 0 : call wrtout(std_out,message,'COLL')
2668 0 : write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,tndim,1)
2669 0 : call wrtout(std_out,message,'COLL')
2670 0 : do m1=1,tndim
2671 0 : write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2)-jdens(m1,m2),m2=1,tndim,1)
2672 0 : call wrtout(std_out,message,'COLL')
2673 : end do ! m1
2674 :
2675 :
2676 0 : ABI_FREE(jdens)
2677 0 : ABI_FREE(udens)
2678 0 : ABI_FREE(app)
2679 0 : ABI_FREE(bpp)
2680 0 : ABI_FREE(a2pp)
2681 0 : ABI_FREE(b2pp)
2682 :
2683 0 : end subroutine udens_inglis_hu
2684 : !!***
2685 :
2686 0 : END MODULE m_hu
2687 : !!***
|
