Branch data 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 : : !!***
|
