TABLE OF CONTENTS


ABINIT/m_hu [ Modules ]

[ Top ] [ Modules ]

NAME

  m_hu

FUNCTION

COPYRIGHT

 Copyright (C) 2006-2022 ABINIT group (BAmadon)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

SOURCE

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 m_abicore
29 
30  use defs_basis
31  use m_pawtab, only : pawtab_type
32 
33  implicit none
34 
35  private
36 
37  public :: init_hu
38  public :: copy_hu
39  public :: destroy_hu
40 ! public :: qmc_hu
41  public :: print_hu
42  public :: vee2udens_hu
43  public :: rotatevee_hu
44  public :: printvee_hu
45  public :: vee2udensatom_hu
46  public :: vee_slm2ylm_hu
47  public :: vee_ndim2tndim_hu
48  public :: vee_ndim2tndim_hu_r
49  public :: udens_slatercondon_hu
50  public :: udens_inglis_hu

ABINIT/udens_inglis_hu [ Functions ]

[ Top ] [ Functions ]

NAME

 udens_inglis_hu

FUNCTION

 For a given angular momentum l and Slater integrals, give the
 density density interactions U(m,m') and J(m,m') from Inglis tables
 in JMJ Basis

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (BA)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  lcor= angular momentum
  fk(lcor+1)= Slater integrals

SIDE EFFECTS

NOTES

SOURCE

2250 subroutine udens_inglis_hu(fk,lcor)
2251 
2252  use defs_basis
2253  use m_errors
2254  use m_abicore
2255  implicit none
2256 
2257 !Arguments ---------------------------------------------
2258 !scalars
2259  integer,intent(in) :: lcor
2260  real(dp), intent(in) :: fk(0:lcor)
2261 
2262 !Local variables ---------------------------------------
2263 !scalars
2264  character(len=500) :: message
2265  integer :: m1,m2,kk,tndim
2266 !arrays
2267  real(dp), allocatable :: udens(:,:)
2268  real(dp), allocatable :: jdens(:,:)
2269  real(dp),allocatable :: app(:,:,:)
2270  real(dp),allocatable :: bpp(:,:,:)
2271  real(dp),allocatable :: a2pp(:,:)
2272  real(dp),allocatable :: b2pp(:,:)
2273 
2274 !*********************************************************************
2275  tndim=2*(2*lcor+1)
2276  ABI_MALLOC(app,(0:lcor,tndim,tndim))
2277  ABI_MALLOC(bpp,(0:lcor,tndim,tndim))
2278  ABI_MALLOC(a2pp,(tndim,tndim))
2279  ABI_MALLOC(b2pp,(tndim,tndim))
2280 
2281  ABI_MALLOC(udens,(tndim,tndim))
2282  ABI_MALLOC(jdens,(tndim,tndim))
2283  udens=zero
2284  jdens=zero
2285  a2pp=zero
2286  b2pp=zero
2287  app=zero
2288  bpp=zero
2289  if(lcor==1) then
2290    app(0,:,:)=one
2291    a2pp(:,:)=RESHAPE((/0._dp,0._dp, 0._dp, 0._dp, 0._dp, 0._dp,&
2292 &                      0._dp,0._dp, 0._dp, 0._dp, 0._dp, 0._dp,&
2293 &                      0._dp,0._dp, 1._dp,-1._dp,-1._dp, 1._dp,&
2294 &                      0._dp,0._dp,-1._dp, 1._dp, 1._dp,-1._dp,&
2295 &                      0._dp,0._dp,-1._dp, 1._dp, 1._dp,-1._dp,&
2296 &                      0._dp,0._dp, 1._dp,-1._dp,-1._dp, 1._dp/),(/6,6/))
2297    app(1,:,:)=a2pp(:,:)
2298    app(1,:,:)=app(1,:,:)/25._dp
2299 
2300    b2pp(:,:)=RESHAPE((/1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2301 &                      0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,&
2302 &                      0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,&
2303 &                      0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,&
2304 &                      0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,&
2305 &                      0._dp,0._dp,0._dp,0._dp,0._dp,1._dp /),(/6,6/))
2306    bpp(0,:,:)=b2pp(:,:)
2307    b2pp(:,:)=RESHAPE((/0._dp,0._dp,1._dp,2._dp,3._dp,4._dp,&
2308 &                      0._dp,0._dp,4._dp,3._dp,2._dp,1._dp,&
2309 &                      1._dp,4._dp,1._dp,2._dp,2._dp,0._dp,&
2310 &                      2._dp,3._dp,2._dp,1._dp,0._dp,2._dp,&
2311 &                      3._dp,2._dp,2._dp,0._dp,1._dp,2._dp,&
2312 &                      4._dp,1._dp,0._dp,2._dp,2._dp,1._dp /),(/6,6/))
2313    bpp(1,:,:)=b2pp(:,:)
2314    bpp(1,:,:)=bpp(1,:,:)/25._dp
2315  endif
2316  if(lcor==2) then
2317    app(0,:,:)=one
2318    a2pp(:,:)=RESHAPE((/ 49._dp,-49._dp,-49._dp, 49._dp, 70._dp,-14._dp,-56._dp,-56._dp,-14._dp, 70._dp,&
2319 &                      -49._dp, 49._dp, 49._dp,-49._dp,-70._dp, 14._dp, 56._dp, 56._dp, 14._dp,-70._dp,&
2320 &                      -49._dp, 49._dp, 49._dp,-49._dp,-70._dp, 14._dp, 56._dp, 56._dp, 14._dp,-70._dp,&
2321 &                       49._dp,-49._dp,-49._dp, 49._dp, 70._dp,-14._dp,-56._dp,-56._dp,-14._dp, 70._dp,&
2322 &                       70._dp,-70._dp,-70._dp, 70._dp,100._dp,-20._dp,-80._dp,-80._dp,-20._dp,100._dp,&
2323 &                      -14._dp, 14._dp, 14._dp,-14._dp,-20._dp,  4._dp, 16._dp, 16._dp,  4._dp,-20._dp,&
2324 &                      -56._dp, 56._dp, 56._dp,-56._dp,-80._dp, 16._dp, 64._dp, 64._dp, 16._dp,-80._dp,&
2325 &                      -56._dp, 56._dp, 56._dp,-56._dp,-80._dp, 16._dp, 64._dp, 64._dp, 16._dp,-80._dp,&
2326 &                      -14._dp, 14._dp, 14._dp,-14._dp,-20._dp,  4._dp, 16._dp, 16._dp,  4._dp,-20._dp,&
2327 &                       70._dp,-70._dp,-70._dp, 70._dp,100._dp,-20._dp,-80._dp,-80._dp,-20._dp,100._dp/),(/10,10/))
2328    app(1,:,:)=a2pp(:,:)/1225._dp
2329    app(2,:,:)=zero
2330 
2331    b2pp(:,:)=RESHAPE((/1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2332 &                      0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2333 &                      0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2334 &                      0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2335 &                      0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,0._dp,&
2336 &                      0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,0._dp,&
2337 &                      0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,0._dp,&
2338 &                      0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,0._dp,&
2339 &                      0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp,0._dp,&
2340 &                      0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,0._dp,1._dp/),(/10,10/))
2341    bpp(0,:,:)=b2pp(:,:)
2342    b2pp(:,:)=RESHAPE((/ 49._dp, 98._dp, 98._dp,  0._dp, 30._dp, 36._dp, 27._dp, 12._dp,  0._dp,  0._dp,&
2343 &                       98._dp, 49._dp,  0._dp, 98._dp, 40._dp,  2._dp,  6._dp, 25._dp, 32._dp,  0._dp,&
2344 &                       98._dp,  0._dp, 49._dp, 98._dp,  0._dp, 32._dp, 25._dp,  6._dp,  2._dp, 40._dp,&
2345 &                        0._dp, 98._dp, 98._dp, 49._dp,  0._dp,  0._dp, 12._dp, 27._dp, 36._dp, 30._dp,&
2346 &                       30._dp, 40._dp,  0._dp,  0._dp,100._dp,120._dp, 60._dp,  0._dp,  0._dp,  0._dp,&
2347 &                       36._dp,  2._dp, 32._dp,  0._dp,120._dp,  4._dp, 48._dp,108._dp,  0._dp,  0._dp,&
2348 &                       27._dp,  6._dp, 25._dp, 12._dp, 60._dp, 48._dp, 64._dp,  0._dp,108._dp,  0._dp,&
2349 &                       12._dp, 25._dp,  6._dp, 27._dp,  0._dp,108._dp,  0._dp, 64._dp, 48._dp, 60._dp,&
2350 &                        0._dp, 32._dp,  2._dp, 36._dp,  0._dp,  0._dp,108._dp, 48._dp,  4._dp,120._dp,&
2351 &                        0._dp,  0._dp, 40._dp, 30._dp,  0._dp,  0._dp,  0._dp, 60._dp,120._dp,100._dp/),(/10,10/))
2352    bpp(1,:,:)=b2pp(:,:)/1225._dp
2353    write(std_out,*) "warning: this test is only valid if f4of2_sla=0"
2354  endif
2355  if(lcor==3) then
2356 !   app(0,:,:)=one
2357 !   a2pp(:,:)=RESHAPE((/ 100.0/1225.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2358 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2359 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2360 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2361 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2362 !&                        0 ,0, 0, 0, 0, 100.0/1225.0, 0, 0, 0, 0, 0, 0, 0, 0,&
2363 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2364 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2365 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2366 !&                        0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,&
2367 !&                        0 ,0, 1,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0,&
2368 !&                        0 ,0,-1, 1, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0,&
2369 !&                        0 ,0,-1, 1, 1,-1, 0, 0, 0, 0, 0, 0, 0, 0,&
2370 !&                        0 ,0, 1,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0,&/),(/14,14/))
2371    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,&
2372 &                        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,&
2373 &                        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,&
2374 &                        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,&
2375 &                        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,&
2376 &                        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,&
2377 &                        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,&
2378 &                        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,&
2379 &                        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,&
2380 &                        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,&
2381 &                        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,&
2382 &                        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,&
2383 &                        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,&
2384 &                        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/))
2385    bpp(0,:,:)=b2pp(:,:)
2386    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,&
2387 &                      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,&
2388 &                       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,&
2389 &                        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,&
2390 &                        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,&
2391 &                        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,&
2392 &                        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,&
2393 &                        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,&
2394 &                        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,&
2395 &                        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,&
2396 &                        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,&
2397 &                        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,&
2398 &                        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,&
2399 &                        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/))
2400    bpp(1,:,:)=b2pp(:,:)/1225._dp
2401    write(std_out,*) "warning: only 5/2 5/2 elements are given and only for exchange and F2 !"
2402 !   app(1,:,:)=a2pp(:,:)
2403 !   app(1,:,:)=app(1,:,:)/25
2404 !
2405 !   bpp(0,:,:)=one
2406 !   b2pp(:,:)=RESHAPE((/0,0,1,2,3,4,&
2407 !&                      0,0,4,3,2,1,&
2408 !&                      1,4,1,2,2,0,&
2409 !&                      2,3,2,1,0,2,&
2410 !&                      3,2,2,0,1,2,&
2411 !&                      4,1,0,2,2,1 /),(/6,6/))
2412 !   bpp(1,:,:)=b2pp(:,:)
2413 !   bpp(1,:,:)=bpp(1,:,:)/25
2414  endif
2415 
2416  do kk=0,lcor
2417    do m1=1,tndim
2418      do m2=1,tndim
2419        !write(6,*) kk,m1,m2
2420        !write(6,*) "--",fk(kk),aklmlmp(kk,m1,m2)
2421        !write(6,*) "--",fk(kk),bklmlmp(kk,m1,m2)
2422        udens(m1,m2)=udens(m1,m2)+fk(kk)*app(kk,m1,m2)
2423        jdens(m1,m2)=jdens(m1,m2)+fk(kk)*bpp(kk,m1,m2)
2424      enddo
2425    enddo
2426  enddo
2427  write(message,'(2x,a,3x,14f10.4)') " Direct Interaction Matrix from Inglis tables (in the JMJ basis) "
2428  call wrtout(std_out,  message,'COLL')
2429  write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,tndim,1)
2430  call wrtout(std_out,  message,'COLL')
2431  do m1=1,tndim
2432    write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2),m2=1,tndim,1)
2433    call wrtout(std_out,  message,'COLL')
2434  enddo
2435 
2436  write(message,'(a,2x,a,3x,14f10.4)') ch10," Exchange Interaction Matrix from Inglis tables (in the JMJ basis) "
2437  call wrtout(std_out,  message,'COLL')
2438  write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,tndim,1)
2439  call wrtout(std_out,  message,'COLL')
2440  do m1=1,tndim
2441    write(message,'(2x,i4,3x,14f10.4)') m1,(jdens(m1,m2),m2=1,tndim,1)
2442    call wrtout(std_out,  message,'COLL')
2443  enddo
2444 
2445  write(message,'(a,2x,a,3x,14f10.4)') ch10, " Density Density interactions from Inglis tables (in the JMJ basis) "
2446  call wrtout(std_out,  message,'COLL')
2447  write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,tndim,1)
2448  call wrtout(std_out,  message,'COLL')
2449  do m1=1,tndim
2450    write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2)-jdens(m1,m2),m2=1,tndim,1)
2451    call wrtout(std_out,  message,'COLL')
2452  enddo
2453 
2454 
2455  ABI_FREE(jdens)
2456  ABI_FREE(udens)
2457  ABI_FREE(app)
2458  ABI_FREE(bpp)
2459  ABI_FREE(a2pp)
2460  ABI_FREE(b2pp)
2461 
2462  end subroutine udens_inglis_hu

ABINIT/udens_slatercondon_hu [ Functions ]

[ Top ] [ Functions ]

NAME

 udens_slatercondon_hu

FUNCTION

 For a given angular momentum l and Slater integrals, give the
 density density interactions U(m,m') and J(m,m') from Slater and
 Condon tables

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (BA)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  lcor= angular momentum
  fk(lcor+1)= Slater integrals

SIDE EFFECTS

NOTES

SOURCE

1981 subroutine udens_slatercondon_hu(fk,lcor)
1982 
1983  use defs_basis
1984  use m_errors
1985  use m_abicore
1986  implicit none
1987 
1988 !Arguments ---------------------------------------------
1989 !scalars
1990  integer,intent(in) :: lcor
1991  real(dp), intent(in) :: fk(0:lcor)
1992 
1993 !Local variables ---------------------------------------
1994 !scalars
1995  character(len=500) :: message
1996  integer :: m1,m2,kk
1997 !arrays
1998  real(dp), allocatable :: aklmlmp(:,:,:)
1999  real(dp), allocatable :: bklmlmp(:,:,:)
2000  real(dp), allocatable :: udens(:,:)
2001  real(dp), allocatable :: jdens(:,:)
2002 
2003 !*********************************************************************
2004  ABI_MALLOC(aklmlmp,(0:lcor,-lcor:lcor,-lcor:lcor)) ! k,m,m'
2005  ABI_MALLOC(bklmlmp,(0:lcor,-lcor:lcor,-lcor:lcor)) ! k,m,m'
2006  ABI_MALLOC(udens,(-lcor:lcor,-lcor:lcor)) ! m,m'
2007  ABI_MALLOC(jdens,(-lcor:lcor,-lcor:lcor)) ! m,m'
2008 ! k=2*(lcor)
2009  aklmlmp=zero
2010  bklmlmp=zero
2011  udens=zero
2012  jdens=zero
2013  if(lcor==0) then
2014    aklmlmp(0,0,0)=1
2015 !
2016    bklmlmp(0,0,0)=0
2017  else if(lcor==1) then
2018    aklmlmp(0, :, :)=1
2019    aklmlmp(1, 1, 1)= one/25._dp
2020    aklmlmp(1,-1,-1)= one/25._dp
2021    aklmlmp(1,-1, 1)= one/25._dp
2022    aklmlmp(1, 1,-1)= one/25._dp
2023    aklmlmp(1, 1, 0)=-two/25._dp
2024    aklmlmp(1,-1, 0)=-two/25._dp
2025    aklmlmp(1, 0,-1)=-two/25._dp
2026    aklmlmp(1, 0, 1)=-two/25._dp
2027    aklmlmp(1, 0, 0)= four/25._dp
2028 !
2029    bklmlmp(0, 1, 1)= one
2030    bklmlmp(0,-1,-1)= one
2031    bklmlmp(0, 0, 0)= one
2032    bklmlmp(1, 1, 1)= one/25._dp
2033    bklmlmp(1,-1,-1)= one/25._dp
2034    bklmlmp(1,-1,+1)= six/25._dp
2035    bklmlmp(1,+1,-1)= six/25._dp
2036    bklmlmp(1,+1, 0)= three/25._dp
2037    bklmlmp(1,-1, 0)= three/25._dp
2038    bklmlmp(1, 0,-1)= three/25._dp
2039    bklmlmp(1, 0, 1)= three/25._dp
2040    bklmlmp(1, 0, 0)= four/25._dp
2041  else if(lcor==2) then
2042    aklmlmp(0, :, :)=1
2043    aklmlmp(1, 2, 2)=  four / 49._dp
2044    aklmlmp(1, 2, 1)=  -two / 49._dp
2045    aklmlmp(1, 2, 0)= -four / 49._dp
2046    aklmlmp(1, 2,-1)=  -two / 49._dp
2047    aklmlmp(1, 2,-2)=  four / 49._dp
2048    aklmlmp(1, 1, 2)=  -two / 49._dp
2049    aklmlmp(1, 1, 1)=   one / 49._dp
2050    aklmlmp(1, 1, 0)=   two / 49._dp
2051    aklmlmp(1, 1,-1)=   one / 49._dp
2052    aklmlmp(1, 1,-2)=  -two / 49._dp
2053    aklmlmp(1, 0, 2)= -four / 49._dp
2054    aklmlmp(1, 0, 1)=   two / 49._dp
2055    aklmlmp(1, 0, 0)=  four / 49._dp
2056    aklmlmp(1, 0,-1)=   two / 49._dp
2057    aklmlmp(1, 0,-2)= -four / 49._dp
2058    aklmlmp(1,-1, 2)=  -two / 49._dp
2059    aklmlmp(1,-1, 1)=   one / 49._dp
2060    aklmlmp(1,-1, 0)=   two / 49._dp
2061    aklmlmp(1,-1,-1)=   one / 49._dp
2062    aklmlmp(1,-1,-2)=  -two / 49._dp
2063    aklmlmp(1,-2, 2)=  four / 49._dp
2064    aklmlmp(1,-2, 1)=  -two / 49._dp
2065    aklmlmp(1,-2, 0)= -four / 49._dp
2066    aklmlmp(1,-2,-1)=  -two / 49._dp
2067    aklmlmp(1,-2,-2)=  four / 49._dp
2068 
2069    aklmlmp(2, 2, 2)=    one  / 441._dp
2070    aklmlmp(2, 2, 1)=  -four  / 441._dp
2071    aklmlmp(2, 2, 0)=    six  / 441._dp
2072    aklmlmp(2, 2,-1)=  -four  / 441._dp
2073    aklmlmp(2, 2,-2)=    one  / 441._dp
2074    aklmlmp(2, 1, 2)=  -four  / 441._dp
2075    aklmlmp(2, 1, 1)=  16._dp / 441._dp
2076    aklmlmp(2, 1, 0)= -24._dp / 441._dp
2077    aklmlmp(2, 1,-1)=  16._dp / 441._dp
2078    aklmlmp(2, 1,-2)=  -four  / 441._dp
2079    aklmlmp(2, 0, 2)=    six  / 441._dp
2080    aklmlmp(2, 0, 1)= -24._dp / 441._dp
2081    aklmlmp(2, 0, 0)=  36._dp / 441._dp
2082    aklmlmp(2, 0,-1)= -24._dp / 441._dp
2083    aklmlmp(2, 0,-2)=    six  / 441._dp
2084    aklmlmp(2,-1, 2)=  -four  / 441._dp
2085    aklmlmp(2,-1, 1)=  16._dp / 441._dp
2086    aklmlmp(2,-1, 0)= -24._dp / 441._dp
2087    aklmlmp(2,-1,-1)=  16._dp / 441._dp
2088    aklmlmp(2,-1,-2)=  -four  / 441._dp
2089    aklmlmp(2,-2, 2)=    one  / 441._dp
2090    aklmlmp(2,-2, 1)=  -four  / 441._dp
2091    aklmlmp(2,-2, 0)=    six  / 441._dp
2092    aklmlmp(2,-2,-1)=  -four  / 441._dp
2093    aklmlmp(2,-2,-2)=    one  / 441._dp
2094    !do m1=lcor,-lcor,-1
2095    ! do m2=lcor,-lcor,-1
2096    !   write(6,*) m1,m2,aklmlmp(2,m1,m2)
2097    ! enddo
2098    !enddo
2099 
2100    !write(message,'(2x,a,3x,14f10.4)') " Slater aklmlmp(2,:,:)"
2101    !call wrtout(std_out,  message,'COLL')
2102    !write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2103    !call wrtout(std_out,  message,'COLL')
2104    !do m1=-lcor,lcor,1
2105    !  write(message,'(2x,i4,3x,14f10.4)') m1,(aklmlmp(2,m1,m2),m2=-lcor,lcor,1)
2106    !  call wrtout(std_out,  message,'COLL')
2107    !enddo
2108 
2109    bklmlmp(0, 2, 2)=1
2110    bklmlmp(0, 1, 1)=1
2111    bklmlmp(0, 0, 0)=1
2112    bklmlmp(0,-1,-1)=1
2113    bklmlmp(0,-2,-2)=1
2114    bklmlmp(1, 2, 2)= four / 49._dp
2115    bklmlmp(1, 2, 1)=  six / 49._dp
2116    bklmlmp(1, 2, 0)= four / 49._dp
2117    bklmlmp(1, 2,-1)=  zero
2118    bklmlmp(1, 2,-2)=  zero
2119    bklmlmp(1, 1, 2)=  six / 49._dp
2120    bklmlmp(1, 1, 1)=  one / 49._dp
2121    bklmlmp(1, 1, 0)=  one / 49._dp
2122    bklmlmp(1, 1,-1)=  six / 49._dp
2123    bklmlmp(1, 1,-2)=  zero
2124    bklmlmp(1, 0, 2)= four / 49._dp
2125    bklmlmp(1, 0, 1)=  one / 49._dp
2126    bklmlmp(1, 0, 0)= four / 49._dp
2127    bklmlmp(1, 0,-1)=  one / 49._dp
2128    bklmlmp(1, 0,-2)= four / 49._dp
2129    bklmlmp(1,-1, 2)=  zero
2130    bklmlmp(1,-1, 1)=  six / 49._dp
2131    bklmlmp(1,-1, 0)=  one / 49._dp
2132    bklmlmp(1,-1,-1)=  one / 49._dp
2133    bklmlmp(1,-1,-2)=  six / 49._dp
2134    bklmlmp(1,-2, 2)=  zero
2135    bklmlmp(1,-2, 1)=  zero
2136    bklmlmp(1,-2, 0)= four / 49._dp
2137    bklmlmp(1,-2,-1)=  six / 49._dp
2138    bklmlmp(1,-2,-2)= four / 49._dp
2139 
2140    bklmlmp(2, 2, 2)=   one  / 441._dp
2141    bklmlmp(2, 2, 1)=  five  / 441._dp
2142    bklmlmp(2, 2, 0)= 15._dp / 441._dp
2143    bklmlmp(2, 2,-1)= 35._dp / 441._dp
2144    bklmlmp(2, 2,-2)= 70._dp / 441._dp
2145    bklmlmp(2, 1, 2)=  five  / 441._dp
2146    bklmlmp(2, 1, 1)= 16._dp / 441._dp
2147    bklmlmp(2, 1, 0)= 30._dp / 441._dp
2148    bklmlmp(2, 1,-1)= 40._dp / 441._dp
2149    bklmlmp(2, 1,-2)= 35._dp / 441._dp
2150    bklmlmp(2, 0, 2)= 15._dp / 441._dp
2151    bklmlmp(2, 0, 1)= 30._dp / 441._dp
2152    bklmlmp(2, 0, 0)= 36._dp / 441._dp
2153    bklmlmp(2, 0,-1)= 30._dp / 441._dp
2154    bklmlmp(2, 0,-2)= 15._dp / 441._dp
2155    bklmlmp(2,-1, 2)= 35._dp / 441._dp
2156    bklmlmp(2,-1, 1)= 40._dp / 441._dp
2157    bklmlmp(2,-1, 0)= 30._dp / 441._dp
2158    bklmlmp(2,-1,-1)= 16._dp / 441._dp
2159    bklmlmp(2,-1,-2)=  five  / 441._dp
2160    bklmlmp(2,-2, 2)= 70._dp / 441._dp
2161    bklmlmp(2,-2, 1)= 35._dp / 441._dp
2162    bklmlmp(2,-2, 0)= 15._dp / 441._dp
2163    bklmlmp(2,-2,-1)=  five  / 441._dp
2164    bklmlmp(2,-2,-2)=   one  / 441._dp
2165 
2166    !write(message,'(2x,a,3x,14f10.4)') " Slater bklmlmp(2,:,:)"
2167    !call wrtout(std_out,  message,'COLL')
2168    !write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2169    !call wrtout(std_out,  message,'COLL')
2170    !do m1=-lcor,lcor,1
2171    !  write(message,'(2x,i4,3x,14f10.4)') m1,(bklmlmp(2,m1,m2),m2=-lcor,lcor,1)
2172    !  call wrtout(std_out,  message,'COLL')
2173    !enddo
2174          ! if f4of2_sla: these data agree with the explicit calculation
2175  else if(lcor==3) then
2176  endif
2177 
2178  do kk=0,lcor
2179    do m1=-lcor,lcor,1
2180      do m2=-lcor,lcor,1
2181        !write(6,*) kk,m1,m2
2182        !write(6,*) "--",fk(kk),aklmlmp(kk,m1,m2)
2183        !write(6,*) "--",fk(kk),bklmlmp(kk,m1,m2)
2184        udens(m1,m2)=udens(m1,m2)+fk(kk)*aklmlmp(kk,m1,m2)
2185        jdens(m1,m2)=jdens(m1,m2)+fk(kk)*bklmlmp(kk,m1,m2)
2186      enddo
2187    enddo
2188  enddo
2189  write(message,'(2x,a,3x,14f10.4)') " Direct Interaction Matrix from Slater tables (in the Ylm basis) "
2190  call wrtout(std_out,  message,'COLL')
2191  write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2192  call wrtout(std_out,  message,'COLL')
2193  do m1=-lcor,lcor,1
2194    write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2),m2=-lcor,lcor,1)
2195    call wrtout(std_out,  message,'COLL')
2196  enddo
2197 
2198  write(message,'(a,2x,a,3x,14f10.4)') ch10," Exchange Interaction Matrix from Slater tables (in the Ylm basis) "
2199  call wrtout(std_out,  message,'COLL')
2200  write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2201  call wrtout(std_out,  message,'COLL')
2202  do m1=-lcor,lcor,1
2203    write(message,'(2x,i4,3x,14f10.4)') m1,(jdens(m1,m2),m2=-lcor,lcor,1)
2204    call wrtout(std_out,  message,'COLL')
2205  enddo
2206 
2207  write(message,'(a,2x,a,3x,14f10.4)') ch10," Density density Interaction Matrix from Slater tables (in the Ylm basis) "
2208  call wrtout(std_out,  message,'COLL')
2209  write(message,'(2x,4x,14(2x,i8))') (m1,m1=-lcor,lcor,1)
2210  call wrtout(std_out,  message,'COLL')
2211  do m1=-lcor,lcor,1
2212    write(message,'(2x,i4,3x,14f10.4)') m1,(udens(m1,m2)-jdens(m1,m2),m2=-lcor,lcor,1)
2213    call wrtout(std_out,  message,'COLL')
2214  enddo
2215 
2216  ABI_FREE(jdens)
2217  ABI_FREE(udens)
2218  ABI_FREE(aklmlmp)
2219  ABI_FREE(bklmlmp)
2220 
2221  end subroutine udens_slatercondon_hu

ABINIT/vee_ylm2jmj_hu [ Functions ]

[ Top ] [ Functions ]

NAME

 vee_ylm2jmj_hu

FUNCTION

 For a given angular momentum lcor, change a matrix  of dimension [2(2*lcor+1)]**4
 from the Ylm basis to the J,M_J basis if option==1

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (BA)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  lcor= angular momentum
  option=  1 matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis
           2 matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis

SIDE EFFECTS

  mat_mlms= Input/Ouput matrix in the Ylm basis, size of the matrix is (2*lcor+1,2*lcor+1,ndij)
  mat_jmj= Input/Output matrix in the J,M_J basis, size is 2*(2*lcor+1),2*(2*lcor+1)

NOTES

  usefull only in ndij==4

SOURCE

1828 subroutine vee_ylm2jmj_hu(lcor,mat_inp_c,mat_out_c,option)
1829 
1830  use defs_basis
1831  use m_errors
1832  use m_abicore
1833  implicit none
1834 
1835 !Arguments ---------------------------------------------
1836 !scalars
1837  integer,intent(in) :: lcor,option
1838 !arrays
1839  complex(dpc),intent(in) :: mat_inp_c(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
1840  complex(dpc),intent(out) :: mat_out_c(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
1841 
1842 !Local variables ---------------------------------------
1843 !scalars
1844  integer :: ii,im,jc1,jj,jm,ll,ml1,ms1,gg,hh,gm,hm
1845  real(dp),parameter :: invsqrt2=one/sqrt2
1846  real(dp) :: invsqrt2lp1,xj,xmj
1847  complex(dpc) :: tmp2
1848  character(len=500) :: message
1849 !arrays
1850  integer, allocatable :: ind_msml(:,:)
1851  complex(dpc),allocatable :: mlms2jmj(:,:)
1852 
1853 !*********************************************************************
1854 
1855  if (option/=1.and.option/=2) then
1856    message=' option=/1 and =/2 !'
1857    ABI_BUG(message)
1858  end if
1859 
1860  if(option==1) then
1861    write(message,'(3a)') ch10,&
1862 &   "matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis"
1863    call wrtout(std_out,message)
1864  else if(option==2) then
1865    write(message,'(3a)') ch10,&
1866 &   "matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis"
1867    call wrtout(std_out,message)
1868  end if
1869 
1870 !--------------- Built indices + allocations
1871  ll=lcor
1872  ABI_MALLOC(mlms2jmj,(2*(2*ll+1),2*(2*ll+1)))
1873  mlms2jmj=czero
1874  ABI_MALLOC(ind_msml,(2,-ll:ll))
1875  mlms2jmj=czero
1876  jc1=0
1877  do ms1=1,2
1878    do ml1=-ll,ll
1879      jc1=jc1+1
1880      ind_msml(ms1,ml1)=jc1
1881    end do
1882  end do
1883 
1884 !--------------- built mlms2jmj
1885 !do jj=ll,ll+1    ! the physical value of j are ll-0.5,ll+0.5
1886 !xj(jj)=jj-0.5
1887  if(ll==0)then
1888    message=' ll should not be equal to zero !'
1889    ABI_BUG(message)
1890  end if
1891  jc1=0
1892  invsqrt2lp1=one/sqrt(float(2*lcor+1))
1893  do jj=ll,ll+1
1894    xj=float(jj)-half !  xj is in {ll-0.5, ll+0.5}
1895    do jm=-jj,jj-1
1896      xmj=float(jm)+half  ! xmj is in {-xj,xj}
1897      jc1=jc1+1           ! Global index for JMJ
1898      if(nint(xj+0.5)==ll+1) then  ! if xj=ll+0.5
1899        if(nint(xmj+0.5)==ll+1)  then
1900          mlms2jmj(ind_msml(2,ll),jc1)=1.0   !  J=L+0.5 and m_J=L+0.5
1901        else if(nint(xmj-0.5)==-ll-1) then
1902          mlms2jmj(ind_msml(1,-ll),jc1)=1.0   !  J=L+0.5 and m_J=-L-0.5
1903        else
1904          mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
1905          mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
1906        end if
1907      end if
1908      if(nint(xj+0.5)==ll) then  ! if xj=ll-0.5
1909        mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
1910        mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
1911      end if
1912    end do
1913  end do
1914  write(message,'(3a)') ch10,"Matrix to go from |M_L,M_S> to |J,M_J>"
1915  call wrtout(std_out,message,"COLL")
1916  do im=1,2*(ll*2+1)
1917    write(message,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mlms2jmj(im,jm),jm=1,2*(ll*2+1))
1918    call wrtout(std_out,message,"COLL")
1919  end do
1920 
1921 !--------------- compute change of basis
1922  do jm=1,2*(2*ll+1)
1923    do im=1,2*(2*ll+1)
1924      do hm=1,2*(2*ll+1)
1925        do gm=1,2*(2*ll+1)
1926          tmp2=czero
1927          do gg=1,2*(2*ll+1)
1928            do hh=1,2*(2*ll+1)
1929              do ii=1,2*(2*ll+1)
1930                do jj=1,2*(2*ll+1)
1931                  if(option==1) then
1932                    tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*CONJG(mlms2jmj(ii,im))*(mlms2jmj(jj,jm))&
1933 &                                                  *CONJG(mlms2jmj(gg,gm))*(mlms2jmj(hh,hm))
1934                  else if(option==2) then
1935 !                   tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*CONJG(mlms2jmj(ii,im))*(mlms2jmj(jj,jm))& ! inv=t*
1936 !&                                                  *CONJG(mlms2jmj(gg,gm))*(mlms2jmj(hh,hm)) ! inv=t*
1937                    tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*(mlms2jmj(im,ii))*CONJG(mlms2jmj(jm,jj))& ! inv=t*
1938 &                                                  *(mlms2jmj(gm,gg))*CONJG(mlms2jmj(hm,hh)) ! inv=t*
1939                  end if
1940                end do
1941              end do
1942            end do
1943          end do
1944          mat_out_c(gm,im,hm,jm)=tmp2
1945        end do
1946      end do
1947    end do
1948  end do
1949  ABI_FREE(mlms2jmj)
1950  ABI_FREE(ind_msml)
1951 
1952  end subroutine vee_ylm2jmj_hu

m_hu/copy_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 copy_hu

FUNCTION

  Allocate variables used in type hu_type.

INPUTS

  hu <type(hu_type)>= U interaction

 OUTPUTS
  hu_new <type(hu_type)>= U interaction

SOURCE

312 subroutine copy_hu(ntypat,hu,hu_new)
313 
314  use defs_basis
315  implicit none
316 
317 !Arguments ------------------------------------
318 !type
319  integer, intent(in) :: ntypat
320  type(hu_type), intent(in) :: hu(ntypat)
321  type(hu_type), intent(out) :: hu_new(ntypat)
322 !Local variables ------------------------------------
323  integer :: itypat,ndim
324 !************************************************************************
325 
326  do itypat=1,ntypat
327    hu_new(itypat)%lpawu      = hu(itypat)%lpawu
328    hu_new(itypat)%jmjbasis   = hu(itypat)%jmjbasis
329    hu_new(itypat)%upawu      = hu(itypat)%upawu
330    hu_new(itypat)%jpawu      = hu(itypat)%jpawu
331    hu_new(itypat)%f2_sla     = hu(itypat)%f2_sla
332    hu_new(itypat)%f4of2_sla  = hu(itypat)%f4of2_sla
333    hu_new(itypat)%f6of2_sla  = hu(itypat)%f6of2_sla
334    hu_new(itypat)%jpawu_zero = hu(itypat)%jpawu_zero
335    ndim=2*hu_new(itypat)%lpawu+1
336    ABI_MALLOC(hu_new(itypat)%uqmc,(ndim*(2*ndim-1)))
337    ABI_MALLOC(hu_new(itypat)%udens,(2*ndim,2*ndim))
338    ABI_MALLOC(hu_new(itypat)%vee,(ndim,ndim,ndim,ndim))
339    hu_new(itypat)%vee        = hu(itypat)%vee
340    hu_new(itypat)%udens      = hu(itypat)%udens
341    hu_new(itypat)%uqmc       = hu(itypat)%uqmc
342  enddo ! itypat
343 
344 end subroutine copy_hu

m_hu/destroy_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 destroy_hu

FUNCTION

  deallocate hu

INPUTS

  ntypat = number of species
  hu <type(hu_type)> = data for the interaction in DMFT.

OUTPUT

SOURCE

362 subroutine destroy_hu(hu,ntypat,t2g,x2my2d)
363 
364  use defs_basis
365  use m_crystal, only : crystal_t
366  implicit none
367 
368 !Arguments ------------------------------------
369 !scalars
370  integer, intent(in) :: ntypat
371  type(hu_type),intent(inout) :: hu(ntypat)
372  integer, intent(in) :: t2g,x2my2d
373 
374 !Local variables-------------------------------
375  integer :: itypat
376 
377 ! *********************************************************************
378 
379  do itypat=1,ntypat
380 !  if ( allocated(hu(itypat)%vee) )  deallocate(hu(itypat)%vee)
381   if ( allocated(hu(itypat)%uqmc) )   then
382     ABI_FREE(hu(itypat)%uqmc)
383   end if
384   if ( allocated(hu(itypat)%udens) )   then
385     ABI_FREE(hu(itypat)%udens)
386   end if
387   if(t2g==1.and.hu(itypat)%lpawu==1) then
388     ABI_FREE(hu(itypat)%vee)
389   else if(x2my2d==1.and.hu(itypat)%lpawu==0) then
390     ABI_FREE(hu(itypat)%vee)
391   else
392     hu(itypat)%vee => null()
393   endif
394  enddo
395 
396 end subroutine destroy_hu

m_hu/hu_type [ Types ]

[ Top ] [ m_hu ] [ Types ]

NAME

  hu_type

FUNCTION

  This structured datatype contains interaction matrices for the correlated subspace

SOURCE

64  type, public :: hu_type ! for each typat
65 
66   integer :: lpawu
67 
68   logical :: jmjbasis
69 
70   real(dp) :: upawu    ! => upaw
71 
72   real(dp) :: jpawu    ! => jpaw
73 
74   real(dp) :: f2_sla    ! => f2_sla
75 
76   real(dp) :: f4of2_sla    ! => f4of2_sla
77 
78   real(dp) :: f6of2_sla    ! => f6of2_sla
79 
80   logical :: jpawu_zero  ! true if all jpawu are zero
81                          ! false if one of the jpaw is not zero
82 
83   real(dp), pointer :: vee(:,:,:,:) => null() ! => vee
84 
85   real(dp), allocatable :: uqmc(:)
86 
87   real(dp), allocatable :: udens(:,:)
88 
89  end type hu_type
90 
91 !----------------------------------------------------------------------
92 
93 
94 CONTAINS  !========================================================================================

m_hu/init_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 init_hu

FUNCTION

  Allocate variables used in type hu_type.

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  pawtab <type(pawtab)>=paw related data

 OUTPUTS
  hu <type(hu_type)>= U interaction

SOURCE

113 subroutine init_hu(cryst_struc,pawtab,hu,t2g,x2my2d)
114 
115  use defs_basis
116  use m_crystal, only : crystal_t
117  implicit none
118 
119 !Arguments ------------------------------------
120 !type
121  type(crystal_t),intent(in) :: cryst_struc
122  type(pawtab_type), target, intent(in)  :: pawtab(cryst_struc%ntypat)
123  type(hu_type), intent(inout) :: hu(cryst_struc%ntypat)
124  integer, intent(in) :: t2g,x2my2d
125 !Local variables ------------------------------------
126  integer :: itypat,i,ij,ij1,ij2,j,lpawu,ms,ms1,m,m1,ndim
127  integer :: ns,ns1,n,n1
128  integer, allocatable :: xij(:,:)
129  real(dp) :: xtemp
130  character(len=500) :: message
131 !************************************************************************
132  write(message,'(2a)') ch10,"  == Compute Interactions for DMFT"
133  call wrtout(std_out,message,'COLL')
134 
135  xtemp=zero
136 
137 ! ====================================
138 !  Compute hu(iatom)%uqmc from vee
139 ! ====================================
140  hu(1)%jpawu_zero=.true.
141  do itypat=1,cryst_struc%ntypat
142    hu(itypat)%lpawu=pawtab(itypat)%lpawu
143    hu(itypat)%jmjbasis=.false.
144    if(t2g==1.and.hu(itypat)%lpawu==2) hu(itypat)%lpawu=1
145    if(x2my2d==1.and.hu(itypat)%lpawu==2) hu(itypat)%lpawu=0
146    lpawu=hu(itypat)%lpawu
147    if(lpawu.ne.-1) then
148      hu(itypat)%upawu=pawtab(itypat)%upawu
149      hu(itypat)%jpawu=pawtab(itypat)%jpawu
150 
151    ! The if below are necessary for an unknown reason with the NAG  compiler.
152     if(pawtab(itypat)%f4of2_sla>=0) then
153       hu(itypat)%f4of2_sla = pawtab(itypat)%f4of2_sla
154     else if(pawtab(itypat)%f4of2_sla==0) then
155       hu(itypat)%f4of2_sla = 0_dp
156     else if(pawtab(itypat)%f4of2_sla<0) then
157       hu(itypat)%f4of2_sla = pawtab(itypat)%f4of2_sla
158     endif
159     if(pawtab(itypat)%f6of2_sla>=0) then
160       hu(itypat)%f6of2_sla = pawtab(itypat)%f6of2_sla
161     else if(pawtab(itypat)%f6of2_sla==0) then
162       hu(itypat)%f6of2_sla = 0_dp
163     else if(pawtab(itypat)%f6of2_sla<0) then
164       hu(itypat)%f6of2_sla = pawtab(itypat)%f6of2_sla
165     endif
166 
167 !    This is copied from pawpuxinit: it would be better not to duplicate
168 !    these lines.
169      if(lpawu==0) then
170        hu(itypat)%f2_sla=zero
171      else if(lpawu==1) then
172        hu(itypat)%f2_sla=pawtab(itypat)%jpawu*5._dp
173      else if(lpawu==2) then
174        hu(itypat)%f2_sla=pawtab(itypat)%jpawu*14._dp/(One+hu(itypat)%f4of2_sla)
175      else if(lpawu==3) then
176        hu(itypat)%f2_sla=pawtab(itypat)%jpawu*6435._dp/(286._dp+&
177 &       195._dp*hu(itypat)%f4of2_sla+250._dp*hu(itypat)%f6of2_sla)
178      endif
179 
180      if(hu(itypat)%jpawu>tol4) hu(1)%jpawu_zero=.false.
181      ndim=2*lpawu+1
182 !     ndim1=2*hu(itypat)%lpawu+1
183      write(message,'(2a,i4)')  ch10,'  -------> For Correlated Species', itypat
184      call wrtout(std_out,  message,'COLL')
185 !     allocate(hu(itypat)%vee(ndim,ndim,ndim,ndim))
186      ABI_MALLOC(hu(itypat)%uqmc,(ndim*(2*ndim-1)))
187      ABI_MALLOC(hu(itypat)%udens,(2*ndim,2*ndim))
188      ABI_MALLOC(xij,(2*ndim,2*ndim))
189      if(t2g==0.and.x2my2d==0) then
190        hu(itypat)%vee => pawtab(itypat)%vee
191 !   t2g case begin
192      else if(t2g==1.and.hu(itypat)%lpawu==1) then
193        ABI_MALLOC(hu(itypat)%vee,(ndim,ndim,ndim,ndim))
194        n=0
195        do m=1,5
196          if((m/=1.and.m/=2.and.m/=4)) cycle
197          n=n+1
198          ns=0
199          do ms=1,5
200            if((ms/=1.and.ms/=2.and.ms/=4)) cycle
201            ns=ns+1
202            n1=0
203            do m1=1,5
204              if((m1/=1.and.m1/=2.and.m1/=4)) cycle
205              n1=n1+1
206              ns1=0
207              do ms1=1,5
208                if((ms1/=1.and.ms1/=2.and.ms1/=4)) cycle
209                ns1=ns1+1
210                hu(itypat)%vee(n,ns,n1,ns1)=pawtab(itypat)%vee(m,ms,m1,ms1)
211              enddo
212            enddo
213          enddo
214        enddo
215 !   t2g case end
216 !   x2my2d case begin
217      else if(x2my2d==1.and.hu(itypat)%lpawu==0) then
218        ABI_MALLOC(hu(itypat)%vee,(ndim,ndim,ndim,ndim))
219        hu(itypat)%vee(1,1,1,1)=pawtab(itypat)%upawu
220      endif
221 !   x2my2d case end
222 
223      hu(itypat)%udens=zero
224      ij=0
225      do ms=1,2*ndim-1
226          xij(ms,ms)=0
227        do ms1=ms+1,2*ndim
228          ij=ij+1
229          xij(ms,ms1)=ij
230          xij(ms1,ms)=ij
231          if(ms<=ndim.and.ms1>ndim) then
232            m1 = ms1 - ndim
233            m  = ms
234            hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)
235            hu(itypat)%udens(ms,ms1)= hu(itypat)%vee(m,m1,m,m1)
236            hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
237          else if(ms<=ndim.and.ms1<=ndim) then
238            m1 = ms1
239            m  = ms
240            hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)-hu(itypat)%vee(m,m1,m1,m)
241            hu(itypat)%udens(ms,ms1)= hu(itypat)%uqmc(ij)
242            hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
243          else
244            m1 = ms1 - ndim
245            m  = ms  - ndim
246            hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)-hu(itypat)%vee(m,m1,m1,m)
247            hu(itypat)%udens(ms,ms1)= hu(itypat)%uqmc(ij)
248            hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
249          endif
250        enddo
251      enddo
252      xij(2*ndim,2*ndim)=0
253      write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation "
254      call wrtout(std_out,  message,'COLL')
255      write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
256      call wrtout(std_out,  message,'COLL')
257 !     xtemp1b=0.d0
258 ! ====================================
259 !  Print hu(iatom)%uqmc
260 ! ====================================
261      ij1=-10
262      ij2=-10
263      ij=0
264      do i=1,2*ndim
265        do j=i+1,2*ndim
266          ij=ij+1
267          if(j==i+1) ij1=ij
268          if(j==2*ndim) ij2=ij
269        enddo
270 !       write(std_out,*) itypat
271 !       do m=1,i
272 !        write(std_out,*) i,m
273 !        write(std_out,*) xij(i,m)
274 !        write(std_out,*) ij1,ij2
275 !       enddo
276        if(i==1)               write(message,'(i3,14f7.3)') &
277 &                              i,xtemp, (hu(itypat)%uqmc(m),m=ij1,ij2)
278        if(i/=2*ndim.and.i/=1) write(message,'(i3,14f7.3)') i, &
279 &        (hu(itypat)%uqmc(xij(i,m)), m=1,i-1),xtemp, (hu(itypat)%uqmc(m),m=ij1,ij2)
280        if(i==2*ndim)          write(message,'(i3,14f7.3)') i, &
281 &                  (hu(itypat)%uqmc(xij(i,m)), m=1,i-1),xtemp
282        call wrtout(std_out,  message,'COLL')
283      enddo
284        write(message,'(5x,a)') "--------------------------------------------------------"
285        call wrtout(std_out,  message,'COLL')
286      ABI_FREE(xij)
287    else
288      hu(itypat)%upawu=zero
289      hu(itypat)%jpawu=zero
290 !     allocate(hu(itypat)%vee(0,0,0,0))
291    endif
292  enddo ! itypat
293 
294 end subroutine init_hu

m_hu/print_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 print_hu

FUNCTION

  print density density interaction (used for DFT+DMFT)

INPUTS

  ntypat = number of species
  prtopt = option for printing
  hu <type(hu_type)> = data for the interaction in DMFT.

OUTPUT

SOURCE

415 subroutine print_hu(hu,ntypat,prtopt)
416 
417  use defs_basis
418  use m_crystal, only : crystal_t
419  implicit none
420 
421 !Arguments ------------------------------------
422 !type
423  integer, intent(in):: ntypat
424  type(hu_type),intent(in) :: hu(ntypat)
425  integer :: prtopt
426 
427 !Local variables-------------------------------
428  integer :: itypat
429  integer :: lpawu,ms,ms1,m,ndim
430  character(len=500) :: message
431 ! *********************************************************************
432 
433  do itypat = 1 , ntypat
434    lpawu=hu(itypat)%lpawu
435    if(lpawu/=-1) then
436      ndim=2*lpawu+1
437      write(message,'(2a,i4)')  ch10,'  -------> For Correlated species'
438      call wrtout(std_out,  message,'COLL')
439      if(prtopt==0) then
440        write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in Slm basis "
441      else if(prtopt==1) then
442        write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in diagonal basis"
443      else if(prtopt==2) then
444        write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in Ylm basis"
445      else if(prtopt==3) then
446        write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation in JMJ basis"
447      endif
448      call wrtout(std_out,  message,'COLL')
449      write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
450      call wrtout(std_out,  message,'COLL')
451        do ms=1,2*ndim
452           write(message,'(i3,14f7.3)') &
453 &          ms, (hu(itypat)%udens(ms,ms1),ms1=1,2*ndim)
454           call wrtout(std_out,  message,'COLL')
455        enddo
456        write(message,'(5x,a)') "--------------------------------------------------------"
457        call wrtout(std_out,  message,'COLL')
458    endif ! lpawu/=1
459  enddo ! ntypat
460 
461 
462 end subroutine print_hu

m_hu/printvee_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 printvee_hu

FUNCTION

  print vee

INPUTS

  vee = number of species
  lpawu = value of l

OUTPUT

SOURCE

1046 subroutine printvee_hu(ndim,vee,prtopt,basis,upawu,f2)
1047 
1048  use defs_basis
1049  use m_crystal, only : crystal_t
1050  implicit none
1051 
1052 !Arguments ------------------------------------
1053 !type
1054  integer,intent(in) :: ndim,prtopt
1055  real(dp), intent(in) :: vee(ndim,ndim,ndim,ndim)
1056  real(dp), optional, intent(in) :: upawu,f2
1057  character(len=*), intent(in) :: basis
1058 
1059 !Local variables-------------------------------
1060  integer :: abcomp,m1,m2,mi,mj,mk,ml
1061  real(dp),allocatable :: b0(:,:)
1062  real(dp),allocatable :: a2pp(:,:)
1063  real(dp),allocatable :: b2pp(:,:)
1064  real(dp):: aver52,aver72,avernondiag,averall,averu,averj,averurestricted
1065  character(len=2000) :: message
1066 ! *********************************************************************
1067   write(message,'(5a)') ch10,&
1068 &  '  Coulomb interaction in the ', trim(basis),' basis'
1069   call wrtout(std_out,message,'COLL')
1070 
1071  if(prtopt>=2) then
1072 
1073    write(message,'(2a)')  ch10," <mi,mi|vee|mi mi> : U1"
1074    call wrtout(std_out,  message,'COLL')
1075    do mi=1,ndim
1076      write(message,'(4i4,3x,e10.3)')   mi,mi,mi,mi,vee(mi,mi,mi,mi)
1077      call wrtout(std_out,  message,'COLL')
1078    enddo
1079 
1080    write(message,'(2a)')  ch10," <mi,mj|vee|mi mj> : U2"
1081    call wrtout(std_out,  message,'COLL')
1082    do mi=1,ndim
1083      do mj=mi+1,ndim
1084        write(message,'(4i4,3x,e10.3)')   mi,mj,mi,mj,vee(mi,mj,mi,mj)
1085        call wrtout(std_out,  message,'COLL')
1086      enddo
1087    enddo
1088 
1089    write(message,'(2a)')  ch10," <mi,mj|vee|mj mi> : J"
1090    call wrtout(std_out,  message,'COLL')
1091    do mi=1,ndim
1092      do mj=mi+1,ndim
1093        write(message,'(4i4,3x,e10.3)')   mi,mj,mj,mi,vee(mi,mj,mj,mi)
1094        call wrtout(std_out,  message,'COLL')
1095      enddo
1096    enddo
1097 
1098    write(message,'(2a)')  ch10," <mi,mi|vee|mj mj> : J"
1099    call wrtout(std_out,  message,'COLL')
1100    do mi=1,ndim
1101      do mj=mi+1,ndim
1102        write(message,'(4i4,3x,e10.3)')   mi,mi,mj,mj,vee(mi,mi,mj,mj)
1103        call wrtout(std_out,  message,'COLL')
1104      enddo
1105    enddo
1106 
1107    write(message,'(2a)')  ch10," vee is non zero also for"
1108    call wrtout(std_out,  message,'COLL')
1109    do mi=1,ndim
1110      do mj=1,ndim
1111        do mk=1,ndim
1112          do ml=1,ndim
1113             if((.not.(mi==mk.and.mj==ml)).and.(.not.(mi==ml.and.mj==mk)).and.&
1114 &               .not.(mi==mj.and.mk==ml)) then
1115               if(vee(mi,mj,mk,ml)>tol8) then
1116                 write(message,'(4i4,3x,e10.3)')   mi,mj,mk,ml,vee(mi,mj,mk,ml)
1117                 call wrtout(std_out,  message,'COLL')
1118               endif
1119             endif
1120          enddo
1121        enddo
1122      enddo
1123    enddo
1124    write(message,'(a)')  ch10
1125    call wrtout(std_out,  message,'COLL')
1126 
1127  endif
1128 
1129  if(prtopt>=1) then
1130 
1131    write(message,'(2x,a,3x,14f10.4)') "Um1m2=Vee(m1,m2,m1,m2)"
1132    call wrtout(std_out,  message,'COLL')
1133    write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1134    call wrtout(std_out,  message,'COLL')
1135    do m1=1,ndim
1136      write(message,'(2x,i4,3x,14f10.4)') m1,(vee(m1,m2,m1,m2),m2=1,ndim)
1137      call wrtout(std_out,  message,'COLL')
1138    enddo
1139    write(message,'(a)') ch10;  call wrtout(std_out,  message,'COLL')
1140 
1141    write(message,'(2x,a,3x,14f10.4)') "Jm1m2=Vee(m1,m2,m2,m1)"
1142    call wrtout(std_out,  message,'COLL')
1143    write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1144    call wrtout(std_out,  message,'COLL')
1145    do m1=1,ndim
1146      write(message,'(2x,i4,3x,14f10.4)') m1,(vee(m1,m2,m2,m1),m2=1,ndim)
1147      call wrtout(std_out,  message,'COLL')
1148    enddo
1149    write(message,'(a)') ch10;  call wrtout(std_out,  message,'COLL')
1150 
1151    write(message,'(2x,a,3x,14f10.4)') "Udens(m1,m2)"
1152    call wrtout(std_out,  message,'COLL')
1153    write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1154    call wrtout(std_out,  message,'COLL')
1155    do m1=1,ndim
1156      write(message,'(2x,i4,3x,14f10.4)') m1,(vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1),m2=1,ndim)
1157      call wrtout(std_out,  message,'COLL')
1158    enddo
1159    write(message,'(a)') ch10;  call wrtout(std_out,  message,'COLL')
1160 
1161    if(prtopt>=3) then
1162      if(ndim==7) then
1163        averu=zero
1164        do m1=1,7
1165          do m2=1,7
1166            averu=vee(m1,m2,m1,m2)+averu
1167          enddo
1168        enddo
1169        averu=averu/49.d0
1170        averurestricted=zero
1171        do m1=1,7
1172          do m2=1,7
1173            if(m1/=m2) averurestricted=vee(m1,m2,m1,m2)+averurestricted
1174          enddo
1175        enddo
1176        averurestricted=averurestricted/42.d0
1177        averj=zero
1178        do m1=1,7
1179          do m2=1,7
1180            if(m1/=m2) averj=vee(m1,m2,m2,m1)+averj
1181          enddo
1182        enddo
1183        averj=averj/42
1184        averall=zero
1185        do m1=1,7
1186          do m2=1,7
1187            if(m1/=m2) averall=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+averall
1188          enddo
1189        enddo
1190        averall=averall/42
1191          !write(6,*) "averages Ylm U, U restricted,J, U-J",averu,averurestricted,averj,averall
1192 
1193      endif
1194 
1195      if(ndim==14.and.trim(basis)=='jmj') then
1196        aver52=zero
1197        do m1=1,6
1198          do m2=1,6
1199            aver52=vee(m1,m2,m1,m2)+aver52
1200          enddo
1201        enddo
1202        aver52=aver52/36.d0
1203        aver72=zero
1204        do m1=7,14
1205          do m2=7,14
1206             aver72=vee(m1,m2,m1,m2)+aver72
1207          enddo
1208        enddo
1209        aver72=aver72/64.d0
1210        avernondiag=zero
1211        do m1=1,6
1212          do m2=7,14
1213            avernondiag=vee(m1,m2,m1,m2)+avernondiag
1214          enddo
1215        enddo
1216        avernondiag=avernondiag/48.d0
1217        averall=zero
1218        do m1=1,14
1219          do m2=1,14
1220            averall=vee(m1,m2,m1,m2)+averall
1221          enddo
1222        enddo
1223        averall=averall/196.d0
1224          !write(6,*) "U averages",aver52,aver72,avernondiag,averall
1225 
1226 
1227 
1228        aver52=zero
1229        do m1=1,6
1230          do m2=1,6
1231            if(m1/=m2) aver52=vee(m1,m2,m2,m1)+aver52
1232          enddo
1233        enddo
1234        aver52=aver52/30.d0
1235        aver72=zero
1236        do m1=7,14
1237          do m2=7,14
1238            if(m1/=m2) aver72=vee(m1,m2,m2,m1)+aver72
1239          enddo
1240        enddo
1241        aver72=aver72/56.d0
1242        avernondiag=zero
1243        do m1=1,6
1244          do m2=7,14
1245            avernondiag=vee(m1,m2,m2,m1)+avernondiag
1246          enddo
1247        enddo
1248        avernondiag=avernondiag/48.d0
1249        averall=zero
1250        do m1=1,14
1251          do m2=1,14
1252            if(m1/=m2) averall=vee(m1,m2,m2,m1)+averall
1253          enddo
1254        enddo
1255        averall=averall/182.d0
1256          !write(6,*) "J averages",aver52,aver72,avernondiag,averall
1257 
1258 
1259 
1260 
1261        aver52=zero
1262        do m1=1,6
1263          do m2=1,6
1264            if(m1/=m2) aver52=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+aver52
1265          enddo
1266        enddo
1267        aver52=aver52/30.d0
1268        aver72=zero
1269        do m1=7,14
1270          do m2=7,14
1271            if(m1/=m2) aver72=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+aver72
1272          enddo
1273        enddo
1274        aver72=aver72/56.d0
1275        avernondiag=zero
1276        do m1=1,6
1277          do m2=7,14
1278            avernondiag=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+avernondiag
1279          enddo
1280        enddo
1281        avernondiag=avernondiag/48.d0
1282        averall=zero
1283        do m1=1,14
1284          do m2=1,14
1285            if(m1/=m2) averall=vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)+averall
1286          enddo
1287        enddo
1288        averall=averall/182.d0
1289        !write(6,*) "U-J averages",aver52,aver72,avernondiag,averall
1290 
1291      endif
1292    endif
1293 
1294    if (present(upawu)) then
1295      ABI_MALLOC(a2pp,(ndim,ndim))
1296      ABI_MALLOC(b2pp,(ndim,ndim))
1297      ABI_MALLOC(b0,(ndim,ndim))
1298 
1299 !     write(message,'(2x,a,3x,14f10.4)') "For check with respect to Slater's paper"
1300 !     call wrtout(std_out,  message,'COLL')
1301 !     ABI_MALLOC(f0,(ndim,ndim,ndim,ndim))
1302 !     write(message,'(2x,a,3x,14f10.4)') "Vee(m1,m2,m1,m2)-F0*ao(m1,m2)"
1303 !     call wrtout(std_out,  message,'COLL')
1304 !     write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1305 !     call wrtout(std_out,  message,'COLL')
1306 !
1307 !     do m1=1,ndim
1308 !       write(message,'(2x,i4,3x,14f10.4)') m1,(vee(m1,m2,m1,m2)-upawu,m2=1,ndim)
1309 !       call wrtout(std_out,  message,'COLL')
1310 !     enddo
1311 !
1312 !     f0=zero
1313 !     do m1=1,ndim
1314 !     f0(m1,m1,m1,m1)=upawu
1315 !     enddo
1316 !     write(message,'(a)') ch10
1317 !     call wrtout(std_out,  message,'COLL')
1318 !     write(message,'(2x,a,3x,14f10.4)') "Vee(m1,m2,m2,m1)-F0*b0(m1,m2)"
1319 !     call wrtout(std_out,  message,'COLL')
1320 !     write(message,'(2x,4x,14(2x,i8))') (m1,m1=1,ndim)
1321 !     call wrtout(std_out,  message,'COLL')
1322 !     do m1=1,ndim
1323 !       write(message,'(2x,i4,3x,14f10.4)') m1,(vee(m1,m2,m2,m1)-f0(m1,m2,m2,m1),m2=1,ndim)
1324 !       call wrtout(std_out,  message,'COLL')
1325 !     enddo
1326 !     ABI_FREE(f0)
1327 
1328 
1329      b0=zero
1330      do m1=1,ndim
1331      b0(m1,m1)=upawu
1332      enddo
1333      abcomp=0
1334      if(ndim==3.and.present(f2).and.(trim(basis)=='slm')) then
1335        a2pp(:,:) = RESHAPE((/1,-2,1,-2,4,-2,-1,-2,-1/),(/3,3/))
1336        a2pp=a2pp/25*f2+upawu
1337        b2pp(:,:) = RESHAPE((/1,3,6,3,4,3,6,3,1/),(/3,3/))
1338        b2pp=b2pp/25*f2+b0
1339        abcomp=1
1340      else if(ndim==6.and.present(f2).and.(trim(basis)=='jmj')) then
1341        ABI_MALLOC(a2pp,(6,6))
1342        ABI_MALLOC(b2pp,(6,6))
1343        a2pp(:,:)=RESHAPE((/0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,-1,-1,1,0,&
1344 &       0,-1,1,1,-1,0,0,-1,1,1,-1,0,0,1,-1,-1,1/),(/6,6/))
1345        a2pp=a2pp/25*f2+upawu
1346        b2pp(:,:)=RESHAPE((/0,0,1,2,3,4,0,0,4,3,2,1,1,4,1,2,2,0,2,3,&
1347 &       2,1,0,2,3,2,2,0,1,2,4,1,0,2,2,1/),(/6,6/))
1348        b2pp=b2pp/25*f2+b0
1349        abcomp=1
1350      endif
1351      if(mod(ndim,3)==0.and.present(f2).and.abcomp==1) then
1352        write(message,'(2x,a)') "Exact result for Umm is"
1353        call wrtout(std_out,message,'COLL')
1354        do m1=1,ndim
1355          write(message,'(2x,i4,3x,14f10.4)') m1,(a2pp(m1,m2),m2=1,ndim)
1356          call wrtout(std_out,  message,'COLL')
1357        enddo
1358        write(message,'(a)') ch10;  call wrtout(std_out,  message,'COLL')
1359        write(message,'(2x,a,3x,14f10.4)') "Exact result for Jmm is"
1360        call wrtout(std_out,message,'COLL')
1361        do m1=1,ndim
1362          write(message,'(2x,i4,3x,14f10.4)') m1,(b2pp(m1,m2),m2=1,ndim)
1363          call wrtout(std_out,message,'COLL')
1364        enddo
1365        write(message,'(a)') ch10;  call wrtout(std_out,  message,'COLL')
1366      endif
1367      ABI_FREE(a2pp)
1368      ABI_FREE(b2pp)
1369      ABI_FREE(b0)
1370    endif
1371 
1372  endif
1373 
1374 end subroutine printvee_hu

m_hu/reddd [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 reddd

FUNCTION

INPUTS

OUTPUT

SOURCE

1482 function reddd(mi,ndim)
1483 
1484  use defs_basis
1485  implicit none
1486 
1487 !Arguments ------------------------------------
1488 !scalars
1489  integer,intent(in) :: mi,ndim
1490  integer :: reddd
1491 ! *************************************************************************
1492 
1493  if(mi<ndim+1)  reddd=mi
1494  if(mi>=ndim+1) reddd=mi-ndim
1495 
1496 end function reddd

m_hu/rotatevee_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 rotatevee_hu

FUNCTION

  interaction udens in recomputed from new vee.

INPUTS

  ntypat = number of species
  cryst_struc <type(crystal_t)>=crystal structure data

OUTPUT

 SIDE EFFECT
  hu <type(hu_type)> = data for the interaction in DMFT.

SOURCE

 567 subroutine rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,rot_mat,udens_atoms,rot_type)
 568 
 569  use defs_basis
 570  use m_errors
 571 
 572  use m_crystal,          only : crystal_t
 573  implicit none
 574 
 575 !Arguments ------------------------------------
 576 !type
 577  type(crystal_t),intent(in) :: cryst_struc
 578  integer, intent(in):: nsppol,nspinor,pawprtvol,rot_type
 579  type(hu_type),intent(in) :: hu(cryst_struc%ntypat)
 580  type(coeff2c_type),optional,intent(in) :: rot_mat(cryst_struc%natom,nsppol)
 581  type(coeff2_type),intent(inout) :: udens_atoms(cryst_struc%natom)
 582 
 583 !Local variables-------------------------------
 584  integer :: dim_vee,iatom,itypat
 585  integer :: lpawu,m1,m2,m3,m4,mi,mj,mk,ml,natom,ndim,tndim
 586  integer :: ms,ms1,nat_correl
 587  character(len=500) :: message
 588  character(len=30) :: basis_vee
 589  real(dp) :: xsum,xsum2,xsumnew,xsum2new,uaver
 590  complex(dpc),allocatable :: temp_mat(:,:)
 591  complex(dpc),allocatable :: temp_mat2(:,:)
 592  real(dp),allocatable :: veetemp(:,:,:,:),fk(:)
 593  complex(dpc),allocatable :: veeylm(:,:,:,:),veeslm(:,:,:,:)
 594  complex(dpc),allocatable :: veetmp(:,:,:),veetmp_slm(:,:,:),veetmp_ylm(:,:,:)
 595  complex(dpc),allocatable :: veejmj(:,:,:,:),veerotated(:,:,:,:),veenew(:,:,:,:)
 596  complex(dpc),allocatable :: veeylm2(:,:,:,:),veeslm2(:,:,:,:)
 597 ! *********************************************************************
 598  natom=cryst_struc%natom
 599 
 600 !================================================
 601 !  NSPINOR = 2 and J=0
 602 !================================================
 603  if(hu(1)%jpawu_zero.and.nspinor==2) then
 604 ! if(3==4) then
 605 
 606 !   call vee2udens_hu(hu,cryst_struc%ntypat,2)
 607    do iatom=1,natom
 608      itypat=cryst_struc%typat(iatom)
 609      lpawu=hu(itypat)%lpawu
 610      if(lpawu.ne.-1) then
 611        ndim=2*lpawu+1
 612        write(message,'(2a,i4)')  ch10,'  -------> For Correlated Species', itypat
 613        call wrtout(std_out,  message,'COLL')
 614     !   write(std_out,*)"ndim",ndim
 615     !   write(std_out,'(a,20e10.4)')"udens before vee2udensaomt", udens_atoms(1)%value
 616 !       write(std_out,*)"vee",hu(cryst_struc%typat(iatom))%vee
 617        call vee2udensatom_hu(ndim,nspinor,udens_atoms(iatom)%value,hu(cryst_struc%typat(iatom))%vee,"Slm")
 618     !   write(std_out,*)"udens after vee2udensatom", udens_atoms(1)%value
 619      endif
 620    enddo
 621   ! write(std_out,*)"udensafter after", udens_atoms(1)%value
 622 
 623 
 624 !================================================
 625 !  NSPINOR = 2 and J/=0
 626 !================================================
 627  else if (nspinor==2) then
 628    write(std_out,*)
 629    write(std_out,*)" == Rotate interaction in the level basis "
 630 
 631 
 632    !rot_type=0 ! keep original Slm basis
 633    !rot_type=2 ! rotation to the Ylm basis ! for tests
 634    !rot_type=4 ! use a rotation from Ylm basis to a new basis
 635    !rot_type=3 ! rotation to the JmJ Basis ! for tests
 636    !rot_type=1 ! use a rotation for diago of dmat,green, levels..
 637 
 638    do iatom=1,natom
 639      itypat=cryst_struc%typat(iatom)
 640      lpawu=hu(itypat)%lpawu
 641      if(lpawu.ne.-1) then
 642        ndim=2*lpawu+1
 643        if(pawprtvol>=3) then
 644 !         write(message,'(2a)')  ch10," VEE INPUT AVANT TRANSFORMATION"
 645 !         call wrtout(std_out,  message,'COLL')
 646 !         call printvee_hu(ndim,hu(itypat)%vee,1,'Slm')
 647        endif
 648        write(message,'(2a,i4)')  ch10,'  -------> For Correlated atom', iatom
 649        call wrtout(std_out,  message,'COLL')
 650        tndim=nspinor*ndim
 651        ABI_MALLOC(veejmj,(tndim,tndim,tndim,tndim))
 652        ABI_MALLOC(veeslm,(ndim,ndim,ndim,ndim))
 653        ABI_MALLOC(veetmp_slm,(ndim,ndim,4))
 654        ABI_MALLOC(veetmp_ylm,(ndim,ndim,4))
 655        ABI_MALLOC(veetmp,(ndim,ndim,4))
 656        ABI_MALLOC(veeslm2,(tndim,tndim,tndim,tndim))
 657        ABI_MALLOC(veeylm,(ndim,ndim,ndim,ndim))
 658        ABI_MALLOC(veeylm2,(tndim,tndim,tndim,tndim))
 659        ABI_MALLOC(veerotated,(tndim,tndim,tndim,tndim))
 660        ABI_MALLOC(fk,(0:lpawu))
 661        fk(0)=hu(itypat)%upawu
 662        if (lpawu==1) then
 663          fk(1)=hu(itypat)%jpawu*5
 664        else if (lpawu==2) then
 665          fk(1)=hu(itypat)%jpawu*14._dp/(One+hu(itypat)%f4of2_sla)
 666          fk(2)=fk(1)*hu(itypat)%f4of2_sla
 667        else if (lpawu==3) then
 668          fk(1)=hu(itypat)%jpawu*6435._dp/(286._dp+195._dp*hu(itypat)%f4of2_sla+250._dp*hu(itypat)%f6of2_sla)
 669          fk(2)=fk(1)*hu(itypat)%f4of2_sla
 670          fk(3)=fk(1)*hu(itypat)%f6of2_sla
 671        endif
 672 
 673 
 674 
 675 
 676 !      ==================================
 677 !      First define veeslm
 678 !      ==================================
 679 
 680 
 681 !!     veeslm2(s1m1,s2m2,s3m3,s4m4)=vee(m1,m2,m3,m4)*delta_s1s3*delta_s2s4
 682        call vee_ndim2tndim_hu(lpawu,hu(itypat)%vee,veeslm2,1)
 683 
 684 !!     veeslm(m1,m2,m3,m4)=cmplx(vee(m1,m2,m3,m4),zero)
 685      !  ABI_FREE(veeslm)
 686        veeslm(:,:,:,:)=cmplx(hu(itypat)%vee(:,:,:,:),zero)
 687        !ABI_FREE(veeslm)! ici cela plante
 688 
 689 !!     build udens in the Slm basis and print it
 690        call vee2udensatom_hu(ndim,nspinor,udens_atoms(iatom)%value,real(veeslm),"slm")
 691        !ABI_FREE(veeslm) ! ici cela plante
 692 
 693        dim_vee=ndim
 694        basis_vee='Slm'
 695 !     first print veeslm
 696        !call printvee_hu(ndim,real(veeslm),1,basis_vee)
 697        call printvee_hu(tndim,real(veeslm2),1,basis_vee)
 698 
 699       ! ABI_FREE(veeslm)
 700 
 701 !      ==================================
 702 !      Then compute veerotated
 703 !      ==================================
 704 
 705 !      In the basis where levels/densitymatrix/green function is  diagonalized
 706 !      ================================================================================
 707        if (rot_type==1) then
 708 !      ---------------------
 709 
 710          veerotated=czero
 711          do m1=1,tndim
 712            do m2=1,tndim
 713              do m3=1,tndim
 714                do m4=1,tndim
 715                  do mi=1,tndim
 716                    do mj=1,tndim
 717                      do mk=1,tndim
 718                        do ml=1,tndim
 719                           veerotated(m1,m2,m3,m4)= veerotated(m1,m2,m3,m4) + &
 720 &                          conjg(rot_mat(iatom,1)%value(mi,m1))* &
 721 &                          conjg(rot_mat(iatom,1)%value(mj,m2))* &
 722 &                                rot_mat(iatom,1)%value(mk,m3)* &
 723 &                                rot_mat(iatom,1)%value(ml,m4)* &
 724 &                              veeslm2(mi,mj,mk,ml)
 725                        enddo
 726                      enddo
 727                    enddo
 728                  enddo
 729                enddo
 730              enddo
 731            enddo
 732          enddo
 733 
 734          dim_vee=tndim
 735          basis_vee='Diagonal basis from Slm basis'
 736 !      In the Ylm basis
 737 !      ================================================================================
 738        else if (rot_type==2.or.rot_type==3.or.rot_type==4) then
 739 !      ---------------------------
 740 
 741 
 742 !      Change basis from slm to ylm basis
 743          call vee_slm2ylm_hu(lpawu,veeslm,veeylm,1,2)
 744 
 745          basis_vee='Ylm'
 746          dim_vee=ndim
 747 
 748 !        print interaction matrix in the ylm basis
 749          call printvee_hu(ndim,real(veeylm),1,basis_vee,hu(itypat)%upawu)
 750 
 751 !        print interaction matrix in the ylm basis from Slater tables
 752          if(pawprtvol>=3) call udens_slatercondon_hu(fk,lpawu)
 753 
 754          if (rot_type==3.or.rot_type==4) then
 755 !        ---------------------------
 756 !
 757 
 758 !          built large matrix
 759            call vee_ndim2tndim_hu(lpawu,real(veeylm),veeylm2,1)
 760 
 761 
 762            call printvee_hu(tndim,real(veeylm2),1,basis_vee)
 763 
 764          endif
 765 
 766 
 767 !      In the JmJ basis
 768 !      ================================================================================
 769          if (rot_type==3) then
 770 
 771 !          apply change of basis
 772            call vee_ylm2jmj_hu(lpawu,veeylm2,veejmj,1)
 773 
 774 !        print interaction matrix in the JMJ basis from Inglis and Julien tables
 775            if(pawprtvol>=3) call udens_inglis_hu(fk,lpawu)
 776 
 777 !          new dimension
 778 
 779 
 780 
 781            dim_vee=tndim
 782            basis_vee='JmJ'
 783 
 784          else if (rot_type==4) then
 785 
 786            call udens_inglis_hu(fk,lpawu)
 787            veerotated=czero
 788 
 789            do m1=1,tndim
 790              do m2=1,tndim
 791                do m3=1,tndim
 792                  do m4=1,tndim
 793                    do mi=1,tndim
 794                      do mj=1,tndim
 795                        do mk=1,tndim
 796                          do ml=1,tndim
 797                             veerotated(m1,m2,m3,m4)= veerotated(m1,m2,m3,m4) + &
 798 &                          conjg(rot_mat(iatom,1)%value(mi,m1))* &
 799 &                          conjg(rot_mat(iatom,1)%value(mj,m2))* &
 800 &                                rot_mat(iatom,1)%value(mk,m3)* &
 801 &                                rot_mat(iatom,1)%value(ml,m4)* &
 802 &                                veeylm2(mi,mj,mk,ml)
 803                          enddo
 804                        enddo
 805                      enddo
 806                    enddo
 807                  enddo
 808                enddo
 809              enddo
 810            enddo
 811 
 812            dim_vee=tndim
 813            basis_vee='Diagonal basis from Ylm'
 814 
 815          endif
 816        endif
 817        ABI_MALLOC(veenew,(dim_vee,dim_vee,dim_vee,dim_vee))
 818        if(rot_type==0) then   ; veenew=veeslm
 819        else if(rot_type==1) then  ; veenew=veerotated
 820        else if(rot_type==2) then  ; veenew=veeylm
 821        else if(rot_type==3) then  ; veenew=veejmj
 822        else if(rot_type==4) then  ; veenew=veerotated
 823        endif
 824 
 825        call printvee_hu(dim_vee,real(veenew),1,basis_vee,hu(itypat)%upawu,hu(itypat)%f2_sla)
 826 !       call printvee_hu(dim_vee,real(veeylm),1,hu(itypat)%upawu)
 827 
 828        uaver=zero
 829        do ms=1,2*ndim
 830          do ms1=1,2*ndim
 831            udens_atoms(iatom)%value(ms,ms1)=veenew(ms,ms1,ms,ms1)-veenew(ms,ms1,ms1,ms)
 832            uaver=uaver+udens_atoms(iatom)%value(ms,ms1)
 833          enddo
 834        enddo
 835 
 836        call vee2udensatom_hu(ndim,nspinor,udens_atoms(iatom)%value,real(veenew),"basis_vee",prtonly=1)
 837 
 838 
 839        ABI_FREE(veenew)
 840        ABI_FREE(veeslm)
 841        ABI_FREE(veeslm2)
 842        ABI_FREE(veeylm)
 843        ABI_FREE(veeylm2)
 844        ABI_FREE(veejmj)
 845        ABI_FREE(veerotated)
 846        ABI_FREE(veetmp_slm)
 847        ABI_FREE(veetmp)
 848        ABI_FREE(veetmp_ylm)
 849        ABI_FREE(fk)
 850      endif
 851    enddo
 852    !ABI_ERROR("Aborting now!")
 853 
 854 !================================================
 855 !  NSPINOR = 1
 856 !================================================
 857  else if (nspinor==1) then
 858 
 859    nat_correl=0
 860    do iatom=1,natom
 861      itypat=cryst_struc%typat(iatom)
 862      lpawu=hu(itypat)%lpawu
 863      if(lpawu.ne.-1) then
 864        nat_correl=nat_correl+1
 865        if(nat_correl>1.and.(hu(itypat)%jpawu>tol4)) then
 866           write(message,'(3a)')  ch10,'  -------> Warning: several atoms: '&
 867 &         ,' not extensively tested '
 868           ABI_WARNING(message)
 869        endif
 870 
 871 !  ! ================================================================
 872 !  !  If rotation for spin 2 and rotation for spin 1 are not equal print
 873 !  !  then print a warning
 874 !  !  useful only for magnetic case
 875 !  ! ================================================================
 876        ndim=2*lpawu+1
 877        tndim=nspinor*ndim
 878        do m1=1,tndim
 879          do m2=1,tndim
 880            if(nsppol==2) then
 881              if(abs(rot_mat(iatom,1)%value(m1,m2)-rot_mat(iatom,2)%value(m1,m2))>tol4.and.pawprtvol>=3) then
 882                write(message,'(2a,i4)')  ch10,' rot_mat differs for value of isppol but value for isppol=2 not used'
 883                call wrtout(std_out,  message,'COLL')
 884                write(message,'(a,4e16.8)')  ch10,rot_mat(iatom,1)%value(m1,m2),rot_mat(iatom,2)%value(m1,m2)
 885                call wrtout(std_out,  message,'COLL', do_flush=.True.)
 886              endif
 887            endif
 888          end do
 889        end do
 890        ABI_MALLOC(temp_mat,(ndim,ndim))
 891        ABI_MALLOC(temp_mat2,(ndim,ndim))
 892        ABI_MALLOC(veetemp,(ndim,ndim,ndim,ndim))
 893        temp_mat(:,:)=czero
 894        temp_mat2(:,:)=czero
 895 
 896 !  ! =================================================
 897 !  !    See if rotation is complex or real
 898 !  ! =================================================
 899        do mi=1,ndim
 900          do m1=1,ndim
 901          if(abs(aimag(rot_mat(iatom,1)%value(mi,m1)))>tol8.and.pawprtvol>=3) then
 902               write(message,'(2a,2i6,2e14.3)')  ch10,"rot_mat is complex for", &
 903 &              mi,m1,rot_mat(iatom,1)%value(mi,m1)
 904               call wrtout(std_out,  message,'COLL')
 905            endif
 906          enddo
 907        enddo
 908 
 909 
 910 
 911 !  !    write vee for information with a classification.
 912        if(pawprtvol>=3) then
 913          call printvee_hu(ndim,hu(itypat)%vee,2,'original')
 914        endif
 915 
 916 !  !    Compute rotated vee.
 917        veetemp=zero
 918        do m1=1,ndim
 919          do m2=1,ndim
 920            do m3=1,ndim
 921              do m4=1,ndim
 922                do mi=1,ndim
 923                  do mj=1,ndim
 924                    do mk=1,ndim
 925                      do ml=1,ndim
 926 !                        if((mi==mk.and.mj==ml).or.(mi==ml.and.mj==mk)) then
 927                         veetemp(m1,m2,m3,m4)= veetemp(m1,m2,m3,m4) + &
 928 &                      real(   &
 929 &                          conjg(rot_mat(iatom,1)%value(mi,m1))* &
 930 &                          conjg(rot_mat(iatom,1)%value(mj,m2))* &
 931 &                                rot_mat(iatom,1)%value(mk,m3)* &
 932 &                                rot_mat(iatom,1)%value(ml,m4)* &
 933 &                            hu(itypat)%vee(mi,mj,mk,ml)&
 934                            )
 935 !                        endif
 936                      enddo
 937                    enddo
 938                  enddo
 939                enddo
 940              enddo
 941            enddo
 942          enddo
 943        enddo
 944        ABI_FREE(temp_mat)
 945        ABI_FREE(temp_mat2)
 946        xsum=zero
 947        xsum2=zero
 948        xsumnew=zero
 949        xsum2new=zero
 950        do m1=1,ndim
 951          do m2=1,ndim
 952            xsum=xsum+hu(itypat)%vee(m1,m2,m1,m2)
 953            xsum2=xsum2+hu(itypat)%vee(m1,m2,m2,m1)
 954            xsumnew=xsumnew+veetemp(m1,m2,m1,m2)
 955            xsum2new=xsum2new+veetemp(m1,m2,m2,m1)
 956          enddo
 957        enddo
 958        if(abs(xsum-xsumnew)>tol5.or.abs(xsum2-xsum2new)>tol5) then
 959          write(message,'(2a)')  ch10," BUG: New interaction after rotation do not respect sum rules"
 960          call wrtout(std_out,  message,'COLL')
 961          write(message,'(2a,2f14.3)')  ch10,' Comparison of \sum_{m1,m3} vee(m1,m3,m1,m3) before and after rotation is',&
 962 &         xsum,xsumnew
 963          call wrtout(std_out,  message,'COLL')
 964          write(message,'(2a,2f14.3)')  ch10,' Comparison of \sum_{m1,m3} vee(m1,m3,m3,m1) before and after rotation is',&
 965 &         xsum2,xsum2new
 966          call wrtout(std_out,  message,'COLL')
 967        endif
 968        if(pawprtvol>=3) then
 969          write(message,'(2a)')  ch10," VEE ROTATED"
 970          call wrtout(std_out,  message,'COLL')
 971          call printvee_hu(tndim,veetemp,2,'diag1')
 972          write(message,'(a)') ch10
 973          call wrtout(std_out,  message,'COLL')
 974        endif
 975 
 976        write(message,'(2a,i4)')  ch10,'  -------> For Correlated Species', itypat
 977        call wrtout(std_out,  message,'COLL')
 978 
 979        call vee2udensatom_hu(ndim,nspinor,udens_atoms(iatom)%value,veetemp,"diag")
 980 
 981 !       udens_atoms(iatom)%value=zero
 982 !       ij=0
 983 !       do ms=1,2*ndim-1
 984 !         do ms1=ms+1,2*ndim
 985 !           ij=ij+1
 986 !           if(ms<=ndim.and.ms1>ndim) then
 987 !             m1 = ms1 - ndim
 988 !             m  = ms
 989 !             hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)
 990 !             udens_atoms(iatom)%value(ms,ms1)= veetemp(m,m1,m,m1)
 991 !             udens_atoms(iatom)%value(ms1,ms)= udens_atoms(iatom)%value(ms,ms1)
 992 !           else if(ms<=ndim.and.ms1<=ndim) then
 993 !             m1 = ms1
 994 !             m  = ms
 995 !             hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
 996 !             udens_atoms(iatom)%value(ms,ms1)= hu(itypat)%uqmc(ij)
 997 !             udens_atoms(iatom)%value(ms1,ms)= udens_atoms(iatom)%value(ms,ms1)
 998 !           else
 999 !             m1 = ms1 - ndim
1000 !             m  = ms  - ndim
1001 !             hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1002 !             udens_atoms(iatom)%value(ms,ms1)= hu(itypat)%uqmc(ij)
1003 !             udens_atoms(iatom)%value(ms1,ms)= udens_atoms(iatom)%value(ms,ms1)
1004 !           endif
1005 !         enddo
1006 !       enddo
1007 !       write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation "
1008 !       call wrtout(std_out,  message,'COLL')
1009 !       write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
1010 !       call wrtout(std_out,  message,'COLL')
1011 !       do ms=1,2*ndim
1012 !          write(message,'(i3,14f7.3)') &
1013 !  &        ms, (udens_atoms(iatom)%value(ms,ms1),ms1=1,2*ndim)
1014 !          call wrtout(std_out,  message,'COLL')
1015 !       enddo
1016 !       write(message,'(5x,a)') "--------------------------------------------------------"
1017 !       call wrtout(std_out,  message,'COLL')
1018        ABI_FREE(veetemp)
1019      endif ! lpawu/=1
1020 !   call print_hu(hu,cryst_struc%ntypat,1)
1021 
1022    enddo ! iatom
1023 !   call print_hu(hu,cryst_struc%ntypat,1)
1024 !   call vee2udens_hu(hu,cryst_struc%ntypat,2)
1025  endif ! nspinor==1
1026 
1027 
1028 end subroutine rotatevee_hu

m_hu/vee2udens_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 print_hu

FUNCTION

  interaction udens in recomputed from new vee.

INPUTS

  ntypat = number of species
  prtopt = option for printing

OUTPUT

 SIDE EFFECT
  hu <type(hu_type)> = data for the interaction in DMFT.

SOURCE

483 subroutine vee2udens_hu(hu,ntypat,prtopt)
484 
485  use defs_basis
486  use m_crystal, only : crystal_t
487  implicit none
488 
489 !Arguments ------------------------------------
490 !type
491  integer, intent(in):: ntypat
492  type(hu_type),intent(inout) :: hu(ntypat)
493  integer :: prtopt
494 
495 !Local variables-------------------------------
496  integer :: ij,itypat
497  integer :: lpawu,m1,ms,ms1,m,ndim
498  character(len=500) :: message
499 ! *********************************************************************
500  do itypat=1,ntypat
501    lpawu=hu(itypat)%lpawu
502    if(lpawu.ne.-1) then
503      ndim=2*lpawu+1
504      write(message,'(2a,i4)')  ch10,'  -------> For Correlated Species', itypat
505      call wrtout(std_out,  message,'COLL')
506 
507      hu(itypat)%udens=zero
508      ij=0
509      do ms=1,2*ndim-1
510 !         xij(ms,ms)=0
511        do ms1=ms+1,2*ndim
512          ij=ij+1
513 !         xij(ms,ms1)=ij
514 !         xij(ms1,ms)=ij
515          if(ms<=ndim.and.ms1>ndim) then
516            m1 = ms1 - ndim
517            m  = ms
518            hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)
519            hu(itypat)%udens(ms,ms1)= hu(itypat)%vee(m,m1,m,m1)
520            hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
521          else if(ms<=ndim.and.ms1<=ndim) then
522            m1 = ms1
523            m  = ms
524            hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)-hu(itypat)%vee(m,m1,m1,m)
525            hu(itypat)%udens(ms,ms1)= hu(itypat)%uqmc(ij)
526            hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
527          else
528            m1 = ms1 - ndim
529            m  = ms  - ndim
530            hu(itypat)%uqmc(ij)=hu(itypat)%vee(m,m1,m,m1)-hu(itypat)%vee(m,m1,m1,m)
531            hu(itypat)%udens(ms,ms1)= hu(itypat)%uqmc(ij)
532            hu(itypat)%udens(ms1,ms)= hu(itypat)%udens(ms,ms1)
533          endif
534        enddo
535      enddo
536 !     xij(2*ndim,2*ndim)=0
537 !     write(message,'(a,5x,a)') ch10,"-------- Interactions in the density matrix representation "
538 !     call wrtout(std_out,  message,'COLL')
539 !     write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
540 !     call wrtout(std_out,  message,'COLL')
541    endif
542  enddo ! itypat
543  call print_hu(hu,ntypat,prtopt)
544 
545 
546 end subroutine vee2udens_hu

m_hu/vee2udensatom_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 vee2udensatom_hu

FUNCTION

  compute density density interaction (used for DFT+DMFT)

INPUTS

  ntypat = number of species
  hu <type(hu_type)> = data for the interaction in DMFT.

OUTPUT

SOURCE

1392 subroutine vee2udensatom_hu(ndim,nspinor,udens_atoms,veetemp,basis,prtonly)
1393 
1394  use defs_basis
1395  use m_crystal, only : crystal_t
1396  implicit none
1397 
1398 !Arguments ------------------------------------
1399 !type
1400  integer, intent(in):: ndim,nspinor
1401  real(dp),intent(out) :: udens_atoms(2*ndim,2*ndim)
1402  !real(dp), intent(in) :: veetemp(nspinor*ndim,nspinor*ndim,nspinor*ndim,nspinor*ndim)
1403  real(dp), intent(in) :: veetemp(ndim,ndim,ndim,ndim)
1404  character(len=*), intent(in) :: basis
1405  integer, intent(in), optional:: prtonly
1406 
1407 !Local variables-------------------------------
1408  integer :: ij,ms,ms1,m,m1,prt
1409  character(len=1000) :: message
1410 ! *********************************************************************
1411  prt=1
1412  if(present(prtonly)) prt=2
1413  if (nspinor==1.or.nspinor==2.and.prt<=1) then
1414    udens_atoms=zero
1415    ij=0
1416    do ms=1,2*ndim-1
1417      do ms1=ms+1,2*ndim
1418        ij=ij+1
1419        if(ms<=ndim.and.ms1>ndim) then
1420          m1 = ms1 - ndim
1421          m  = ms
1422 !         hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)
1423          udens_atoms(ms,ms1)= veetemp(m,m1,m,m1)
1424          udens_atoms(ms1,ms)= udens_atoms(ms,ms1)
1425         ! write(6,*)"A", ms,ms1,udens_atoms(ms,ms1)
1426        else if(ms<=ndim.and.ms1<=ndim) then
1427          m1 = ms1
1428          m  = ms
1429 !         hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1430          udens_atoms(ms,ms1)= veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1431          udens_atoms(ms1,ms)= udens_atoms(ms,ms1)
1432         ! write(6,*)"B", ms,ms1,udens_atoms(ms,ms1)
1433        else
1434          m1 = ms1 - ndim
1435          m  = ms  - ndim
1436 !         hu(itypat)%uqmc(ij)=veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1437          udens_atoms(ms,ms1)= veetemp(m,m1,m,m1)-veetemp(m,m1,m1,m)
1438          udens_atoms(ms1,ms)= udens_atoms(ms,ms1)
1439         ! write(6,*)"C", ms,ms1,udens_atoms(ms,ms1)
1440        endif
1441      enddo
1442    enddo
1443 
1444 ! else if(nspinor==2) then
1445 !
1446 !   do ms=1,2*ndim
1447 !     do ms1=1,2*ndim
1448 !       udens_atoms(ms,ms1)=veetemp(ms,ms1,ms,ms1)-veetemp(ms,ms1,ms1,ms)
1449 !     enddo
1450 !   enddo
1451 
1452  endif
1453 
1454 
1455  write(message,'(4a)') ch10,"-------- Interactions in the ",trim(basis)," basis "
1456  call wrtout(std_out,  message,'COLL')
1457  write(message,'(1x,14(2x,i5))') (m,m=1,2*ndim)
1458  call wrtout(std_out,  message,'COLL')
1459  do ms=1,2*ndim
1460     write(message,'(i3,14f7.3)') &
1461 &    ms, (udens_atoms(ms,ms1),ms1=1,2*ndim)
1462     call wrtout(std_out,  message,'COLL')
1463  enddo
1464  write(message,'(5x,a)') "--------------------------------------------------------"
1465  call wrtout(std_out,  message,'COLL')
1466 
1467 end subroutine vee2udensatom_hu

m_hu/vee_ndim2tndim_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 vee_ndim2tndim_hu

FUNCTION

 Change a matrix  of interaction of dimension [(2*lcor+1)]**2
 into a full spin and orbital interaction matrix of dimension [2*(2l+1)]**4

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (BA)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  lcor= angular momentum, size of the matrix is 2(2*lcor+1)
  mat_inp_c= real input matrix
  prtvol=printing volume
  option= 1 : Vout_(s1m1,s2m2,s3m3,s4m4)=Vinp_(m1,m2,m3,m4)*delta_s1s3*delta_s2s4

OUTPUT

  mat_out_c= Complex output matrix

NOTES

SOURCE

1752 subroutine vee_ndim2tndim_hu(lcor,mat_inp_c,mat_out_c,option)
1753 
1754  use defs_basis
1755  use m_errors
1756  use m_abicore
1757  implicit none
1758 
1759 !Arguments ---------------------------------------------
1760 !scalars
1761  integer,intent(in) :: lcor,option
1762 !arrays
1763  real(dp), intent(in) :: mat_inp_c(2*lcor+1,2*lcor+1,2*lcor+1,2*lcor+1)
1764  complex(dpc), intent(out) :: mat_out_c(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
1765 
1766 !Local variables ---------------------------------------
1767 !scalars
1768  integer :: m1,m2,m3,m4,is1,is2,is3,is4,ndim,s1,s2,s3,s4
1769 
1770 ! *********************************************************************
1771   ndim=2*lcor+1
1772   mat_out_c=czero
1773   if(option==1) then
1774     do m1=1,ndim
1775       do m2=1,ndim
1776         do m3=1,ndim
1777           do m4=1,ndim
1778             do is1=1,2
1779               do is2=1,2
1780 
1781                 is3=is1 ; is4=is2
1782 
1783                 s1=(is1-1)*ndim ; s2=(is2-1)*ndim ; s3=(is3-1)*ndim ; s4=(is4-1)*ndim
1784 
1785                 mat_out_c(m1+s1,m2+s2,m3+s3,m4+s4)=  cmplx(mat_inp_c(m1,m2,m3,m4),zero)
1786 
1787               enddo
1788             enddo
1789           enddo
1790         enddo
1791       enddo
1792     enddo
1793   endif
1794 
1795 
1796 end subroutine vee_ndim2tndim_hu

m_hu/vee_ndim2tndim_hu_r [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 vee_ndim2tndim_hu_r

FUNCTION

 Change a matrix  of interaction of dimension [(2*lcor+1)]**2
 into a full spin and orbital interaction matrix of dimension [2*(2l+1)]**4

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (BA)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  lcor= angular momentum, size of the matrix is 2(2*lcor+1)
  mat_inp_c= real input matrix
  prtvol=printing volume
  option= 1 : Vout_(s1m1,s2m2,s3m3,s4m4)=Vinp_(m1,m2,m3,m4)*delta_s1s3*delta_s2s4

OUTPUT

  mat_out_c= real output matrix

NOTES

SOURCE

1675 subroutine vee_ndim2tndim_hu_r(lcor,mat_inp_c,mat_out_c,option)
1676 
1677  use defs_basis
1678  use m_errors
1679  use m_abicore
1680  implicit none
1681 
1682 !Arguments ---------------------------------------------
1683 !scalars
1684  integer,intent(in) :: lcor,option
1685 !arrays
1686  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(:,:,:,:)
1687  real(dp), intent(out)       :: mat_out_c(:,:,:,:) !(2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1),2*(2*lcor+1))
1688 
1689 !Local variables ---------------------------------------
1690 !scalars
1691  integer :: m1,m2,m3,m4,is1,is2,is3,is4,ndim,s1,s2,s3,s4
1692 
1693 ! *********************************************************************
1694   ndim=2*lcor+1
1695   mat_out_c=czero
1696 
1697   if(option==1) then
1698     do m1=1,ndim
1699       do m2=1,ndim
1700         do m3=1,ndim
1701           do m4=1,ndim
1702             do is1=1,2
1703               do is2=1,2
1704 
1705                 is3=is1 ; is4=is2
1706 
1707                 s1=(is1-1)*ndim ; s2=(is2-1)*ndim ; s3=(is3-1)*ndim ; s4=(is4-1)*ndim
1708 
1709                 mat_out_c(m1+s1,m2+s2,m3+s3,m4+s4)=  mat_inp_c(m1,m2,m3,m4)
1710 
1711               enddo
1712             enddo
1713           enddo
1714         enddo
1715       enddo
1716     enddo
1717   endif
1718 
1719 
1720 end subroutine vee_ndim2tndim_hu_r

m_hu/vee_slm2ylm_hu [ Functions ]

[ Top ] [ m_hu ] [ Functions ]

NAME

 vee_slm2ylm_hu

FUNCTION

 For a given angular momentum lcor, change a matrix  of interaction of dimension (2*lcor+1)
 from the Slm to the Ylm basis if option==1 or from Ylm to Slm if !option==2

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (BA)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  lcor= angular momentum, size of the matrix is 2(2*lcor+1)
  mat_inp_c= Input matrix
  option= 1  Change matrix from Slm to Ylm basis
          2  Change matrix from Ylm to Slm basis
  optspin=  1  Spin up are first
            2  Spin dn are first
  prtvol=printing volume

OUTPUT

  mat_inp_c= Output matrix in Ylm or Slm basis according to option

NOTES

SOURCE

1530 subroutine vee_slm2ylm_hu(lcor,mat_inp_c,mat_out_c,option,prtvol)
1531 
1532  use defs_basis
1533  use m_errors
1534  use m_abicore
1535  implicit none
1536 
1537 !Arguments ---------------------------------------------
1538 !scalars
1539  integer,intent(in) :: lcor,option,prtvol
1540 !arrays
1541  complex(dpc), intent(in) :: mat_inp_c(2*lcor+1,2*lcor+1,2*lcor+1,2*lcor+1)
1542  complex(dpc), intent(out) :: mat_out_c(2*lcor+1,2*lcor+1,2*lcor+1,2*lcor+1)
1543 
1544 !Local variables ---------------------------------------
1545 !scalars
1546  integer :: gm,hm,jm,gg,hh,ii,jj,ll,mm,im
1547  real(dp),parameter :: invsqrt2=one/sqrt2
1548  real(dp) :: onem
1549  complex(dpc) :: tmp2
1550  character(len=500) :: message
1551 !arrays
1552  complex(dpc),allocatable :: slm2ylm(:,:)
1553 
1554 ! *********************************************************************
1555 
1556 
1557  if (option/=1.and.option/=2) then
1558    message=' option=/1 or 2 !'
1559    ABI_BUG(message)
1560  end if
1561 
1562  if(abs(prtvol)>2) then
1563    write(message,'(3a)') ch10, "   vee_slm2ylm_hu"
1564    call wrtout(std_out,message,'COLL')
1565  end if
1566 
1567  if(abs(prtvol)>2) then
1568    if(option==1) then
1569      write(message,'(3a)') ch10,"matrix in Slm basis is changed into Ylm basis"
1570      call wrtout(std_out,message,'COLL')
1571    else if(option==2) then
1572      write(message,'(3a)') ch10,"matrix in Ylm basis is changed into Slm basis"
1573      call wrtout(std_out,message,'COLL')
1574    end if
1575  end if
1576 
1577  ll=lcor
1578  ABI_MALLOC(slm2ylm,(2*ll+1,2*ll+1))
1579  slm2ylm=czero
1580  mat_out_c=czero
1581 
1582 !  ===== Definitions of slm2ylm
1583  do im=1,2*ll+1
1584    mm=im-ll-1;jm=-mm+ll+1   ! mmj=-mm
1585 !! im is in {1,....2*ll+1}
1586 !! mm is in {-ll,....+ll}
1587 !! jm is in {2*ll+1,....,1}
1588    onem=dble((-1)**mm)
1589    if (mm> 0) then ! im in {ll+1,2ll+1} and jm in {ll+1,1}
1590      slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
1591      slm2ylm(jm,im)= cmplx(invsqrt2,     zero,kind=dp)
1592    end if
1593    if (mm==0) then
1594      slm2ylm(im,im)=cone
1595    end if
1596    if (mm< 0) then
1597      slm2ylm(im,im)= cmplx(zero,     invsqrt2,kind=dp)
1598      slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
1599    end if
1600  end do
1601 ! do im=1,2*ll+1
1602 !   write(message,'(7(2f14.5))') (slm2ylm(im,jm),jm=1,2*ll+1)
1603 !   call wrtout(std_out,message,'COLL')
1604 ! end do
1605 
1606 !  ===== Definitions of slm2ylm
1607 !!!!  pawtab(itypat)%vee(m11,m31,m21,m41)= <m11 m31| vee| m21 m41 >
1608 !!!!  pawtab(itypat)%vee(m11,m21,m31,m41)= <m11 m21| vee| m31 m41 >
1609 
1610  do jm=1,2*ll+1
1611    do im=1,2*ll+1
1612      do hm=1,2*ll+1
1613        do gm=1,2*ll+1
1614          tmp2=czero
1615          do gg=1,2*ll+1
1616            do hh=1,2*ll+1
1617              do ii=1,2*ll+1
1618                do jj=1,2*ll+1
1619                  if(option==1) then
1620                    tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*(slm2ylm(im,ii))*CONJG(slm2ylm(jm,jj))&
1621 &                                                   *(slm2ylm(gm,gg))*CONJG(slm2ylm(hm,hh))
1622 !                   if(gm==1.and.hm==1.and.im==1.and.jm==1) then
1623 !                      write(6,'(4i4,2f10.5,2f10.5)') gg,hh,ii,jj,tmp2,mat_inp_c(gg,hh,ii,jj)
1624 !                      write(6,*) "i1"
1625 !                   endif
1626                  else if(option==2) then
1627                    tmp2=tmp2+mat_inp_c(gg,ii,hh,jj)*CONJG(slm2ylm(ii,im))*(slm2ylm(jj,jm))&
1628 &                                                   *CONJG(slm2ylm(gg,gm))*(slm2ylm(hh,hm))
1629                  end if
1630                end do
1631              end do
1632            end do
1633          end do
1634 !         mat_out_c(gm,hm,im,jm)=tmp2
1635          mat_out_c(gm,im,hm,jm)=tmp2
1636        end do
1637      end do
1638    end do
1639  end do
1640 
1641  ABI_FREE(slm2ylm)
1642 
1643 end subroutine vee_slm2ylm_hu