TABLE OF CONTENTS


ABINIT/m_forctqmc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_forctqmc

FUNCTION

 Prepare CTQMC and call CTQMC

COPYRIGHT

 Copyright (C) 2006-2022 ABINIT group (BAmadon, VPlanes)
 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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_forctqmc
27 
28  use defs_basis
29  use m_errors
30  use m_xmpi
31  use m_abicore
32  use m_Ctqmc
33  use m_CtqmcInterface
34  use m_Ctqmcoffdiag
35  use m_CtqmcoffdiagInterface
36  use m_GreenHyb
37  use m_data4entropyDMFT
38 
39  use m_pawang, only : pawang_type
40  use m_crystal, only : crystal_t
41  use m_green, only : green_type,occup_green_tau,print_green,printocc_green,spline_fct,copy_green,init_green,destroy_green,&
42 & int_fct,greendftcompute_green,fourier_green
43  use m_paw_dmft, only : paw_dmft_type
44  use m_hide_lapack,         only : xginv
45  use m_oper, only : oper_type,destroy_oper,init_oper,inverse_oper
46  use m_self, only : self_type
47  use m_matlu, only : matlu_type,sym_matlu, print_matlu, &
48 & diag_matlu,init_matlu,destroy_matlu,rotate_matlu,checkdiag_matlu,checkreal_matlu, &
49 & copy_matlu, diff_matlu, slm2ylm_matlu, shift_matlu, prod_matlu,fac_matlu,&
50 & add_matlu,printplot_matlu,identity_matlu,zero_matlu
51  use m_hu, only : hu_type,rotatevee_hu,vee_ndim2tndim_hu_r,copy_hu,destroy_hu
52  use m_io_tools, only : flush_unit, open_file
53  use m_datafordmft, only : hybridization_asymptotic_coefficient,compute_levels
54  use m_special_funcs, only : sbf8
55  use m_paw_numeric, only : jbessel=>paw_jbessel
56  use m_paw_correlations, only : calc_vee
57 #ifdef HAVE_NETCDF
58   use netcdf !If calling TRIQS via python invocation, write a .nc file
59 #endif
60 
61  implicit none
62 
63  private
64 
65  public :: qmc_prep_ctqmc
66  public :: testcode_ctqmc
67  public :: testcode_ctqmc_b
68  public :: ctqmcoutput_to_green
69  public :: ctqmcoutput_printgreen
70  public :: ctqmc_calltriqs

m_forctqmc/ctqmc_calltriqs [ Functions ]

[ Top ] [ m_forctqmc ] [ Functions ]

NAME

 ctqmc_calltriqs

FUNCTION

  Call TRIQS solver and perform calculation of Green's function using
  Legendre coefficients.

INPUTS

  paw_dmft <type(paw_dmft_type)>= DMFT data structure
  cryst_struc <type(crystal_t)>=crystal structure data
  hu <type(hu_type)>= U interaction
  levels_ctqmc(nflavor) = atomic levels
  gw_tmp_nd(nb_of_frequency,nflavor,nflavor) = Green's fct in imag freq (with off diag terms)
  gtmp_nd(dmftqmc_l,nflavor,nflavor) = Green's fct in imag time (with off diag terms)
  fw1_nd(dmft_nwlo,nflavor,nflavor) = Hybridization fct in imag time (with off diag terms)
  leg_measure = logical, true is legendre measurement is activated
  iatom= index of atom

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

2432 subroutine ctqmc_calltriqs(paw_dmft,cryst_struc,hu,levels_ctqmc,gtmp_nd,gw_tmp_nd,fw1_nd,leg_measure,iatom)
2433 
2434 #if defined HAVE_TRIQS_v2_0 || defined HAVE_TRIQS_v1_4
2435  use TRIQS_CTQMC !Triqs module
2436 #endif
2437 #if defined HAVE_PYTHON_INVOCATION
2438  use INVOKE_PYTHON
2439 #endif
2440  use, intrinsic :: iso_c_binding
2441 
2442 !Arguments ------------------------------------
2443 !scalars
2444  type(paw_dmft_type), intent(in)  :: paw_dmft
2445  type(crystal_t),intent(in) :: cryst_struc
2446  type(hu_type), intent(in) :: hu(cryst_struc%ntypat)
2447  real(dp), allocatable, target, intent(inout) :: gtmp_nd(:,:,:)
2448  complex(dpc), allocatable, target, intent(inout) :: gw_tmp_nd(:,:,:)
2449  complex(dpc), allocatable, target, intent(in) :: fw1_nd(:,:,:)
2450  real(dp), allocatable, target, intent(inout) ::  levels_ctqmc(:)
2451  logical(kind=1), intent(in) :: leg_measure
2452  integer, intent(in) :: iatom
2453 
2454 !Local variables ------------------------------
2455  complex(dpc), allocatable, target ::fw1_nd_tmp(:,:,:)
2456  complex(dpc), allocatable, target :: g_iw(:,:,:)
2457  real(dp), allocatable, target :: u_mat_ij(:,:)
2458  real(dp), allocatable, target :: u_mat_ijkl(:,:,:,:)
2459  real(dp), allocatable, target :: u_mat_ijkl_tmp(:,:,:,:)
2460  real(dp), allocatable, target :: gl_nd(:,:,:)
2461  type(c_ptr) :: levels_ptr, fw1_nd_ptr, u_mat_ij_ptr, u_mat_ijkl_ptr, g_iw_ptr, gtau_ptr, gl_ptr
2462  real(dp), allocatable :: jbes(:)
2463  character(len=500) :: message
2464  integer :: itau, ifreq, iflavor1
2465  integer :: iflavor2,iflavor,nflavor,iflavor3,itypat
2466  integer :: nfreq,ntau,nleg,ileg
2467  integer :: verbosity_solver ! min 0 -> max 3
2468  logical(kind=1) :: rot_inv = .false.
2469  logical(kind=1) :: hist = .false.
2470  logical(kind=1) :: wrt_files = .true.
2471  logical(kind=1) :: tot_not = .true.
2472  real(dp) :: beta,besp,bespp,xx
2473  complex(dpc) :: u_nl
2474 
2475 !----------
2476 !Variables for writing out the NETCDF file
2477 !----------
2478  integer(kind=4) :: ncid
2479  integer(kind=4) :: dim_one_id, dim_nflavor_id, dim_nwlo_id, dim_nwli_id
2480  integer(kind=4) :: dim_qmc_l_id, dim_nleg_id
2481  integer(kind=4), dimension(2) :: dim_u_mat_ij_id
2482  integer(kind=4), dimension(3) :: dim_fw1_id, dim_g_iw_id, dim_gl_id, dim_gtau_id
2483  integer(kind=4), dimension(4) :: dim_u_mat_ijkl_id
2484  integer(kind=4) :: var_rot_inv_id, var_leg_measure_id, var_hist_id, var_wrt_files_id
2485  integer(kind=4) :: var_tot_not_id, var_n_orbitals_id, var_n_freq_id, var_n_tau_id, var_n_l_id, var_n_cycles_id
2486  integer(kind=4) :: var_cycle_length_id, var_ntherm_id, var_verbo_id, var_seed_id, var_beta_id
2487  integer(kind=4) :: var_levels_id, var_u_mat_ij_id, var_u_mat_ijkl_id, var_real_fw1_nd_id, var_imag_fw1_nd_id
2488  integer(kind=4) :: var_real_g_iw_id, var_imag_g_iw_id, var_gtau_id, var_gl_id, var_spacecomm_id
2489 
2490  integer(kind=4) :: varid
2491  logical :: file_exists
2492  complex :: i
2493  character(len=100) :: filename
2494 
2495  real(dp), allocatable, target :: new_re_g_iw(:,:,:), new_im_g_iw(:,:,:)
2496  real(dp), allocatable, target :: new_g_tau(:,:,:), new_gl(:,:,:)
2497 !----------
2498 ! ************************************************************************
2499 
2500  ! fw1_nd: Hybridation
2501  ! levels_ctqmc: niveaux
2502  ! hu(itypat)%udens(:,:) : U_ij
2503  ! hu(itypat)%u(:,:,:,:) : uijkl
2504  ! temperature : paw_dmft%temp
2505  ! paw_dmft%dmftqmc_l: nombre de points en temps -1
2506  ! paw_dmft%dmftqmc_n: nombre de cycles
2507  ! ?? Quelles sorties: Les fonctions de Green
2508  ! frequence/temps/Legendre.
2509  ! Double occupations ?? <n_i n_j>
2510  ! test n_tau > 2*nfreq => ntau = 2*nfreq + 1
2511  !   for non diagonal code:
2512  !   call CtqmcInterface_run(hybrid,fw1_nd(1:paw_dmft%dmftqmc_l,:,:),Gtau=gtmp_nd,&
2513  !&  Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),&
2514  !&  Noise=Noise,matU=hu(itypat)%udens,opt_levels=levels_ctqmc,hybri_limit=hybri_limit)
2515  !Check choice of user to fix model bool var for the solver
2516  if (paw_dmft%dmft_solv==6) then
2517    rot_inv = .false.
2518  else !obviously paw_dmft%dmft_solv==7 with rot invariant terms
2519    rot_inv = .true.
2520  end if
2521 
2522  nfreq = paw_dmft%dmft_nwli
2523  !paw_dmft%dmft_nwlo = paw_dmft%dmft_nwli !transparent for user
2524  ntau  = paw_dmft%dmftqmc_l !(2*paw_dmft%dmftqmc_l)+1 !nfreq=paw_dmft%dmft_nwli
2525  nleg  = paw_dmft%dmftctqmc_triqs_nleg
2526  nflavor=2*(2*paw_dmft%lpawu(iatom)+1)
2527  itypat=cryst_struc%typat(iatom)
2528 
2529 
2530  verbosity_solver = paw_dmft%prtvol
2531  beta = 1.0/(paw_dmft%temp*Ha_eV)
2532 
2533  !Allocation in/output array phase:
2534  ABI_MALLOC(fw1_nd_tmp,(1:nflavor,1:nflavor,1:nfreq)) !column major
2535  ABI_MALLOC(g_iw,(1:nflavor,1:nflavor,1:nfreq)) !column major
2536  ABI_MALLOC(u_mat_ij,(1:nflavor,1:nflavor)) !column major
2537  ABI_MALLOC(u_mat_ijkl,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major
2538  ABI_MALLOC(u_mat_ijkl_tmp,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major
2539 
2540  if ( leg_measure ) then !only if functionality is enabled
2541    ABI_MALLOC(gl_nd,(1:nleg,1:nflavor,1:nflavor)) !column major !nl = 30 by default
2542  end if
2543 
2544  !Conversion datas Ha -> eV (some duplications for test...)
2545  !fw1_nd_tmp = fw1_nd(1:paw_dmft%dmftqmc_l,:,:) * Ha_eV !fw1_nd = fw1_nd * Ha_eV !Ok?
2546 
2547  do iflavor=1,nflavor
2548    do iflavor1=1,nflavor
2549      do ifreq=1,nfreq
2550        fw1_nd_tmp(iflavor,iflavor1,ifreq) = fw1_nd(ifreq,iflavor,iflavor1) * Ha_eV
2551 !        WRITE(500,*) "[IN Fortran] F[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",fw1_nd(ifreq,iflavor,iflavor1)
2552      end do
2553    end do
2554  end do
2555 
2556     !Report test
2557 !    WRITE(502,*) hu(itypat)%udens
2558 !    do ifreq=1,paw_dmft%dmftqmc_l
2559 !      write(501,*) ((fw1_nd(ifreq,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor)
2560 !    enddo
2561     !write(866,*)paw_dmft%dmft_nwlo,paw_dmft%dmftqmc_l
2562     !write(866,*) u_mat_ij
2563 !   do iflavor=1,nflavor+1
2564 !     do iflavor1=1,nflavor+1
2565 !       WRITE(502,*) "[OUT Fortran] U(i,j)[ l= ",iflavor," l_= ",iflavor1,"] = ",hu(itypat)%udens(iflavor,iflavor1)
2566 !     enddo
2567 !   enddo
2568 
2569 !          if(paw_dmft%myproc==0) then
2570 !          do iflavor=1,nflavor
2571 !            do iflavor1=1,nflavor
2572 !               do iflavor2=1,nflavor
2573 !                  do iflavor3=1,nflavor
2574 !                    write(490,*), hu(itypat)%vee(iflavor,iflavor1,iflavor2,iflavor3)
2575 !                  enddo
2576 !                 enddo
2577 !                enddo
2578 !              enddo
2579 !          endif
2580 
2581 !          if(paw_dmft%myproc==0) then
2582 !          do iflavor=1,nflavor
2583 !            do iflavor1=1,nflavor
2584 !            write(491,*), hu(itypat)%udens(iflavor,iflavor1) !(1,1,1,1)
2585 !            enddo
2586 !          enddo
2587 !          endif
2588 
2589 !          do iflavor=1,nflavor
2590 !            do iflavor1=1,nflavor
2591 !          do iflavor2=1,nflavor
2592 !            do iflavor3=1,nflavor
2593                  ! WRITE(552,*), hu(itypat)%vee!(iflavor,iflavor1,iflavor2,iflavor3)
2594 !            enddo
2595 !          enddo
2596 !            enddo
2597 !          enddo
2598 
2599  call vee_ndim2tndim_hu_r(paw_dmft%lpawu(iatom),hu(itypat)%vee,u_mat_ijkl_tmp,1)
2600  do iflavor=1,nflavor
2601    do iflavor1=1,nflavor
2602      do iflavor2=1,nflavor
2603        do iflavor3=1,nflavor
2604          u_mat_ijkl(iflavor,iflavor1,iflavor2,iflavor3)   =  Ha_eV * u_mat_ijkl_tmp(iflavor,iflavor1,iflavor2,iflavor3)
2605        end do
2606      end do
2607    end do
2608  end do
2609 
2610  !u_mat_ijkl   =  Ha_eV * reshape( u_mat_ijkl , [nflavor,nflavor,nflavor,nflavor] )  !column -> row major + conversion
2611  u_mat_ij     = transpose( hu(itypat)%udens ) * Ha_eV !column -> row major + conversion
2612  levels_ctqmc = levels_ctqmc * Ha_eV
2613 
2614  !Location array in memory for C++ pointer args to pass
2615  !----------------------------------------------------
2616  g_iw_ptr       = C_LOC( gw_tmp_nd ) !C_LOC( g_iw )
2617  gtau_ptr       = C_LOC( gtmp_nd ) !C_LOC( gtau )
2618  gl_ptr         = C_LOC( gl_nd )
2619  fw1_nd_ptr     = C_LOC( fw1_nd_tmp )
2620  u_mat_ij_ptr   = C_LOC( u_mat_ij )
2621  u_mat_ijkl_ptr = C_LOC( u_mat_ijkl )
2622  levels_ptr     = C_LOC( levels_ctqmc )
2623 
2624  !Calling interfaced TRIQS solver subroutine from src/67_triqs_ext package
2625  if (paw_dmft%dmft_solv==9) then
2626 #ifndef HAVE_NETCDF
2627   write(message,'(2a)') ch10,' NETCDF requiered! ABINIT communicates with the python script through netcdf.'
2628   call wrtout(std_out,message,'COLL')
2629   ABI_ERROR(message)
2630 #else
2631 #ifndef HAVE_PYTHON_INVOCATION
2632   write(message,'(23a)') ch10,' Python invocation flag requiered! You need to install ABINIT with ',&
2633    'enable_python_invocation = yes" in your "configure.ac" file.'
2634   call wrtout(std_out,message,'COLL')
2635   ABI_ERROR(message)
2636 #else
2637   ! Creating the NETCDF file
2638   ! write(std_out, "(a)") trim(paw_dmft%filapp)
2639   write(filename, '(a, a)') trim(paw_dmft%filnamei), "_abinit_output_for_py.nc"
2640   write(std_out, '(3a)') ch10, "    Creating NETCDF file: ", trim(filename)
2641   NCF_CHECK(nf90_create(filename, NF90_CLOBBER, ncid))
2642  
2643   ! Defining the dimensions of the variables to write in the NETCDF file
2644   NCF_CHECK(nf90_def_dim(ncid, "one", 1, dim_one_id))
2645   NCF_CHECK(nf90_def_dim(ncid, "nflavor", nflavor, dim_nflavor_id))
2646   NCF_CHECK(nf90_def_dim(ncid, "nwlo", paw_dmft%dmft_nwlo, dim_nwlo_id))
2647   NCF_CHECK(nf90_def_dim(ncid, "nwli", paw_dmft%dmft_nwli, dim_nwli_id))
2648   NCF_CHECK(nf90_def_dim(ncid, "qmc_l", paw_dmft%dmftqmc_l, dim_qmc_l_id))
2649   NCF_CHECK(nf90_def_dim(ncid, "nleg", nleg, dim_nleg_id))
2650  
2651   dim_u_mat_ij_id = (/ dim_nflavor_id, dim_nflavor_id /)
2652   dim_u_mat_ijkl_id = (/ dim_nflavor_id, dim_nflavor_id, dim_nflavor_id, dim_nflavor_id /)
2653   dim_fw1_id = (/ dim_nflavor_id, dim_nflavor_id, dim_nwli_id /)
2654   dim_g_iw_id = (/ dim_nwli_id, dim_nflavor_id, dim_nflavor_id /)
2655   dim_gtau_id = (/ dim_qmc_l_id, dim_nflavor_id, dim_nflavor_id /)
2656   dim_gl_id = (/ dim_nleg_id, dim_nflavor_id, dim_nflavor_id /)
2657  
2658   ! Defining the variables
2659   NCF_CHECK(nf90_def_var(ncid, "rot_inv",         NF90_INT, dim_one_id,           var_rot_inv_id))
2660   NCF_CHECK(nf90_def_var(ncid, "leg_measure",     NF90_INT, dim_one_id,           var_leg_measure_id))
2661   NCF_CHECK(nf90_def_var(ncid, "hist",            NF90_INT, dim_one_id,           var_hist_id))
2662   NCF_CHECK(nf90_def_var(ncid, "wrt_files",       NF90_INT, dim_one_id,           var_wrt_files_id))
2663   NCF_CHECK(nf90_def_var(ncid, "tot_not",         NF90_INT, dim_one_id,           var_tot_not_id))
2664   NCF_CHECK(nf90_def_var(ncid, "n_orbitals",      NF90_INT, dim_one_id,           var_n_orbitals_id))
2665   NCF_CHECK(nf90_def_var(ncid, "n_freq",          NF90_INT, dim_one_id,           var_n_freq_id))
2666   NCF_CHECK(nf90_def_var(ncid, "n_tau",           NF90_INT, dim_one_id,           var_n_tau_id))
2667   NCF_CHECK(nf90_def_var(ncid, "n_l",             NF90_INT, dim_one_id,           var_n_l_id))
2668   NCF_CHECK(nf90_def_var(ncid, "n_cycles",        NF90_INT, dim_one_id,           var_n_cycles_id))
2669   NCF_CHECK(nf90_def_var(ncid, "cycle_length",    NF90_INT, dim_one_id,           var_cycle_length_id))
2670   NCF_CHECK(nf90_def_var(ncid, "ntherm",          NF90_INT, dim_one_id,           var_ntherm_id))
2671   NCF_CHECK(nf90_def_var(ncid, "verbo",           NF90_INT, dim_one_id,           var_verbo_id))
2672   NCF_CHECK(nf90_def_var(ncid, "seed",            NF90_INT, dim_one_id,           var_seed_id))
2673   NCF_CHECK(nf90_def_var(ncid, "beta",            NF90_FLOAT, dim_one_id,         var_beta_id))
2674   NCF_CHECK(nf90_def_var(ncid, "levels",          NF90_DOUBLE, dim_nflavor_id,    var_levels_id))
2675   NCF_CHECK(nf90_def_var(ncid, "u_mat_ij",        NF90_DOUBLE, dim_u_mat_ij_id,   var_u_mat_ij_id))
2676   NCF_CHECK(nf90_def_var(ncid, "u_mat_ijkl",      NF90_DOUBLE, dim_u_mat_ijkl_id, var_u_mat_ijkl_id))
2677   NCF_CHECK(nf90_def_var(ncid, "real_fw1_nd",     NF90_DOUBLE, dim_fw1_id,        var_real_fw1_nd_id))
2678   NCF_CHECK(nf90_def_var(ncid, "imag_fw1_nd",     NF90_DOUBLE, dim_fw1_id,        var_imag_fw1_nd_id))
2679   NCF_CHECK(nf90_def_var(ncid, "real_g_iw",       NF90_DOUBLE, dim_g_iw_id,       var_real_g_iw_id))
2680   NCF_CHECK(nf90_def_var(ncid, "imag_g_iw",       NF90_DOUBLE, dim_g_iw_id,       var_imag_g_iw_id))
2681   NCF_CHECK(nf90_def_var(ncid, "gtau",            NF90_DOUBLE, dim_gtau_id,       var_gtau_id))
2682   NCF_CHECK(nf90_def_var(ncid, "gl",              NF90_DOUBLE, dim_gl_id,         var_gl_id))
2683   NCF_CHECK(nf90_def_var(ncid, "spacecomm",       NF90_INT, dim_one_id,           var_spacecomm_id))
2684   NCF_CHECK(nf90_enddef(ncid))
2685  
2686   ! Filling the variables with actual data
2687   if (rot_inv) then 
2688    NCF_CHECK(nf90_put_var(ncid, var_rot_inv_id,       1))  
2689   else 
2690    NCF_CHECK(nf90_put_var(ncid, var_rot_inv_id,       0))  
2691   end if
2692   if (leg_measure) then
2693    NCF_CHECK(nf90_put_var(ncid, var_leg_measure_id,   1))
2694   else
2695    NCF_CHECK(nf90_put_var(ncid, var_leg_measure_id,   0))
2696   end if
2697   if (hist) then
2698    NCF_CHECK(nf90_put_var(ncid, var_hist_id,          1))
2699   else
2700    NCF_CHECK(nf90_put_var(ncid, var_hist_id,          0))
2701   end if
2702   if (wrt_files) then
2703    NCF_CHECK(nf90_put_var(ncid, var_wrt_files_id,     1))
2704   else
2705    NCF_CHECK(nf90_put_var(ncid, var_wrt_files_id,     0))
2706   end if
2707   if (tot_not) then
2708    NCF_CHECK(nf90_put_var(ncid, var_tot_not_id,       1))
2709   else
2710    NCF_CHECK(nf90_put_var(ncid, var_tot_not_id,       0))
2711   end if
2712   NCF_CHECK(nf90_put_var(ncid, var_n_orbitals_id,         nflavor))
2713   NCF_CHECK(nf90_put_var(ncid, var_n_freq_id,             nfreq))
2714   NCF_CHECK(nf90_put_var(ncid, var_n_tau_id,              ntau))
2715   NCF_CHECK(nf90_put_var(ncid, var_n_l_id,                nleg))
2716   NCF_CHECK(nf90_put_var(ncid, var_n_cycles_id,           int(paw_dmft%dmftqmc_n/paw_dmft%nproc)))
2717   NCF_CHECK(nf90_put_var(ncid, var_cycle_length_id,       paw_dmft%dmftctqmc_meas*2*2*nflavor))
2718   NCF_CHECK(nf90_put_var(ncid, var_ntherm_id,             paw_dmft%dmftqmc_therm))
2719   NCF_CHECK(nf90_put_var(ncid, var_verbo_id,              verbosity_solver))
2720   NCF_CHECK(nf90_put_var(ncid, var_seed_id,               paw_dmft%dmftqmc_seed))
2721   NCF_CHECK(nf90_put_var(ncid, var_beta_id,               beta))
2722   NCF_CHECK(nf90_put_var(ncid, var_levels_id,             levels_ctqmc))
2723   NCF_CHECK(nf90_put_var(ncid, var_u_mat_ij_id,           u_mat_ij))
2724   NCF_CHECK(nf90_put_var(ncid, var_u_mat_ijkl_id,         u_mat_ijkl))
2725   NCF_CHECK(nf90_put_var(ncid, var_real_fw1_nd_id,        real(fw1_nd_tmp)))
2726   NCF_CHECK(nf90_put_var(ncid, var_imag_fw1_nd_id,        aimag(fw1_nd_tmp)))
2727   NCF_CHECK(nf90_put_var(ncid, var_real_g_iw_id,          real(gw_tmp_nd)))
2728   NCF_CHECK(nf90_put_var(ncid, var_imag_g_iw_id,          aimag(gw_tmp_nd)))
2729   NCF_CHECK(nf90_put_var(ncid, var_gtau_id,               gtmp_nd))
2730   NCF_CHECK(nf90_put_var(ncid, var_gl_id,                 gl_nd))
2731   NCF_CHECK(nf90_put_var(ncid, var_spacecomm_id,          paw_dmft%spacecomm))
2732   NCF_CHECK(nf90_close(ncid))
2733 
2734   write(std_out, '(4a)') ch10, "    NETCDF file ", trim(filename), " written; Launching python invocation"
2735  
2736   ! Invoking python to execute the script
2737   ! call Invoke_python_triqs (paw_dmft%myproc, trim(paw_dmft%filnamei)//c_null_char, paw_dmft%spacecomm)
2738   call Invoke_python_triqs (paw_dmft%myproc, trim(paw_dmft%filnamei)//c_null_char)
2739   call xmpi_barrier(paw_dmft%spacecomm)
2740   call flush_unit(std_out)
2741  
2742   ! Allocating the fortran variables for the results
2743   ABI_MALLOC(new_re_g_iw,(nflavor,nflavor, paw_dmft%dmft_nwli))
2744   ABI_MALLOC(new_im_g_iw,(nflavor,nflavor, paw_dmft%dmft_nwli))
2745   ABI_MALLOC(new_g_tau,(nflavor,nflavor, paw_dmft%dmftqmc_l))
2746   ABI_MALLOC(new_gl,(nflavor,nflavor, nleg))
2747   i = (0, 1)
2748   
2749   ! Check if file exists
2750   write(filename, '(a, a)') trim(paw_dmft%filnamei), "_py_output_for_abinit.nc"
2751 
2752   INQUIRE(FILE=filename, EXIST=file_exists)
2753   if(.not. file_exists) then
2754    write(message,'(4a)') ch10,' Cannot find file ', trim(filename), '! Make sure the python script writes it with the right name and at the right place!'
2755    call wrtout(std_out,message,'COLL')
2756    ABI_ERROR(message)
2757   endif
2758 
2759   write(std_out, '(3a)') ch10, "    Reading NETCDF file ", trim(filename)
2760 
2761   ! Opening the NETCDF file
2762   NCF_CHECK(nf90_open(filename, nf90_nowrite, ncid))
2763  
2764   ! Read from the file
2765   ! Re{G_iw}
2766   write(std_out, '(2a)') ch10, "    -- Re[G(iw_n)]"
2767   NCF_CHECK(nf90_inq_varid(ncid, "re_g_iw", varid))
2768   NCF_CHECK(nf90_get_var(ncid, varid, new_re_g_iw))
2769   ! Im{G_iw}
2770   write(std_out, '(2a)') ch10, "    -- Im[G(iw_n)]"
2771   NCF_CHECK(nf90_inq_varid(ncid, "im_g_iw", varid))
2772   NCF_CHECK(nf90_get_var(ncid, varid, new_im_g_iw))
2773   ! G_tau
2774   write(std_out, '(2a)') ch10, "    -- G(tau)"
2775   NCF_CHECK(nf90_inq_varid(ncid, "g_tau", varid))
2776   NCF_CHECK(nf90_get_var(ncid, varid, new_g_tau))
2777   ! G_l
2778   write(std_out, '(2a)') ch10, "    -- G_l"
2779   NCF_CHECK(nf90_inq_varid(ncid, "gl", varid))
2780   NCF_CHECK(nf90_get_var(ncid, varid, new_gl))
2781 
2782   ! Assigning data
2783   do iflavor1=1, nflavor
2784    do iflavor2=1, nflavor
2785     do ifreq=1, paw_dmft%dmft_nwli
2786      gw_tmp_nd(ifreq, iflavor1, iflavor2) = new_re_g_iw(iflavor1, iflavor2, ifreq) &
2787 &               + i*new_im_g_iw(iflavor1, iflavor2, ifreq)
2788     end do
2789     do itau=1, paw_dmft%dmftqmc_l
2790      gtmp_nd(itau, iflavor1, iflavor2) = new_g_tau(iflavor1, iflavor2, itau)
2791     end do
2792     do ileg=1, nleg
2793      gl_nd(ileg, iflavor1, iflavor2) = new_gl(iflavor1, iflavor2, ileg)
2794     end do
2795    end do
2796   end do
2797  
2798   ! Deallocating
2799   ABI_FREE(new_re_g_iw)
2800   ABI_FREE(new_im_g_iw)
2801   ABI_FREE(new_g_tau)
2802   ABI_FREE(new_gl) 
2803 #endif
2804 #endif
2805  elseif(paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
2806   !Calling interfaced TRIQS solver subroutine from src/01_triqs_ext package
2807   !----------------------------------------------------
2808 #if defined HAVE_TRIQS_v2_0 || defined HAVE_TRIQS_v1_4
2809  call Ctqmc_triqs_run (     rot_inv, leg_measure, hist, wrt_files, tot_not,   &
2810 &  nflavor, nfreq, ntau , nleg, int(paw_dmft%dmftqmc_n/paw_dmft%nproc),       &
2811 &  paw_dmft%dmftctqmc_meas*2*2*nflavor, paw_dmft%dmftqmc_therm,               &
2812 &  verbosity_solver, paw_dmft%dmftqmc_seed,beta,                              &
2813 &  levels_ptr,  u_mat_ij_ptr, u_mat_ijkl_ptr, fw1_nd_ptr,                     &
2814 !&  g_iw_ptr, gtau_ptr, gl_ptr, paw_dmft%spacecomm                             )
2815 &  g_iw_ptr, gtau_ptr, gl_ptr, paw_dmft%myproc                             )
2816 #endif
2817  endif
2818 
2819   !WRITE(*,*) "Hello Debug"
2820   !call xmpi_barrier(paw_dmft%spacecomm) !Resynch all processus after calling Impurity solver from TRIQS
2821 
2822   !Report output datas from TRIQS to Abinit
2823   !Interacting G(iw)
2824  do ifreq=1,nfreq
2825    do iflavor1=1,nflavor
2826      do iflavor=1,nflavor
2827     !   gw_tmp_nd(ifreq,iflavor,iflavor1) = g_iw(iflavor,iflavor1,ifreq) !* Ha_eV !because 1/ G0(eV)
2828     !  WRITE(503,*) "[OUT Fortran] G(iw)[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",gw_tmp_nd(ifreq,iflavor,iflavor1)!g_iw(iflavor,iflavor1,ifreq)
2829      end do
2830    end do
2831  end do
2832 
2833 ! Convert in Ha
2834  gw_tmp_nd = gw_tmp_nd*Ha_eV
2835 
2836 !     do iflavor1=1,nflavor
2837 !       do iflavor=1,nflavor
2838 !
2839 !        WRITE(510,*) "[OUT Fortran] U[ l= ",iflavor," l_= ",iflavor1,"] = ",u_mat_ij(iflavor,iflavor1)
2840 !       enddo
2841 !     enddo
2842 
2843 ! if(paw_dmft%myproc==0) write(6,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1)
2844 ! if(paw_dmft%myproc==1) write(6,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1)
2845 ! if(paw_dmft%myproc==0) write(621,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1)
2846 ! if(paw_dmft%myproc==1) write(622,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1)
2847 ! call flush_unit(621)
2848 ! call flush_unit(622)
2849 ! write(message,*) ch10, "essai",paw_dmft%myproc, paw_dmft%myproc,paw_dmft%dmftqmc_seed!gw_tmp_nd(2,1,1)
2850 ! call wrtout(555,message,'PERS',.true.)
2851 ! if(paw_dmft%myproc==0) write(499,*) "essai",paw_dmft%myproc, paw_dmft%dmftqmc_seed
2852 ! if(paw_dmft%myproc==1) write(498,*) "essai",paw_dmft%myproc,paw_dmft%dmftqmc_seed
2853 
2854   !Its associated G(tau): Problem of compatibility => paw_dmft%dmftqmc_l < (2*paw_dmft%dmftqmc_l)+1 => We report only  paw_dmft%dmftqmc_l =  first values of G(tau)...
2855 !   do iflavor=1,nflavor
2856 !     do iflavor1=1,nflavor
2857 !       do itau=1,ntau
2858 !         if ( modulo(itau,2) == 1 ) then !Problem of binding: paw_dmft%dmftqmc_l =! ntau => We take one value by 2 and Write in file all the G(tau) out function from TRIQS
2859           !gtmp_nd(itau,iflavor,iflavor1) = gtau(iflavor,iflavor1,itau)
2860 !         endif
2861 !         if(paw_dmft%myproc==0) then
2862 !           WRITE(504,*) "[OUT Fortran] G[ tau= ",itau," l= ",iflavor," l_= ",iflavor1,"] = ",gtmp_nd(itau,iflavor,iflavor1) !gtmp_nd(itau,iflavor,iflavor1) !passage ok avec ntau/iflavor1/iflavor (iflavor,iflavor1,ntau)
2863 !         endif
2864 !       enddo
2865 !     enddo
2866 !   enddo
2867 
2868   ! Write Legendre Polynoms G(L) for extrapolation of Interacting G(iw) by FT, only if leg_measure == TRUE
2869   ! -------------------------------------------------------------------------------------------
2870  if (leg_measure) then
2871    do ileg=1,nleg
2872      WRITE(505,*) ileg,((gl_nd(ileg,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor)
2873    end do
2874    close(505)
2875  end if
2876 ! f(paw_dmft%myproc==0) then
2877 !  do itau=1,paw_dmft%dmftqmc_l
2878 !    write(490,*) ((gtmp_nd(itau,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor)
2879 !  enddo
2880 ! ndif
2881  ABI_FREE( fw1_nd_tmp )
2882  ABI_FREE( g_iw )
2883  ABI_FREE( u_mat_ijkl )
2884  ABI_FREE( u_mat_ijkl_tmp )
2885  ABI_FREE( u_mat_ij )
2886 
2887 
2888   !  Compute Green's function in imaginary freq using Legendre coefficients
2889   ! -----------------------------------------------------------------------
2890  if (leg_measure) then
2891    call xmpi_barrier(paw_dmft%spacecomm)
2892    call flush_unit(std_out)
2893    write(message,'(2a)') ch10,"    ==  Compute G(iw_n) from Legendre coefficients"
2894    call wrtout(std_out,message,'COLL')
2895    ABI_MALLOC( jbes, (nleg))
2896    gw_tmp_nd=czero
2897 
2898   !   write(77,*) " TEST OF BESSEL S ROUTINES 0 0"
2899 
2900   !   xx=0_dp
2901   !   ileg=0
2902   !   call sbf8(ileg+1,xx,jbes)
2903   !   write(77,*) "T0 A",jbes(ileg+1)
2904   !   call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx)
2905   !   write(77,*) "T0 B",jbes(ileg+1)
2906   !   write(77,*) "T0 C",bessel_jn(ileg,xx)
2907 
2908   !   write(77,*) " TEST OF BESSEL S ROUTINES 1.5 0"
2909 
2910   !   xx=1.5_dp
2911   !   ileg=0
2912   !   call sbf8(ileg+1,xx,jbes)
2913   !   write(77,*) "T1 A",jbes(ileg+1)
2914   !   call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx)
2915   !   write(77,*) "T1 B",jbes(ileg+1)
2916   !   write(77,*) "T1 C",bessel_jn(ileg,xx)
2917 
2918   !   write(77,*) " TEST OF BESSEL S ROUTINES 1.5 1"
2919 
2920   !   xx=1.5_dp
2921   !   ileg=1
2922   !   call sbf8(ileg+1,xx,jbes)
2923   !   write(77,*) "T2 A",jbes(ileg+1)
2924   !   call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx)
2925   !   write(77,*) "T2 B",jbes(ileg+1)
2926   !   write(77,*) "T2 C",bessel_jn(ileg,xx)
2927 
2928 
2929    do ifreq=1,paw_dmft%dmft_nwli
2930      xx=real(2*ifreq-1,kind=dp)*pi/two
2931      if(xx<=100_dp) call sbf8(nleg,xx,jbes)
2932      do ileg=1,nleg
2933     ! write(77,*) "A",ifreq,jbes(ileg),xx
2934 
2935        if(xx>=99) call jbessel(jbes(ileg),besp,bespp,ileg-1,1,xx)
2936     ! write(77,*) "B",ifreq,jbes(ileg),xx
2937 
2938      !write(77,*) "C",ifreq,jbes(ileg),xx
2939 
2940        u_nl=sqrt(float(2*ileg-1))*(-1)**(ifreq-1)*cmplx(0_dp,one)**(ileg)*jbes(ileg)
2941       write(77,*) "----------",ileg,jbes(ileg), u_nl,gl_nd(ileg,1,1)
2942 
2943        do iflavor=1,nflavor
2944          do iflavor1=1,nflavor
2945            gw_tmp_nd(ifreq,iflavor,iflavor1)= gw_tmp_nd(ifreq,iflavor,iflavor1) + &
2946 &           u_nl*gl_nd(ileg,iflavor,iflavor1)
2947          end do
2948        end do
2949 
2950   !    write(77,*) "------------------", gw_tmp_nd(ifreq,1,1)
2951 
2952      end do
2953   !  write(77,*) "------------------ sum ", gw_tmp_nd(ifreq,1,1)
2954    end do
2955    ABI_FREE( jbes )
2956    call xmpi_barrier(paw_dmft%spacecomm)
2957    call flush_unit(std_out)
2958  end if
2959  gw_tmp_nd = gw_tmp_nd*Ha_eV
2960 
2961 
2962  if ( leg_measure ) then !only if functionality is enabled
2963    ABI_FREE(gl_nd)
2964  end if
2965 
2966 
2967 end subroutine ctqmc_calltriqs

m_forctqmc/ctqmcoutput_printgreen [ Functions ]

[ Top ] [ m_forctqmc ] [ Functions ]

NAME

 ctqmcoutput_printgreen

FUNCTION

  Print values of green function in files.
  Symetrize imaginary time Green's function in a peculiar case
  (dmft_solv=8 and natom=1). Should be moved later.

INPUTS

  paw_dmft <type(paw_dmft_type)>= DMFT data structure
  gtmp_nd(dmftqmc_l,nflavor,nflavor) = Green's fct in imag time (with off diag terms)
  gw_tmp_nd(nb_of_frequency,nflavor,nflavor) = Green's fct in imag freq (with off diag terms)
  gtmp(dmftqmc_l,nflavor) = Green's fct in imag time (diag)
  gw_tmp(nb_of_frequency,nflavor+1) =Green's fct in imag freq (diag)
  iatom = atoms on which the calculation has been done

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

2261 subroutine ctqmcoutput_printgreen(paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom)
2262 
2263 !Arguments ------------------------------------
2264 !scalars
2265  type(paw_dmft_type), intent(in)  :: paw_dmft
2266  real(dp), allocatable, intent(inout) :: gtmp_nd(:,:,:)
2267  complex(dpc), allocatable, intent(in) :: gw_tmp(:,:)
2268  complex(dpc), allocatable, intent(in) :: gw_tmp_nd(:,:,:)
2269  real(dp), allocatable, intent(in) :: gtmp(:,:)
2270  integer, intent(in) :: iatom
2271 
2272 !Local variables ------------------------------
2273  character(len=500) :: message
2274  integer :: ifreq, itau,iflavor1
2275  integer :: tndim,iflavor,nflavor
2276  character(len=2) :: gtau_iter,iatomnb
2277  integer :: unt
2278 ! ************************************************************************
2279  tndim=2*paw_dmft%lpawu(iatom)+1
2280  nflavor=2*(tndim)
2281  !----------------------------------------
2282  ! <DEBUG>
2283  !----------------------------------------
2284  ! Construct UNIT
2285  if(paw_dmft%idmftloop < 10) then
2286    write(gtau_iter,'("0",i1)') paw_dmft%idmftloop
2287  elseif(paw_dmft%idmftloop >= 10 .and. paw_dmft%idmftloop < 100) then
2288    write(gtau_iter,'(i2)') paw_dmft%idmftloop
2289  else
2290    gtau_iter="xx"
2291  end if
2292  if(iatom < 10) then
2293    write(iatomnb,'("0",i1)') iatom
2294  elseif(iatom >= 10 .and. iatom < 100) then
2295    write(iatomnb,'(i2)') iatom
2296  else
2297    iatomnb='xx'
2298  end if
2299 
2300  if(paw_dmft%myproc .eq. mod(paw_dmft%nproc+1,paw_dmft%nproc)) then
2301 ! < HACK >
2302   if(paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7) then
2303     if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /=0) then
2304       ABI_ERROR(message)
2305     end if
2306     do ifreq=1,paw_dmft%dmft_nwli
2307       write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),((gw_tmp_nd(ifreq,iflavor,iflavor)), iflavor=1, nflavor)
2308     end do
2309     close(unt)
2310   else
2311     if(paw_dmft%dmft_solv==5) then
2312       if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_"//gtau_iter//".dat", message, newunit=unt) /= 0) then
2313         ABI_ERROR(message)
2314       end if
2315       do itau=1,paw_dmft%dmftqmc_l
2316         write(unt,'(29f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,&
2317         (gtmp(itau,iflavor), iflavor=1, nflavor)
2318       end do
2319       write(unt,'(29f21.14)') 1/paw_dmft%temp, (-1_dp-gtmp(1,iflavor), iflavor=1, nflavor)
2320       close(unt)
2321     endif
2322     if(paw_dmft%dmft_solv==8) then
2323       if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_offdiag_unsym_"//gtau_iter//".dat",&
2324 &      message, newunit=unt) /= 0) then
2325         ABI_ERROR(message)
2326       end if
2327       do itau=1,paw_dmft%dmftqmc_l
2328         write(unt,'(196f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,&
2329         ((gtmp_nd(itau,iflavor,iflavor1), iflavor=1, nflavor),iflavor1=1, nflavor)
2330       end do
2331       close(unt)
2332 !      if(paw_dmft%natom==1) then ! If natom>1, it should be moved outside the loop over atoms
2333 !        ABI_MALLOC(matlu1,(paw_dmft%natom))
2334 !        call init_matlu(paw_dmft%natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,matlu1)
2335 !        do itau=1,paw_dmft%dmftqmc_l
2336 !          do isppol=1,paw_dmft%nsppol
2337 !            do ispinor1=1,paw_dmft%nspinor
2338 !              do im1=1,tndim
2339 !                iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
2340 !                do ispinor2=1,paw_dmft%nspinor
2341 !                  do im2=1,tndim
2342 !                    iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
2343 !                    matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
2344 !&                     gtmp_nd(itau,iflavor1,iflavor2)
2345 !                  end do  ! im2
2346 !                end do  ! ispinor2
2347 !              end do  ! im1
2348 !            end do  ! ispinor
2349 !          end do ! isppol
2350 !          call rotate_matlu(matlu1,eigvectmatlu,paw_dmft%natom,3,0)
2351 !          call slm2ylm_matlu(matlu1,paw_dmft%natom,2,0)
2352 !          call sym_matlu(cryst_struc,matlu1,pawang,paw_dmft)
2353 !          call slm2ylm_matlu(matlu1,paw_dmft%natom,1,0)
2354 !          call rotate_matlu(matlu1,eigvectmatlu,paw_dmft%natom,3,1)
2355 !          do isppol=1,paw_dmft%nsppol
2356 !            do ispinor1=1,paw_dmft%nspinor
2357 !              do im1=1,tndim
2358 !                iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
2359 !                do ispinor2=1,paw_dmft%nspinor
2360 !                  do im2=1,tndim
2361 !                    iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
2362 !                    gtmp_nd(itau,iflavor1,iflavor2)=&
2363 !                     matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
2364 !                  end do  ! im2
2365 !                end do  ! ispinor2
2366 !              end do  ! im1
2367 !            end do  ! ispinor
2368 !          end do ! isppol
2369 !        end do  !itau
2370 !        call destroy_matlu(matlu1,paw_dmft%natom)
2371 !        ABI_FREE(matlu1)
2372 !      endif ! if natom=1
2373       if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_offdiag_"//gtau_iter//".dat",&
2374 &      message, newunit=unt) /= 0) then
2375         ABI_ERROR(message)
2376       end if
2377       do itau=1,paw_dmft%dmftqmc_l
2378         write(unt,'(196f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,&
2379         ((gtmp_nd(itau,iflavor,iflavor1), iflavor=1, nflavor),iflavor1=1, nflavor)
2380       end do
2381       close(unt)
2382     endif
2383     !open(unit=4243, file=trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_F_"//gtau_iter//".dat")
2384     !call BathOperator_printF(paw_dmft%hybrid(iatom)%hybrid%bath,4243) !Already comment here
2385     !close(4243)
2386     if(paw_dmft%dmft_solv==5) then
2387       if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /= 0) then
2388         ABI_ERROR(message)
2389       end if
2390       do ifreq=1,paw_dmft%dmft_nwlo
2391         write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq), &
2392 &        (gw_tmp(ifreq,iflavor), iflavor=1, nflavor)
2393       end do
2394     endif
2395     close(unt)
2396   end if
2397 ! </ HACK >
2398   end if
2399 
2400 
2401 end subroutine ctqmcoutput_printgreen

m_forctqmc/ctqmcoutput_to_green [ Functions ]

[ Top ] [ m_forctqmc ] [ Functions ]

NAME

 ctqmcoutput_to_green

FUNCTION

  Put values of green function from ctqmc into green datatype
  Symetrize over spin if calculation is non magnetic

INPUTS

  paw_dmft <type(paw_dmft_type)>= DMFT data structure
  gtmp_nd(dmftqmc_l,nflavor,nflavor) = Green's fct in imag time (with off diag terms)
  gw_tmp_nd(nb_of_frequency,nflavor,nflavor) = Green's fct in imag freq (with off diag terms)
  gtmp(dmftqmc_l,nflavor) = Green's fct in imag time (diag)
  gw_tmp(nb_of_frequency,nflavor+1) =Green's fct in imag freq (diag)
  iatom = atoms on which the calculation has been done
  leg_measure = logical, to Legendre Measurement or not (if done Green function is frequency is computed)
  opt_nondiag = integer, it activated, then

OUTPUT

  green <type(green_type)>= green's function

SIDE EFFECTS

NOTES

SOURCE

2130 subroutine ctqmcoutput_to_green(green,paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom,leg_measure,opt_nondiag)
2131 
2132 !Arguments ------------------------------------
2133 !scalars
2134  type(paw_dmft_type), intent(in)  :: paw_dmft
2135  type(green_type), intent(inout) :: green
2136  real(dp), allocatable, intent(in) :: gtmp_nd(:,:,:)
2137  complex(dpc), allocatable, intent(in) :: gw_tmp(:,:)
2138  complex(dpc), allocatable, intent(in) :: gw_tmp_nd(:,:,:)
2139  real(dp), allocatable, intent(in) :: gtmp(:,:)
2140  integer, intent(in) :: iatom,opt_nondiag
2141  logical(kind=1), intent(in) :: leg_measure
2142  character(len=500) :: message
2143 
2144 !Local variables ------------------------------
2145  integer :: ifreq, itau,im1,im2,isppol,ispinor1,ispinor2,iflavor1
2146  integer :: iflavor2,tndim,ispinor,iflavor,im,nflavor
2147 ! ************************************************************************
2148  tndim=2*paw_dmft%lpawu(iatom)+1
2149  nflavor=2*(tndim)
2150 
2151  do itau=1,paw_dmft%dmftqmc_l
2152    green%oper_tau(itau)%matlu(iatom)%mat(:,:,:,:,:)=czero
2153  end do
2154  green%occup_tau%matlu(iatom)%mat(nflavor:,:,:,:,:)=czero
2155 
2156  do ifreq=1,paw_dmft%dmft_nwlo
2157    green%oper(ifreq)%matlu(iatom)%mat(:,:,:,:,:)=czero
2158  end do
2159  green%occup%matlu(iatom)%mat(:,:,:,:,:)=czero
2160 
2161 !   built time and frequency green's function from output of CTQMC
2162 ! =================================================================
2163  if(opt_nondiag==1) then
2164    do isppol=1,paw_dmft%nsppol
2165      do ispinor1=1,paw_dmft%nspinor
2166        do im1=1,tndim
2167          iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
2168          do ispinor2=1,paw_dmft%nspinor
2169            do im2=1,tndim
2170              iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
2171              do itau=1,paw_dmft%dmftqmc_l
2172                green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
2173 &               gtmp_nd(itau,iflavor1,iflavor2)
2174                ! symetrize over spin if nsppol=nspinor=1
2175                if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then
2176                  green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
2177 &                 (gtmp_nd(itau,iflavor1,iflavor2)+gtmp_nd(itau,iflavor1+tndim,iflavor2+tndim))/two
2178                end if
2179              end do  !itau
2180              if(paw_dmft%dmft_solv<6.or.leg_measure) then
2181                do ifreq=1,paw_dmft%dmft_nwlo
2182                  green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
2183 &                 gw_tmp_nd(ifreq,iflavor1,iflavor2)
2184                ! symetrize over spin if nsppol=nspinor=1
2185                  if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then
2186                    green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
2187 &                   (gw_tmp_nd(ifreq,iflavor1,iflavor2)+&
2188 &                   gw_tmp_nd(ifreq,iflavor1+tndim,iflavor2+tndim))/two
2189                  end if
2190                end do ! ifreq
2191              end if
2192            end do  ! im2
2193          end do  ! ispinor2
2194        end do  ! im1
2195      end do  ! ispinor
2196    end do ! isppol
2197  else
2198    iflavor=0
2199    do isppol=1,paw_dmft%nsppol
2200      do ispinor=1,paw_dmft%nspinor
2201        do im=1,tndim
2202          iflavor=iflavor+1
2203          do itau=1,paw_dmft%dmftqmc_l
2204            green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gtmp(itau,iflavor)
2205            ! symetrize over spin if nsppol=paw_dmft%nspinor=1
2206            if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then
2207              green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=&
2208 &             (gtmp(itau,iflavor)+gtmp(itau,iflavor+tndim))/two
2209            end if
2210          end do
2211 !       ifreq2=0
2212          do ifreq=1,paw_dmft%dmft_nwlo
2213 !         if(paw_dmft%select_log(ifreq)==1) then
2214 !           ifreq2=ifreq2+1
2215            green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gw_tmp(ifreq,iflavor)
2216            ! symetrize over spin if nsppol=paw_dmft%nspinor=1
2217            if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then
2218              green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=&
2219 &             (gw_tmp(ifreq,iflavor)+gw_tmp(ifreq,iflavor+tndim))/two
2220            end if
2221          end do
2222        end do
2223      end do
2224    end do
2225  end if
2226  if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) then
2227    write(message,'(a,2x,a,f13.5)') ch10,&
2228 &   " == nsppol==1 and nspden==1: Green functions from CTQMC have been symetrized over spin"
2229    call wrtout(std_out,message,'COLL')
2230  end if
2231 
2232 end subroutine ctqmcoutput_to_green

m_forctqmc/qmc_prep_ctqmc [ Functions ]

[ Top ] [ m_forctqmc ] [ Functions ]

NAME

 qmc_prep_ctqmc

FUNCTION

 Prepare and call the qmc subroutines

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  hu <type(hu_type)>= U interaction
  paw_dmft <type(paw_dmft_type)>= DMFT data structure
  pawang <type(pawang)>=paw angular mesh and related data
  pawprtvol = drive the amount of writed data.
  weiss <type(green_type)>= weiss function

OUTPUT

  green <type(green_type)>= green function

NOTES

SOURCE

  97 subroutine qmc_prep_ctqmc(cryst_struc,green,self,hu,paw_dmft,pawang,pawprtvol,weiss)
  98 
  99 !Arguments ------------------------------------
 100 !scalars
 101 ! type(pawang_type), intent(in) :: pawang
 102  type(crystal_t),intent(in) :: cryst_struc
 103  type(green_type), intent(inout) :: green  ! MGNAG: This fix the problem with v7[27:29] on nag@petrus
 104  type(hu_type), intent(in) :: hu(cryst_struc%ntypat)
 105  type(paw_dmft_type), intent(inout)  :: paw_dmft
 106  type(pawang_type), intent(in) :: pawang
 107  integer, intent(in) :: pawprtvol
 108  type(green_type), intent(inout) :: weiss
 109  type(self_type), intent(in) :: self
 110 
 111 !Local variables ------------------------------
 112  character(len=500) :: message
 113  integer :: iatom,ierr,if1,if2,iflavor1,iflavor2,ifreq,im1,ispinor,ispinor1,isppol,itau,itypat,im2,ispinor2
 114  integer :: lpawu,master,mbandc,natom,nflavor,nkpt,nspinor,nsppol,nsppol_imp,tndim,ispa,ispb,ima,imb
 115  integer :: nproc,opt_diag,opt_nondiag,testcode,testrot,dmft_nwlo,opt_fk,useylm,nomega,opt_rot
 116  integer :: ier,rot_type_vee
 117  real(dp) :: f4of2_sla,f6of2_sla
 118  complex(dpc) :: omega_current,integral(2,2)
 119  real(dp) :: doccsum,noise,omega,EE
 120  logical :: nondiaglevels
 121 ! arrays
 122  real(dp), allocatable :: docc(:,:)
 123  real(dp), allocatable, target :: gtmp(:,:), levels_ctqmc(:) !modif
 124  complex(dpc), allocatable :: levels_ctqmc_nd(:,:)
 125  complex(dpc), allocatable :: hybri_limit(:,:)
 126  real(dp), allocatable, target :: gtmp_nd(:,:,:)
 127  real(dp) :: umod(2,2)
 128  complex(dpc), allocatable :: fw1(:,:),gw_tmp(:,:)
 129  complex(dpc), allocatable, target :: gw_tmp_nd(:,:,:) !modif
 130  complex(dpc), allocatable, target :: fw1_nd(:,:,:) !modif
 131  complex(dpc), allocatable :: gw1_nd(:,:,:)
 132  complex(dpc), allocatable :: shift(:)
 133  integer,parameter :: optdb=0
 134  type(coeff2_type), allocatable :: udens_atoms(:)
 135  type(coeff2_type), allocatable :: udens_atoms_for_s(:)
 136 ! Type    -----------------------------------------
 137  type(coeff2c_type), allocatable :: eigvectmatlu(:,:)
 138  type(green_type)  :: weiss_for_rot
 139  type(matlu_type), allocatable :: dmat_diag(:)
 140  type(matlu_type), allocatable :: matlu1(:)
 141  type(matlu_type), allocatable :: matlu2(:)
 142  type(matlu_type), allocatable :: matlu3(:)
 143  type(matlu_type), allocatable :: matlu4(:)
 144  type(matlu_type), allocatable :: identity(:)
 145  type(matlu_type), allocatable :: level_diag(:)
 146  type(oper_type)  :: energy_level
 147  type(hu_type), allocatable :: hu_for_s(:)
 148  !type(self_type) :: self
 149 ! type(green_type) :: gw_loc
 150  type(CtqmcInterface) :: hybrid   !!! WARNING THIS IS A BACKUP PLAN
 151  type(CtqmcoffdiagInterface) :: hybridoffdiag   !!! WARNING THIS IS A BACKUP PLAN
 152  type(green_type) :: greendft
 153  type(matlu_type), allocatable  :: hybri_coeff(:)
 154  integer :: unt,unt2
 155 ! Var added to the code for TRIQS_CTQMC test and default value -----------------------------------------------------------
 156  logical(kind=1) :: leg_measure = .true.
 157 ! ************************************************************************
 158  mbandc=paw_dmft%mbandc
 159  nkpt=paw_dmft%nkpt
 160  nsppol=paw_dmft%nsppol
 161  natom=paw_dmft%natom
 162  nspinor=paw_dmft%nspinor
 163  greendft%whichgreen="DFT"
 164 
 165  call init_green(weiss_for_rot,paw_dmft,opt_oper_ksloc=2)
 166 ! weiss_for_rot=>weiss
 167 ! call init_green(gw_loc,paw_dmft)
 168  call copy_green(weiss,weiss_for_rot,opt_tw=2)
 169 !=======================================================================
 170 !== Use one QMC solver   ===============================================
 171 !=======================================================================
 172  write(message,'(2a)') ch10,'  ===  CT-QMC solver === '
 173  call wrtout(std_out,message,'COLL')
 174 
 175 ! Initialise for compiler
 176  omega_current=czero
 177 
 178 ! Initialise nproc
 179  nproc=paw_dmft%nproc
 180 
 181 ! ======================================
 182 ! Allocations: diagonalization and eigenvectors
 183 ! ======================================
 184  ABI_MALLOC(udens_atoms,(natom))
 185  ABI_MALLOC(eigvectmatlu,(natom,nsppol))
 186  if(paw_dmft%ientropy==1) then
 187    ABI_MALLOC(udens_atoms_for_s,(natom))
 188  endif
 189  ABI_MALLOC(dmat_diag,(natom))
 190  ABI_MALLOC(identity,(natom))
 191  call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,dmat_diag)
 192  call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,identity)
 193  call identity_matlu(identity,natom)
 194  do iatom=1,cryst_struc%natom
 195    lpawu=paw_dmft%lpawu(iatom)
 196    if(lpawu/=-1) then
 197      tndim=nspinor*(2*lpawu+1)
 198      do isppol=1,nsppol
 199        ABI_MALLOC(eigvectmatlu(iatom,isppol)%value,(tndim,tndim))
 200      end do
 201      ABI_MALLOC(udens_atoms(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1)))
 202      if(paw_dmft%ientropy==1) then
 203        ABI_MALLOC(udens_atoms_for_s(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1)))
 204      endif
 205      dmat_diag(iatom)%mat=czero
 206    end if
 207  end do
 208 
 209 ! ___________________________________________________________________________________
 210 !
 211 !  FIRST PART: DIAGONALISATION AND ROTATIONS.
 212 ! ___________________________________________________________________________________
 213 !
 214 
 215 ! =================================================================
 216 ! Choose to diagonalize and how to do it
 217 ! =================================================================
 218 
 219 ! =================================================================
 220 ! Impose diago of density matrix
 221 ! =================================================================
 222 
 223 ! =================================================================
 224 ! Impose diago of levels and Ylm basis if opt_nondiag=1
 225 ! =================================================================
 226 ! opt_diag=1 ! 1: diago the levels (The best choice).
 227 ! opt_diag=2 ! 2: diago density matrix (can be used for historical reasons)
 228 
 229 !  Need in the general case of two input variable for opt_diag and
 230 !  opt_nondiag!
 231 !  The default value of opt_diag should be 2 for historical reasons (or
 232 !  we decide to change the automatic tests)
 233 !  opt_nondiag should be 0 by default
 234  opt_diag    = 1
 235  if(paw_dmft%dmft_solv>=6)  then
 236    opt_nondiag = 1 ! Use cthyb in triqs or ctqmc in abinit with offdiag terms in F
 237  else
 238    opt_nondiag = 0 ! use fast ctqmc in ABINIT without off diagonal terms in F
 239  end if
 240 
 241  useylm=0
 242  if(nspinor==2) then
 243    useylm=1      ! to avoid complex G(tau)
 244  end if
 245 
 246  !write(6,*) "nspinor,useylm",nspinor,useylm
 247  if(useylm==0) then
 248    write(std_out,*) " Slm basis is used (before rotation)"
 249    rot_type_vee=1 ! for rotatevee_hu
 250  else if(useylm==1) then
 251    write(std_out,*) " Ylm basis is used (before rotation)"
 252    rot_type_vee=4 ! for rotatevee_hu
 253  end if
 254 
 255 
 256 ! if(useylm==1.and.opt_diag/=1) ABI_ERROR("useylm==1 and opt_diag/=0 is not possible")
 257  if(hu(1)%jpawu_zero.and.nsppol==2) nsppol_imp=2 ! J=0 and nsppol=2
 258  if(.not.hu(1)%jpawu_zero.or.nsppol/=2) nsppol_imp=1  ! J/=0 ou nsppol=1
 259 ! =================================================================
 260 ! Compute DFT Green's function to compare to weiss_for_rot (check)
 261 ! =================================================================
 262 ! call init_green(greendft,paw_dmft,opt_oper_ksloc=3)
 263 ! call greendftcompute_green(cryst_struc,greendft,pawang,paw_dmft)
 264 !! call copy_green(greendft,weiss_for_rot,2)
 265 
 266 ! =================================================================
 267 ! Compute atomic levels
 268 ! =================================================================
 269  call init_oper(paw_dmft,energy_level,opt_ksloc=3)
 270 
 271  ! Compute atomic levels in Slm basis
 272  ! ----------------------------------
 273  call compute_levels(cryst_struc,energy_level,self%hdc,pawang,paw_dmft,nondiag=nondiaglevels)
 274 
 275  ! If levels are not diagonal, then diagonalize it (according to
 276  ! dmftctqmc_basis)
 277  ! ------------------------------------------------
 278  if(paw_dmft%dmftctqmc_basis==1) then
 279    if(nondiaglevels.or.useylm==1) then
 280      opt_diag=1
 281      write(message,'(3a)') ch10, "   == Hamiltonian in local basis is non diagonal: diagonalise it",ch10
 282    else
 283      opt_diag=0
 284      write(message,'(5a)') ch10, "   == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 &
 285 &     ,"      CTQMC will use this basis",ch10
 286    end if
 287  else if (paw_dmft%dmftctqmc_basis==2) then
 288    if(nondiaglevels.or.useylm==1) then
 289      write(message,'(7a)') ch10, "   == Hamiltonian in local basis is non diagonal",ch10, &
 290 &     "   == According to dmftctqmc_basis: diagonalise density matrix",ch10, &
 291 &     "   == Warning : Check that the Hamiltonian is diagonal !",ch10
 292      opt_diag=2
 293    else
 294      write(message,'(5a)') ch10, "   == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 &
 295 &     ,"      CTQMC will use this basis",ch10
 296      opt_diag=0
 297    end if
 298  else if (paw_dmft%dmftctqmc_basis==0) then
 299    if(nondiaglevels) then
 300      write(message,'(4a)') ch10, "   == Hamiltonian in local basis is non diagonal",ch10, &
 301 &     "   == According to dmftctqmc_basis: keep this non diagonal basis for the calculation"
 302    else
 303      write(message,'(5a)') ch10, "   == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 &
 304 &     ,"      CTQMC will use this basis",ch10
 305    end if
 306    opt_diag=0
 307  end if
 308  call wrtout(std_out,message,'COLL')
 309  if(opt_diag==1) then
 310    write(std_out,*) "  ==  The atomic levels are diagonalized"
 311  else if(opt_diag==2) then
 312    write(std_out,*) "  ==  The correlated occupation matrix is diagonalized"
 313  end if
 314 
 315 ! =================================================================
 316 ! Now, check if diagonalisation is necessary
 317 ! =================================================================
 318 
 319 
 320 ! =================================================================
 321 ! First rotate to Ylm basis the atomic levels
 322 ! =================================================================
 323 
 324  if(useylm==1) then
 325 
 326    ! Rotate from Slm to Ylm the atomic levels
 327    ! ----------------------------------------
 328    call slm2ylm_matlu(energy_level%matlu,natom,1,pawprtvol)
 329 
 330    ! Print atomic energy levels in Ylm basis
 331    ! --------------------------------
 332    if(pawprtvol>=3) then
 333      write(message,'(a,a)') ch10, " == Print Energy levels in Ylm basis"
 334      call wrtout(std_out,message,'COLL')
 335      call print_matlu(energy_level%matlu,natom,1)
 336    end if
 337 
 338  end if ! useylm
 339 
 340 ! ===========================================================================================
 341 ! Start for diagonalization of levels/density matrix according to opt_diag
 342 ! ===========================================================================================
 343  !opt_rot=2 ! do it one time before CTQMC
 344  opt_rot=1 ! do all the rotations successively on all different quantities.
 345  if(opt_diag==1.or.opt_diag==0) then
 346 
 347 
 348    if(opt_diag==1) then
 349 ! =================================================================
 350 ! Diagonalize atomic levels
 351 ! =================================================================
 352      ABI_MALLOC(level_diag,(natom))
 353      call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,level_diag)
 354 
 355      ! Diagonalise atomic levels (opt_real is necessary, because
 356      ! rotation must be real in order for the occupations and Green's
 357      ! function to be real)
 358      ! ---------------------------------------------------------------
 359      call diag_matlu(energy_level%matlu,level_diag,natom,&
 360 &     prtopt=pawprtvol,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,optreal=1,&
 361 &     test=paw_dmft%dmft_solv)  ! temporary: test should be extended to all cases.
 362 
 363 !     call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1)
 364 !       write(message,'(a,2x,a,f13.5)') ch10,&
 365 !&       " == Print first Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie
 366 !       call wrtout(std_out,message,'COLL')
 367 !       call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1)
 368 
 369      if(opt_rot==1) call copy_matlu(level_diag,energy_level%matlu,natom)
 370 
 371 
 372      call destroy_matlu(level_diag,natom)
 373      ABI_FREE(level_diag)
 374 
 375      ! Print diagonalized levels
 376      ! --------------------------
 377      if(pawprtvol>=3) then
 378        write(message,'(a,2x,a,f13.5)') ch10,&
 379 &       " == Print Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie
 380        call wrtout(std_out,message,'COLL')
 381        call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1)
 382      else
 383        write(message,'(a,2x,a,f13.5)') ch10,&
 384 &       " == Energy levels Diagonalized for Fermi Level=",paw_dmft%fermie
 385        call wrtout(std_out,message,'COLL')
 386      end if
 387 
 388      call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee)
 389 
 390 
 391    else if (opt_diag==0) then
 392      do iatom=1,cryst_struc%natom
 393        lpawu=paw_dmft%lpawu(iatom)
 394        itypat=cryst_struc%typat(iatom)
 395        if(lpawu/=-1) then
 396        !  write(6,*) size(udens_atoms(iatom)%value)
 397        !  write(6,*) size(hu(itypat)%udens)
 398        !  write(6,*) udens_atoms(iatom)%value
 399        !  write(6,*) hu(itypat)%udens
 400          udens_atoms(iatom)%value=hu(itypat)%udens
 401        end if
 402      end do
 403    end if
 404   ! call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms)
 405 
 406  else if(opt_diag==2) then
 407 ! =================================================================
 408 ! Diagonalizes density matrix and keep eigenvectors in eigvectmatlu
 409 ! =================================================================
 410 
 411    ! Print density matrix before diagonalization
 412    ! -------------------------------------------
 413    if(pawprtvol>=3) then
 414      write(message,'(a,2x,a)') ch10,        " == Density Matrix before diagonalisation ="
 415      call wrtout(std_out,message,'COLL')
 416      !MGNAG: This call is wrong if green has intent(out), now we use intent(inout)
 417      call print_matlu(green%occup%matlu,natom,1)
 418    end if
 419 
 420 !!  checkstop: we can have two different diagonalisation basis for the up and dn
 421 !!  but one use the same basis, unless the error is really to large(>0.1)
 422 
 423    ! Diagonalize density matrix
 424    ! ---------------------------
 425    call diag_matlu(green%occup%matlu,dmat_diag,natom,&
 426 &   prtopt=4,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,checkstop=.false.)
 427 
 428    ! Print diagonalized density matrix
 429    ! ----------------------------------
 430    if(pawprtvol>=3) then
 431      write(message,'(a,2x,a)') ch10,&
 432 &     " == Diagonalized Density Matrix in the basis used for QMC ="
 433      call wrtout(std_out,message,'COLL')
 434      call print_matlu(dmat_diag,natom,1)
 435 
 436      !write(message,'(2a,i3,13x,a)') ch10,'    ==  Rotation of interaction matrix =='
 437      !call wrtout(std_out,message,'COLL')
 438    end if
 439 
 440    !if (.not.hu(1)%jpawu_zero) &
 441    !ABI_WARNING("In qmc_prep_ctqmc J/=0 and rotation matrix not rotated")
 442 !  Rotate interaction.
 443 !   call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms)
 444    call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee)
 445 
 446 
 447  end if
 448 ! ===========================================================================================
 449 ! END Of diagonalization
 450 ! ===========================================================================================
 451      if(paw_dmft%ientropy==1) then
 452        ABI_MALLOC(hu_for_s,(cryst_struc%ntypat))
 453        ! Usefull to compute interaction energy for U=1 J=J/U when U=0.
 454        call copy_hu(cryst_struc%ntypat,hu,hu_for_s)
 455        f4of2_sla=-1_dp
 456        f6of2_sla=-1_dp
 457        do itypat=1,cryst_struc%ntypat
 458          call calc_vee(f4of2_sla,f6of2_sla,paw_dmft%j_for_s/paw_dmft%u_for_s,hu_for_s(itypat)%lpawu,pawang,one,hu_for_s(itypat)%vee)
 459        enddo
 460        call rotatevee_hu(cryst_struc,hu_for_s,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms_for_s,rot_type_vee)
 461        call destroy_hu(hu_for_s,cryst_struc%ntypat,paw_dmft%dmftqmc_t2g,paw_dmft%dmftqmc_x2my2d)
 462 !      udens_atoms_for_s will be used later.
 463        ABI_FREE(hu_for_s)
 464      endif
 465 
 466 
 467  call flush_unit(std_out)
 468 
 469 ! ===========================================================================================
 470 ! Broadcast matrix of rotation from processor 0 to the other
 471 ! In case of degenerate levels, severals rotations are possible. Here we
 472 ! choose the rotation of proc 0. It is arbitrary.
 473 ! ===========================================================================================
 474  do iatom=1,cryst_struc%natom
 475    lpawu=paw_dmft%lpawu(iatom)
 476    if(lpawu/=-1) then
 477      tndim=nspinor*(2*lpawu+1)
 478      do isppol=1,nsppol
 479        call xmpi_bcast(eigvectmatlu(iatom,isppol)%value,0,paw_dmft%spacecomm,ier)
 480      end do
 481    end if
 482  end do
 483 
 484 
 485      !unitnb=300000+paw_dmft%myproc
 486      !call int2char4(paw_dmft%myproc,tag_proc)
 487      !tmpfil = 'eigvectmatluaftermpi'//tag_proc
 488      !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
 489      !do iflavor1=1,14
 490      !  do iflavor2=1,14
 491      !    write(unitnb,*) iflavor1,iflavor2,eigvectmatlu(1,1)%value(iflavor1,iflavor2)
 492      !  enddo
 493      !enddo
 494 
 495 ! ===========================================================================================
 496 ! Now rotate various quantities in the new basis
 497 ! ===========================================================================================
 498 
 499 !=======================================================
 500 ! Allocate, Compute, and Rotate atomic levels for CTQMC
 501 !=======================================================
 502 
 503    ! If levels not rotated, rotate them
 504    ! -----------------------------------
 505  if(opt_diag==2.and.opt_rot==1) call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1)
 506 
 507    ! Print atomic levels
 508    ! -------------------
 509  if(pawprtvol>=3) then
 510    write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels after rotation"
 511    call wrtout(std_out,message,'COLL')
 512    call print_matlu(energy_level%matlu,natom,1)
 513  else
 514    write(message,'(a,2x,a,f13.5)') ch10," == CT-QMC Energy levels rotated"
 515    call wrtout(std_out,message,'COLL')
 516  end if
 517 
 518 !====================================================================
 519 ! If levels were diagonalized before, then rotate density matrix for
 520 ! information.
 521 !====================================================================
 522  if(opt_diag==1) then
 523    ABI_MALLOC(matlu1,(natom))
 524    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
 525    call copy_matlu(green%occup%matlu,matlu1,natom)
 526    if(pawprtvol>=3) then
 527      write(message,'(a,2x,a)') ch10,&
 528 &     " == Occupations before rotations"
 529      call wrtout(std_out,message,'COLL')
 530      call print_matlu(green%occup%matlu,natom,1)
 531    end if
 532 
 533    ! 1) rotate density matrix to Ylm basis
 534    ! --------------------------------------
 535    if(useylm==1) then
 536      call slm2ylm_matlu(matlu1,natom,1,pawprtvol)
 537      if(pawprtvol>=3) then
 538        write(message,'(a,a)') ch10, " == Print occupations in Ylm basis"
 539        call wrtout(std_out,message,'COLL')
 540        call print_matlu(matlu1,natom,1)
 541      end if
 542    end if
 543 
 544    ! 2) rotate density matrix to rotated basis
 545    ! -------------------------------------------
 546    if(opt_rot==1.or.opt_rot==2) call rotate_matlu(matlu1,eigvectmatlu,natom,3,1)
 547    write(message,'(a,2x,a,f13.5)') ch10," == Rotated occupations (for information)"
 548    call wrtout(std_out,message,'COLL')
 549    call print_matlu(matlu1,natom,1,compl=1)
 550    call checkreal_matlu(matlu1,natom,tol10)
 551    call destroy_matlu(matlu1,natom)
 552    ABI_FREE(matlu1)
 553 
 554  end if
 555 
 556  call flush_unit(std_out)
 557 
 558 
 559 ! =================================================================
 560 ! Rotate weiss function according to eigenvectors.
 561 ! =================================================================
 562 !!!stop
 563   ! Rotate Weiss function first in Ylm basis
 564   ! -------------------------------------------------------------------
 565  if(useylm==1) then
 566    write(message,'(a,2x,a)') ch10, " == Rotation of weiss and greendft in the Ylm Basis="
 567    call wrtout(std_out,message,'COLL')
 568    do ifreq=1,paw_dmft%dmft_nwlo
 569      call slm2ylm_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,1,0)
 570      call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,1,0)
 571      ! call slm2ylm_matlu(greendft%oper(ifreq)%matlu,natom,1,0)
 572    end do
 573  end if
 574 
 575  if(pawprtvol>=3) then
 576    !   write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 577    !   " == Print weiss for small freq 1 before rot" ! debug
 578    !   call wrtout(std_out,message,'COLL') ! debug
 579    !   call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) !  debug
 580 
 581     ! Print Weiss function
 582     ! --------------------
 583    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 584 &  " == Print weiss for 1st freq before rot" ! debug
 585    call wrtout(std_out,message,'COLL') ! debug
 586    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) !  debug
 587    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 588 &  " == Print weiss for last freq before rot" ! debug
 589    call wrtout(std_out,message,'COLL') ! debug
 590    call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) !  debug
 591 !    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 592 !&   " == Print DFT G for 1st freq before rot" ! debug
 593 !    call wrtout(std_out,message,'COLL') ! debug
 594 !    call print_matlu(greendft%oper(1)%matlu,natom,1,compl=1,opt_exp=2) !  debug
 595 !    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 596 !&   " == Print DFT G for last freq before rot" ! debug
 597 !    call wrtout(std_out,message,'COLL') ! debug
 598 !    call print_matlu(greendft%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) !  debug
 599  end if
 600 
 601  if(opt_diag/=0) then
 602    ! Rotate Weiss function from the Slm (or Ylm) to the basis of diagonalisation
 603    ! -------------------------------------------------------------------
 604    write(message,'(a,2x,a)') ch10, " == Rotation of weiss ="
 605    call wrtout(std_out,message,'COLL')
 606    do ifreq=1,paw_dmft%dmft_nwlo
 607      if(opt_rot==1) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 608      if(opt_rot==1) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 609 !    call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6)
 610    end do
 611 
 612    if(paw_dmft%myproc .eq. mod(nproc+1,nproc)) then
 613      if (open_file(trim(paw_dmft%filapp)//"_atom__G0w_.dat", message, newunit=unt) /= 0) then
 614        ABI_ERROR(message)
 615      end if
 616      do ifreq=1,paw_dmft%dmft_nwlo
 617        write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),&
 618 &       (((weiss_for_rot%oper(ifreq)%matlu(1)%mat(im1,im1,isppol,ispinor,ispinor),&
 619 &       im1=1,3),ispinor=1,nspinor),isppol=1,nsppol)
 620      end do
 621      close(unt)
 622    end if
 623 
 624    call flush_unit(std_out)
 625    if(pawprtvol>=3) then
 626      write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 627 &    " == Print weiss for small freq 1 after rot" ! debug
 628      call wrtout(std_out,message,'COLL') ! debug
 629      call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) !  debug
 630      write(message,'(a,2x,a,f13.5)') ch10,&   ! debug
 631 &    " == Print weiss for last freq after rot"   ! debug
 632      call wrtout(std_out,message,'COLL')   ! debug
 633      call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug
 634    end if
 635 
 636 !   ! Rotate DFT Green's function first in Ylm basis then in the rotated basis and compare to weiss_for_rot
 637 !   ! -----------------------------------------------------------------------------------------------------
 638 !   write(message,'(a,2x,a)') ch10, " == Rotation of greendft ="
 639 !   call wrtout(std_out,message,'COLL')
 640 !   do ifreq=1,paw_dmft%dmft_nwlo
 641 !     if(opt_rot==1) call rotate_matlu(greendft%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 642 !     call diff_matlu("Weiss_for_rot","greendft",weiss_for_rot%oper(ifreq)%matlu,greendft%oper(ifreq)%matlu,natom,1,tol14)
 643 !!    call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6)
 644 !   end do
 645 !   if(pawprtvol>=3) then
 646 !     write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 647 !&    " == Print greendft for small freq 1 after rot" ! debug
 648 !     call wrtout(std_out,message,'COLL') ! debug
 649 !     call print_matlu(greendft%oper(1)%matlu,natom,1,compl=1,opt_exp=2) !  debug
 650 !     write(message,'(a,2x,a,f13.5)') ch10,&   ! debug
 651 !&    " == Print greendft for last freq after rot"   ! debug
 652 !     call wrtout(std_out,message,'COLL')   ! debug
 653 !     call print_matlu(greendft%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) ! debug
 654 !   end if
 655 !   call flush_unit(std_out)
 656  end if
 657 
 658 ! =================================================================
 659 ! Compute analytic limit of hybridization and rotate it
 660 ! =================================================================
 661  ABI_MALLOC(hybri_coeff,(paw_dmft%natom))
 662  call init_matlu(paw_dmft%natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,hybri_coeff)
 663  !write(6,*)"hybri1",hybri_coeff(1)%mat(1,1,1,1,1),paw_dmft%natom,cryst_struc%natom
 664 
 665  ! Compute analytical C_ij such that F_ij -> C_ij/iw_n
 666  ! ---------------------------------------
 667  call hybridization_asymptotic_coefficient(cryst_struc,paw_dmft,pawang,hybri_coeff)
 668  write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n for large frequency"
 669  call wrtout(std_out,message,'COLL')
 670 
 671  ! Print analytical C_ij (not rotated)
 672  ! ---------------------------------------
 673  call print_matlu(hybri_coeff,natom,1)
 674 
 675  ! Rotate analytical C_ij in Ylm basis
 676  ! ---------------------------------------
 677  if(useylm==1) call slm2ylm_matlu(hybri_coeff,natom,1,pawprtvol)
 678  if(opt_diag/=0)  then
 679 
 680  ! Rotate analytical C_ij in rotated basis
 681  ! ---------------------------------------
 682    if(opt_rot==1.or.opt_rot==2) call rotate_matlu(hybri_coeff,eigvectmatlu,natom,3,1)
 683 
 684  ! Print analytical C_ij (rotated)
 685  ! ---------------------------------------
 686    write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n after rotation"
 687    call wrtout(std_out,message,'COLL')
 688    call print_matlu(hybri_coeff,natom,1,compl=1,opt_exp=1)
 689  end if
 690 
 691 ! =================================================================
 692 ! Check if rotation is properly done.
 693 ! =================================================================
 694  if(3==4) then
 695    write(message,'(a,2x,a)') ch10,&
 696 &   " == Print  dmat before rot"
 697    call wrtout(std_out,message,'COLL')
 698    call print_matlu(green%occup%matlu,natom,1)
 699    if(useylm==1) call slm2ylm_matlu(green%occup%matlu,natom,1,pawprtvol)
 700    if(opt_rot==1) call rotate_matlu(green%occup%matlu,eigvectmatlu,natom,3,1)
 701    write(message,'(a,2x,a)') ch10,&
 702 &   " == Print  dmat after rot"
 703    call wrtout(std_out,message,'COLL')
 704    call print_matlu(green%occup%matlu,natom,1)
 705 
 706    write(message,'(2a)') ch10,' QMC STOP: DEBUG'
 707    call wrtout(std_out,message,'COLL')
 708    ABI_ERROR(message)
 709  end if
 710 ! =================================================================
 711 ! Check
 712 ! =================================================================
 713 
 714 ! write(message,'(a,2x,a,f13.5)') ch10,&
 715 !&   " == Print weiss for small tau"
 716 ! call wrtout(std_out,message,'COLL')
 717 ! call print_matlu(weiss%oper(1)%matlu,natom,1)
 718 ! write(message,'(a,2x,a,f13.5)') ch10,&
 719 !&   " == Print weiss for large tau"
 720 ! call wrtout(std_out,message,'COLL')
 721 ! call print_matlu(weiss%oper(paw_dmft%dmft_nwlo)%matlu,natom,1)
 722 ! call flush_unit(std_out)
 723 ! write(message,'(2a)') ch10,' Check weiss_for_rot(last freq)'
 724 ! call wrtout(std_out,message,'COLL')
 725 ! call checkdiag_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,tol6,opt=nspinor)
 726 ! call flush_unit(std_out)
 727 ! write(message,'(2a)') ch10,' Check weiss_for_rot(ifreq=1)'
 728 ! call wrtout(std_out,message,'COLL')
 729 ! call checkdiag_matlu(weiss_for_rot%oper(1)%matlu,natom,tol6,opt=nspinor)
 730 ! call flush_unit(std_out)
 731 
 732  master=0
 733 
 734 ! =================================================================
 735 ! Print out
 736 ! =================================================================
 737 
 738 ! Print Weiss
 739 ! -------------
 740  if(paw_dmft%dmft_prgn==1) then
 741    call print_green('Weiss_diag',weiss_for_rot,1,paw_dmft,pawprtvol=1,opt_wt=1,opt_decim=1)
 742  end if
 743 
 744  write(message,'(a,2x,a,f13.5)') ch10,&
 745 & " == Preparing data for CTQMC"
 746  call wrtout(std_out,message,'COLL')
 747 
 748 ! Print Rotate Weiss for 1st and last frequencies
 749 ! ------------------------------------------------
 750  if (pawprtvol>=3) then
 751    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 752 &  " == Print rotated weiss function for small freq in the rotated basis"  ! debug
 753    call wrtout(std_out,message,'COLL')  ! debug
 754    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1)  ! debug
 755    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 756 &  " == Print rotated weiss function for largest freq in the rotated basis"  ! debug
 757    call wrtout(std_out,message,'COLL')  ! debug
 758    call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1)  ! debug
 759  end if
 760 
 761 ! =================================================================
 762 !  VARIABLES FOR CTQMC TESTS
 763  testcode = 0
 764  testrot  = 0
 765  opt_fk=0 ! for developpers to check Fourier transform and computes G0(tau)
 766  opt_fk=1 ! usual case: for real calculations
 767 ! =================================================================
 768 
 769 ! ___________________________________________________________________________________
 770 !
 771 !  SECOND PART : BUILT HYBRIDIZATION FROM G0
 772 ! ___________________________________________________________________________________
 773 !
 774 ! ===========================================================================================
 775 ! Compute inverse of weiss  and compute hybridization
 776 ! ===========================================================================================
 777 
 778 ! Compute inverse of weiss  for each Frequency
 779 ! ----------------------------------------------
 780  do ifreq=1,paw_dmft%dmft_nwlo
 781    ABI_MALLOC(matlu1,(natom))
 782    ABI_MALLOC(matlu2,(natom))
 783    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
 784    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2)
 785 
 786    call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,natom)
 787 
 788    ! Print G_0(iw_n)
 789    ! ----------------
 790    if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"go",60000,imre=1)
 791 
 792    ! Compute G_0^-1
 793    ! -------------------------------------------
 794    ! if opt_fk=1 or testcode/=0  Do the inversion
 795    ! if opt_fk=0                 Do not inverse.
 796    ! If testcode=2 and opt_fk=0  Do the inversion
 797    ! If testcode=1 and opt_fk=0  Do the inversion but no effect, because it will nevertheless be erased
 798    ! If opt_fk=1                 Do the inversion
 799    ! -------------------------------------------
 800    if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"weiss",12000,imre=1)
 801    if(opt_fk==1.or.testcode/=0) call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1)
 802 
 803    ! Print G_0^-1(iw_n)
 804    ! ----------------
 805    if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"goinv",70000,imre=1)
 806 
 807    if(pawprtvol>=4.or.ifreq==paw_dmft%dmft_nwlo) then
 808      if(opt_fk==1.or.testcode/=0) then
 809       ! Check inversion : do the product
 810       ! ----------------------------------------------
 811        call prod_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,matlu2,natom)
 812        write(message,'(a,2x,a,i7)') ch10,&  ! debug
 813 &      " == Print product of  weiss times invers for freq",ifreq
 814        call wrtout(std_out,message,'COLL')  ! debug
 815        call print_matlu(matlu2,natom,1)  ! debug
 816      end if
 817    end if
 818 
 819    call destroy_matlu(matlu1,natom)
 820    call destroy_matlu(matlu2,natom)
 821    ABI_FREE(matlu1)
 822    ABI_FREE(matlu2)
 823  end do
 824 
 825  ! Copy weiss_for_rot into weiss
 826  ! -------------------------------
 827  !call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,weiss%oper(ifreq)%matlu,natom)
 828 
 829 
 830  ! Print G_0^-1 for 1st and last frequencies.
 831  ! -----------------------------------------
 832  if(pawprtvol>=3) then
 833    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 834 &  " == Print G_0^-1 for small freq in the rotated basis"  ! debug
 835    call wrtout(std_out,message,'COLL')  ! debug
 836    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1)  ! debug
 837    write(message,'(a,2x,a,e18.10,a)') ch10,&   ! debug
 838 &  " == Print G_0^-1 for last freq in the rotated basis (last freq=", paw_dmft%omega_lo(paw_dmft%dmft_nwlo),")"  ! debug
 839    call wrtout(std_out,message,'COLL')   ! debug
 840    call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug
 841  end if
 842 
 843 ! Substract frequency from diagonal part
 844 ! ======================================
 845 
 846  ABI_MALLOC(shift,(natom))
 847  do ifreq=1,paw_dmft%dmft_nwlo
 848    shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
 849 
 850 !  write(5555,'(400e17.4)') paw_dmft%omega_lo(ifreq),((((((weiss_for_rot%oper(ifreq)%matlu(1)%mat&
 851 !  & (im,im1,isppol,ispinor,ispinor1)-cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)),im=1,2*3+1),&
 852 !&      im1=1,2*3+1),isppol=1,nsppol),ispinor=1,nspinor),ispinor1=1,nspinor)
 853 
 854   ! Compute G_0^-1-iw_n
 855   ! --------------------
 856    if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift)
 857 
 858   ! Compute -G_0^-1+iw_n
 859   ! --------------------
 860    if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone)
 861 
 862   ! Print -G_0^-1+iw_n
 863   ! --------------------
 864    if(optdb==1) then
 865      call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"G0inv_minus_omega",20000,imre=1)
 866    end if
 867  end do
 868 
 869  ! Print -G_0^+1-iw_n=(F-levels) for last freq in the rotated basis"
 870  ! ------------------------------------------------------------------
 871  ABI_FREE(shift)
 872  if(pawprtvol>=3) then
 873    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 874 &  " == Print G_0^-1-iw_n=-(F-levels) for last freq in the rotated basis"  ! debug
 875    call wrtout(std_out,message,'COLL')   ! debug
 876    call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug
 877  end if
 878 
 879 ! Check numerical limit of F(i_wn)*iw_n (can be used also to compute F )
 880 ! ======================================
 881 
 882  if(opt_nondiag==1) then
 883    ABI_MALLOC(matlu1,(natom))
 884    ABI_MALLOC(matlu2,(natom))
 885    ABI_MALLOC(matlu3,(natom))
 886    ABI_MALLOC(matlu4,(natom))
 887    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
 888    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2)
 889    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu3)
 890    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu4)
 891 
 892    write(message,'(a,2x,a,f13.5)') ch10,  " == energy_levels"
 893    call wrtout(std_out,message,'COLL')
 894    call print_matlu(energy_level%matlu,natom,1,opt_exp=2,compl=1)
 895 
 896    do ifreq=paw_dmft%dmft_nwlo,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency.
 897    !do ifreq=paw_dmft%dmftqmc_l,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency.
 898       ! Compute F (substract levels) for max frequency
 899       ! -----------------------------------------------
 900      call add_matlu(weiss_for_rot%oper(ifreq)%matlu,energy_level%matlu,matlu1,natom,-1)
 901 
 902       ! Print F(iw_n)=-(G_0^-1-iw_n+levels)  for last frequency.
 903       ! --------------------------------------------------------
 904      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
 905        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
 906 &       " == Print F(iw_n)=-(G_0^-1-iw_n+levels) for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
 907        call wrtout(std_out,message,'COLL')
 908        call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 909      end if
 910      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"Hybridization",10000,imre=1)
 911 
 912       ! Put F in weiss_for_rot -> CTQMC
 913       ! -------------------------------
 914      if(opt_rot==2) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 915 !   The following line will produce directly the weiss function for the CTQMC code
 916      if(opt_fk==1) call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom)
 917 
 918 
 919       ! Multiply F by frequency
 920       ! ------------------------
 921      call copy_matlu(matlu1,matlu2,natom)
 922      call fac_matlu(matlu1,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp))
 923      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
 924        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
 925 &       " == Print numerical C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
 926        call wrtout(std_out,message,'COLL')
 927        call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 928      end if
 929      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij",72800,imre=1)
 930      !call rotate_matlu(matlu1,eigvectmatlu,natom,3,1)
 931 
 932      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
 933        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
 934 &       " == Print numerical after back rotation C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
 935        call wrtout(std_out,message,'COLL')
 936        call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 937      end if
 938      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_rotated",72900,imre=1)
 939 
 940       ! Built C_ij/iw_n
 941       ! ------------------------
 942      call copy_matlu(hybri_coeff,matlu1,natom)
 943      call fac_matlu(matlu1,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp))
 944      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_over_omega",72000)
 945     ! if(ifreq==paw_dmft%dmft_nwlo) then
 946     !   write(message,'(a,2x,a,f13.5)') ch10,  " == Print numerical C_ij/iw_n for frequency",paw_dmft%omega_lo(ifreq)
 947     !   call wrtout(std_out,message,'COLL')
 948     !   call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 949     ! endif
 950 
 951       ! For test: put C_ij/i_wn into weiss_for_rot
 952       ! --------------------------------------------
 953      !call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1)
 954 
 955       ! Compute Hybri - C_ij/iw_n
 956       ! ------------------------
 957      call add_matlu(matlu2,matlu1,matlu3,natom,-1)
 958 
 959       ! Print Hybri - C_ij/iw_n
 960       ! ------------------------
 961      if(optdb==1) call printplot_matlu(matlu3,natom,paw_dmft%omega_lo(ifreq),"hybri_minus_asymp",74000,imre=1)
 962 
 963       ! Multiply (F-C_ij/i_wn) by (iw_n)**2 to find D_ij such that (F-C_ij/i_wn) -> D_ij/(iw_n)^2 only for last frequency.
 964       ! ------------------------------------------------------------------------------------------------------------------
 965      call copy_matlu(matlu3,matlu2,natom)
 966      call fac_matlu(matlu2,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2)
 967      if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"fminuscijtimesw2",75000,imre=1)
 968      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
 969        call copy_matlu(matlu2,matlu4,natom)
 970        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
 971 &       " == Print numerical (F(iw_n)-C_ij/iw_n)%iw_n^2 for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
 972        call wrtout(std_out,message,'COLL')
 973        call print_matlu(matlu4,natom,1)
 974      end if
 975 
 976       ! Built C_ij/iw_n+D_ij/(iw_n)^2
 977       ! ------------------------
 978      call copy_matlu(matlu4,matlu3,natom,opt_re=1)
 979      call fac_matlu(matlu3,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2)
 980      call add_matlu(matlu1,matlu3,matlu2,natom,1)
 981      if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"cij_w_plus_dij_w2",72700,imre=1)
 982       ! For test: put C_ij/i_wn +D_ij/(iw_n)^2 into weiss_for_rot
 983       ! --------------------------------------------
 984      !call copy_matlu(matlu2,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1)
 985 
 986 
 987    end do
 988 
 989    call destroy_matlu(matlu1,natom)
 990    call destroy_matlu(matlu2,natom)
 991    call destroy_matlu(matlu3,natom)
 992    call destroy_matlu(matlu4,natom)
 993    ABI_FREE(matlu1)
 994    ABI_FREE(matlu2)
 995    ABI_FREE(matlu3)
 996    ABI_FREE(matlu4)
 997  end if ! if opt_nondiag=1
 998 
 999 ! =========================================================================================
1000 ! Start big loop over atoms to compute hybridization and do the CTQMC
1001 ! =========================================================================================
1002  do iatom=1,cryst_struc%natom
1003    green%ecorr_qmc(iatom)=zero
1004    itypat=cryst_struc%typat(iatom)
1005    lpawu=paw_dmft%lpawu(iatom)
1006    tndim=2*lpawu+1
1007    if(lpawu/=-1) then
1008 
1009      nflavor=2*(tndim)
1010      if(testcode>=1) then
1011        nflavor=2
1012        if(testcode==2) then
1013          ispa=1
1014          ispb=2
1015          if(nspinor==1) ispb=1
1016          ima=1
1017          imb=1
1018          if(tndim>4) then
1019            ima=5 ! row
1020            imb=4 ! column
1021          end if
1022        end if
1023      end if
1024 
1025      ABI_MALLOC(fw1,(paw_dmft%dmft_nwlo,nflavor))
1026      ABI_MALLOC(fw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor))
1027      ABI_MALLOC(levels_ctqmc,(nflavor))
1028      ABI_MALLOC(levels_ctqmc_nd,(nflavor,nflavor))
1029      levels_ctqmc_nd=czero
1030      ABI_MALLOC(hybri_limit,(nflavor,nflavor))
1031      hybri_limit=czero
1032      fw1_nd=czero
1033      fw1=czero
1034 
1035      ! =================================================================
1036      ! Put hybridization in new arrays for CTQMC
1037      ! =================================================================
1038      if (testcode==0) then
1039        iflavor1=0
1040        iflavor2=0
1041 
1042        do isppol=1,nsppol
1043          do ispinor1=1,nspinor
1044            do ispinor2=1,nspinor
1045              do im1=1,tndim
1046                do im2=1,tndim
1047 
1048                  ! first diagonal terms whatever opt_nondiag
1049                  iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
1050                  iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
1051 
1052                  if ( iflavor1==iflavor2 ) then
1053 
1054                    ! Put weiss_for_rot in fw1
1055                    do ifreq=1,paw_dmft%dmft_nwlo
1056                      if(opt_fk==1) then
1057                        fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1)
1058                      else if (opt_fk==0) then
1059                        fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1)
1060                      end if
1061                    end do
1062                    fw1_nd(:,iflavor1,iflavor1)=fw1(:,iflavor1)
1063 
1064                    levels_ctqmc(iflavor1)=real(energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1),kind=dp)
1065                    hybri_limit(iflavor1,iflavor1)=hybri_coeff(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1)
1066 
1067 
1068                    ! case nsppol=nspinor=1
1069                    if(nsppol==1.and.nspinor==1) then
1070                      fw1(:,iflavor1+tndim)=fw1(:,iflavor1)
1071                      fw1_nd(:,iflavor1+tndim,iflavor1+tndim)=fw1(:,iflavor1)
1072                      levels_ctqmc(iflavor1+tndim)=levels_ctqmc(iflavor1)
1073                      hybri_limit(iflavor1+tndim,iflavor1+tndim)=hybri_limit(iflavor1,iflavor1)
1074                    end if
1075 
1076                  ! off diagonal terms
1077                  else
1078 
1079                    ! Put weiss_for_rot in fw1_nd
1080                    do ifreq=1,paw_dmft%dmft_nwlo
1081                      if(opt_fk==1) then
1082                        fw1_nd(ifreq,iflavor1,iflavor2)= &
1083 &                       weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1084                      else if (opt_fk==0) then
1085                        fw1_nd(ifreq,iflavor1,iflavor2)= &
1086 &                       weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1087                      end if
1088                    end do
1089                    hybri_limit(iflavor1,iflavor2)=hybri_coeff(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1090 
1091                    ! case nsppol=nspinor=1
1092                    if(nsppol==1.and.nspinor==1) then
1093                      fw1_nd(:,iflavor1+tndim,iflavor2+tndim) = fw1_nd(:,iflavor1,iflavor2)
1094                      hybri_limit(iflavor1+tndim,iflavor2+tndim)=hybri_limit(iflavor1,iflavor2)
1095                    end if
1096 
1097                  end if
1098 
1099 ! <  / HACK >
1100                end do !im2
1101              end do !im1
1102            end do  !ispinor2
1103          end do  !ispinor1
1104        end do  !isppol
1105 ! < HACK >
1106        ! JB. On 1000 cpus this can not work since all CPU try to open/write the files
1107        ! Action : Don't print it or check only one cpu does it.
1108 
1109        if(pawprtvol>=10000000) then
1110          write(message,'(a,2x,a)') ch10,  " == Hybri for all flavors for CTQMC "
1111          call wrtout(std_out,message,'COLL')
1112          do iflavor1=1,nflavor
1113            write(message,'(4x,14(2e14.5,2x))') (hybri_limit(iflavor1,iflavor2),iflavor2=1,nflavor)
1114            call wrtout(std_out,message,'COLL')
1115          end do
1116 
1117          if (open_file('Hybri_cijoveromega',message, newunit=unt, status='unknown', form='formatted') /= 0) then
1118            ABI_ERROR(message)
1119          end if
1120          if (open_file('Hybri',message,newunit=unt2,status='unknown',form='formatted') /= 0) then
1121            ABI_ERROR(message)
1122          end if
1123          do ifreq=1,paw_dmft%dmft_nwlo
1124            !  weiss_for_rot is G_0^-1-iw_n=-(F-levels)
1125            if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"weissbefore112",30000)
1126          end do
1127          do iflavor1=1,nflavor
1128            do iflavor2=1,nflavor
1129              do ifreq=1,paw_dmft%dmft_nwlo
1130                omega=pi*paw_dmft%temp*(two*float(ifreq)-1)
1131                ! fw1_nd is -G_0^+1-iw_n=(F-levels)
1132                write(unt,'(300e16.5)') paw_dmft%omega_lo(ifreq)&
1133 &               ,fw1_nd(ifreq,iflavor1,iflavor2)-hybri_limit(iflavor1,iflavor2)/cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)
1134                write(unt2,'(300e16.5)') paw_dmft%omega_lo(ifreq),fw1_nd(ifreq,iflavor1,iflavor2)
1135              end do
1136              write(unt,*)
1137              write(unt2,*)
1138            end do
1139          end do
1140          close(unt)
1141          close(unt2)
1142        end if
1143      end if ! testcode
1144    ! </ HACK >
1145 
1146      ! ====================================================================================
1147      !  TEST
1148      !  For testing purpose, built ultra simple hybridization (constant in
1149      !  imaginary time or very simple) or extract some part of the calculated hybridization
1150      ! ====================================================================================
1151      if(testcode>=1) then
1152        dmft_nwlo=paw_dmft%dmft_nwlo
1153        paw_dmft%dmft_nwlo=paw_dmft%dmftqmc_l
1154        ABI_MALLOC(gw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor))
1155        gw1_nd=czero
1156 
1157        !  Call testcode_ctqmc: built simple hybridization
1158        !--------------------------------------------------
1159        if (testcode==1) then
1160          call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,&
1161 &         levels_ctqmc,hybri_limit,nflavor,1,paw_dmft%temp,testrot,testcode,umod)
1162        !  Select 2x2 hybridization matrix from the current larger matrix
1163        !  ima and imb are defined above.
1164        !----------------------------------------------------------------
1165        else if (testcode==2) then
1166          !close(unt)
1167          !close(unt2)
1168          call testcode_ctqmc_b(energy_level,hybri_coeff,weiss_for_rot,paw_dmft%dmftqmc_l,fw1_nd,&
1169 &         levels_ctqmc,levels_ctqmc_nd,hybri_limit,paw_dmft%temp,umod,opt_diag,opt_fk)
1170        end if
1171 
1172        ! Calculation of Inverse Green's function from hybridization
1173        !-------------------------------------------------------------
1174        do if1=1,2
1175          do if2=1,2
1176            do ifreq=1,paw_dmft%dmftqmc_l
1177              omega=pi*paw_dmft%temp*(two*float(ifreq)-1)
1178              if(if1==if2) then
1179                gw1_nd(ifreq,if1,if2) =  (cmplx(0.d0,omega,kind=dp)-fw1_nd(ifreq,if1,if2))
1180              else
1181                gw1_nd(ifreq,if1,if2) =  (-fw1_nd(ifreq,if1,if2))
1182              end if
1183            end do
1184          end do
1185        end do
1186        ! Calculation of Green's function (call inverse)
1187        !-------------------------------------------------------------
1188        do ifreq=1,paw_dmft%dmftqmc_l
1189          call xginv(gw1_nd(ifreq,:,:),2)
1190        end do
1191        write(std_out,*) " testctqmc high frequency limit of hybridization",fw1_nd(paw_dmft%dmftqmc_l,:,:)
1192 
1193        ! Integrate Green's function
1194        !-------------------------------------------------------------
1195        do if1=1,2
1196          do if2=1,2
1197            call int_fct(gw1_nd(:,if1,if2),(if1==if2),2,paw_dmft,integral(if1,if2))  ! test_1
1198          end do
1199        end do
1200        ! Write Occupations
1201        write(std_out,*) "Occupation of model in matrix form"
1202        do if1=1,2
1203          write(std_out,'(2(2f13.5,3x))') ((integral(if1,if2)+conjg(integral(if2,if1)))/two,if2=1,2)
1204        end do
1205        write(std_out,*) "Limit of hybridization "
1206        do if1=1,2
1207          write(std_out,'(2(2f13.5,3x))') (hybri_limit(if1,if2),if2=1,2)
1208        end do
1209 
1210        ! If opt_fk=0, give Green's function to CTQMC code instead of
1211        ! hybridization
1212        !-------------------------------------------------------------
1213        if(opt_fk==0) then
1214          fw1_nd=gw1_nd
1215        end if
1216 
1217        ABI_FREE(gw1_nd)
1218        paw_dmft%dmft_nwlo=dmft_nwlo
1219 
1220      ! and testcode>1
1221      end if
1222 
1223 
1224      call flush_unit(std_out)
1225 ! =================================================================
1226 
1227 ! ___________________________________________________________________________________
1228 !
1229 !  THIRD PART : CALL CTQMC
1230 ! ___________________________________________________________________________________
1231 
1232 ! =================================================================
1233 !    Main calls to CTQMC code in ABINIT (INITIALIZATION and OPTIONS)
1234 ! =================================================================
1235      if(paw_dmft%dmft_solv==5.or.paw_dmft%dmft_solv==8) then
1236 
1237        write(message,'(a,2x,a)') ch10,&
1238 &       " == Initializing CTQMC"
1239        call wrtout(std_out,message,'COLL')
1240 
1241 !    Initialisation
1242 ! =================================================================
1243      if(paw_dmft%dmft_solv==5) then
1244        nomega=paw_dmft%dmftqmc_l
1245        call CtqmcInterface_init(hybrid,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, &
1246 &       paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,&
1247 &       std_out,paw_dmft%spacecomm)
1248 !    options
1249 ! =================================================================
1250        call CtqmcInterface_setOpts(hybrid,&
1251        opt_Fk      =opt_fk,&
1252 &       opt_order   =paw_dmft%dmftctqmc_order ,&
1253 &       opt_histo   =paw_dmft%dmftctqmc_config ,&
1254 &       opt_movie   =paw_dmft%dmftctqmc_mov   ,&
1255 &       opt_analysis=paw_dmft%dmftctqmc_correl,&
1256 &       opt_check   =paw_dmft%dmftctqmc_check ,&
1257 &       opt_noise   =paw_dmft%dmftctqmc_grnns ,&
1258 &       opt_spectra =paw_dmft%dmftctqmc_mrka  ,&
1259 &       opt_gmove   =paw_dmft%dmftctqmc_gmove )
1260      endif
1261 
1262      if(paw_dmft%dmft_solv==8) then
1263        nomega=paw_dmft%dmftqmc_l
1264        call CtqmcoffdiagInterface_init(hybridoffdiag,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, &
1265 &        paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,&
1266 &        paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,&
1267 &        std_out,paw_dmft%spacecomm,opt_nondiag)
1268 !    options
1269 ! =================================================================
1270        call CtqmcoffdiagInterface_setOpts(hybridoffdiag,&
1271        opt_Fk      =opt_fk,&
1272 &       opt_order   =paw_dmft%dmftctqmc_order ,&
1273 &       opt_histo   =paw_dmft%dmftctqmc_config ,&
1274 &       opt_movie   =paw_dmft%dmftctqmc_mov   ,&
1275 &       opt_analysis=paw_dmft%dmftctqmc_correl,&
1276 &       opt_check   =paw_dmft%dmftctqmc_check ,&
1277 &       opt_noise   =paw_dmft%dmftctqmc_grnns ,&
1278 &       opt_spectra =paw_dmft%dmftctqmc_mrka  ,&
1279 &       opt_gmove   =paw_dmft%dmftctqmc_gmove )
1280      endif
1281 
1282        write(message,'(a,2x,2a)') ch10, " == Initialization CTQMC done", ch10
1283        call wrtout(std_out,message,'COLL')
1284      end if
1285 
1286      if(paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7.or.paw_dmft%dmft_solv==9) then
1287        ABI_MALLOC(gw_tmp_nd,(paw_dmft%dmft_nwli,nflavor,nflavor))
1288        !because size allocation problem with TRIQS paw_dmft%dmft_nwlo must be >= paw_dmft%dmft_nwli
1289          open(unit=505,file=trim(paw_dmft%filapp)//"_Legendre_coefficients.dat", status='unknown',form='formatted')
1290      else
1291        if(paw_dmft%dmft_solv==5) then
1292          ABI_MALLOC(gw_tmp,(paw_dmft%dmft_nwlo,nflavor+1))
1293        end if
1294        ABI_MALLOC(gw_tmp_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor+1))
1295        !use  gw_tmp to put freq
1296        do ifreq=1,paw_dmft%dmft_nwlo
1297          if(paw_dmft%dmft_solv==5) gw_tmp(ifreq,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
1298          gw_tmp_nd(ifreq,nflavor,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
1299        end do
1300      end if
1301 
1302      ABI_MALLOC(gtmp,(paw_dmft%dmftqmc_l,nflavor))
1303      ! THIS IS A BACKUP PLAN. USING paw_dmft%hybrid makes a segfault on TIKAL
1304      ! PSC with MPI only (and max2_open64). paw_dmf%hybrid is corrupted
1305      ! somewhere but I could not find the place in all DMFT routines
1306      ABI_MALLOC(gtmp_nd,(paw_dmft%dmftqmc_l,nflavor,nflavor))
1307      call flush_unit(std_out)
1308 
1309      ! =================================================================
1310      !    BEGIN CALL TO CTQMC SOLVERS
1311      ! =================================================================
1312 
1313      if(testcode==0) then
1314 
1315        ! =================================================================
1316        !    CTQMC run Abinit
1317        ! =================================================================
1318        if(paw_dmft%dmft_solv==5) then
1319 
1320          ABI_MALLOC(docc,(1:nflavor,1:nflavor))
1321          docc(:,:) = zero
1322          call CtqmcInterface_run(hybrid,fw1(1:paw_dmft%dmftqmc_l,:),Gtau=gtmp,&
1323 &         Gw=gw_tmp,D=docc(:,:),E=green%ecorr_qmc(iatom),&
1324 !&       matU=hu(itypat)%udens,opt_levels=levels_ctqmc)
1325 &         matU=udens_atoms(iatom)%value,opt_levels=levels_ctqmc)
1326          call data4entropyDMFT_setDocc(paw_dmft%forentropyDMFT,iatom,docc)
1327          ABI_FREE(docc)
1328          !DO iflavor = 1, nflavor
1329          !  hybrid%Hybrid%Greens(iflavor)%oper(1:this%samples) = gtmp(1:this%samples,iflavor)
1330          !  CALL GreenHyb_forFourier(this%Greens(iflavor), Gomega=Gw(:,iflavor), omega=Gw(:,this%flavors+1))
1331          !END DO
1332 
1333        ! =================================================================
1334        !    CTQMC run Abinit off diagonal terms in hybridization
1335        ! =================================================================
1336        else if (paw_dmft%dmft_solv==8) then
1337        ! =================================================================
1338 
1339          ABI_MALLOC(docc,(1:nflavor,1:nflavor))
1340          docc(:,:) = zero
1341          call CtqmcoffdiagInterface_run(hybridoffdiag,fw1_nd(1:paw_dmft%dmftqmc_l,:,:),Gtau=gtmp_nd,&
1342 &        Gw=gw_tmp_nd,D=doccsum,E=green%ecorr_qmc(iatom),&
1343 &        Noise=noise,matU=udens_atoms(iatom)%value,Docc=docc,opt_levels=levels_ctqmc,hybri_limit=hybri_limit)
1344          ! For entropy (alternative formulation)
1345          if(paw_dmft%ientropy==1) then
1346            EE=zero
1347            do if1=1,nflavor
1348              do if2=if1+1,nflavor
1349                EE=EE+docc(if1,if2)*udens_atoms_for_s(iatom)%value(if1,if2)
1350          !      write(std_out,*) udens_atoms_for_s(iatom)%value(if1,if2),docc(if1,if2)
1351              enddo
1352            enddo
1353            ! Here in udens U=1, J=J/U, so we need to multiply bu U/Ha_eV
1354            write(message,'(a,3(f14.10,3x))') "For entropy calculation E_corr_qmc, u_for_s, j_for,s", &
1355 &          paw_dmft%u_for_s*EE/Ha_eV,paw_dmft%u_for_s,paw_dmft%j_for_s
1356            call wrtout(std_out,message,'COLL')
1357            EE=zero
1358            do if1=1,nflavor
1359              do if2=if1+1,nflavor
1360                EE=EE+docc(if1,if2)*udens_atoms(iatom)%value(if1,if2)
1361          !      write(std_out,*) udens_atoms(iatom)%value(if1,if2),docc(if1,if2)
1362              enddo
1363            enddo
1364            ! Here in udens U=U, J=J, so we obtain directly the results
1365            write(message,'(a,3(f14.10,3x))') "Reference   calculation E_corr_qmc, upawu  , jpawu  ", &
1366 &             EE,hu(itypat)%upawu*Ha_eV,hu(itypat)%jpawu*Ha_eV
1367            call wrtout(std_out,message,'COLL')
1368          endif
1369          ABI_FREE(docc)
1370        ! TODO: Handle de luj0 case for entropy
1371 
1372        ! =================================================================
1373        !    CTQMC run TRIQS
1374        ! =================================================================
1375        else if ((paw_dmft%dmft_solv>=6.and.paw_dmft%dmft_solv<=7).or.paw_dmft%dmft_solv==9) then
1376        ! =================================================================
1377 
1378          if ( paw_dmft%dmft_solv==9.or.(paw_dmft%dmftqmc_l >= (2*paw_dmft%dmft_nwli)+1) ) then
1379 
1380            call ctqmc_calltriqs(paw_dmft,cryst_struc,hu,levels_ctqmc,gtmp_nd,gw_tmp_nd,fw1_nd,leg_measure,iatom)
1381 
1382          else
1383            write(message,'(2a)') ch10," Can't launch TRIQS CTHYB solver because dmftqmc_l must be >= 2*dmft_nwli + 1"
1384            ABI_ERROR(message)
1385          end if
1386 
1387        end if
1388 
1389      ! =================================================================
1390      !    CTQMC run for tests
1391      ! =================================================================
1392      else if (testcode>=1) then
1393        call CtqmcInterface_run(hybrid,fw1(1:nomega,:),Gtau=gtmp,&
1394 &       Gw=gw_tmp,E=green%ecorr_qmc(iatom),&
1395 &       matU=umod,opt_levels=levels_ctqmc)
1396 
1397       ! for non diagonal code
1398       !       call CtqmcInterface_run(hybrid,fw1_nd(1:nomega,:,:),Gtau=gtmp_nd,&
1399       !&       Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),&
1400       !&       Noise=Noise,matU=umod,opt_levels=levels_ctqmc,hybri_limit=hybri_limit)
1401 
1402       !  If test of the code is activated, and testrot =1 rotate back green's function   and stop the code.
1403       ! --------------------------------------------------------------------------------------------------
1404        if(testcode==1) then
1405 
1406          call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,&
1407 &         levels_ctqmc,hybri_limit,nflavor,2,paw_dmft%temp,testrot,testcode,umod)
1408 
1409          write(message,'(2a)') ch10,' testcode end of test calculation'
1410          ABI_ERROR(message)
1411        end if
1412        if(testcode==2) then
1413          write(message,'(2a)') ch10,' testcode 2 end of test calculation'
1414          ABI_ERROR(message)
1415        end if
1416 
1417      end if
1418      ! =================================================================
1419      !    END CALL TO CTQMC SOLVERS
1420      ! =================================================================
1421 
1422 
1423      ! Print green function is files directly from CTQMC
1424      ! --------------------------------------------------
1425      call ctqmcoutput_printgreen(paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom)
1426 
1427 
1428      ! If the CTQMC code in ABINIT was used, then destroy it and deallocate arrays
1429      ! ----------------------------------------------------------------------------
1430      if(paw_dmft%dmft_solv<6.and.paw_dmft%dmft_solv>7) then
1431      !Nothing just hybrid var problem
1432      else
1433        write(message,'(a,2x,a)') ch10,&
1434 &       " == Destroy CTQMC"
1435        call wrtout(std_out,message,'COLL')
1436        if(paw_dmft%dmft_solv==5) call CtqmcInterface_finalize(hybrid)
1437        if(paw_dmft%dmft_solv==8) call CtqmcoffdiagInterface_finalize(hybridoffdiag)
1438        write(message,'(a,2x,a)') ch10,&
1439 &       " == Destroy CTQMC done"
1440        call wrtout(std_out,message,'COLL')
1441      end if
1442      ABI_FREE(hybri_limit)
1443      ABI_FREE(levels_ctqmc_nd)
1444      ABI_FREE(levels_ctqmc)
1445      ABI_FREE(fw1)
1446      ABI_FREE(fw1_nd)
1447 
1448 ! ___________________________________________________________________________________
1449 !
1450 !  FOURTH PART : USE OUTPUT OF CTQMC AND THEN DO BACK ROTATION
1451 ! ___________________________________________________________________________________
1452 !
1453 
1454      ! Put green's function values from CTQMC into green structure
1455      !-------------------------------------------------------------
1456      call ctqmcoutput_to_green(green,paw_dmft,gtmp_nd,gw_tmp_nd,gtmp,gw_tmp,iatom,leg_measure,opt_nondiag)
1457 
1458      ! Deallocate arrays for CTQMC
1459      !-----------------------------
1460      if(paw_dmft%dmft_solv<6) then
1461        ABI_FREE(gw_tmp)
1462      endif
1463      ABI_FREE(gw_tmp_nd)
1464      ABI_FREE(gtmp)
1465      ABI_FREE(gtmp_nd)
1466 
1467 
1468 
1469      ! Do Fourier transform if it was not done (ie if TRIQS is used without legendre measurement)
1470      !----------------------------------------------------------------------------------------------
1471      if(opt_nondiag==1) then  ! (As leg_measure is activated by defautl, this fourier is never done).
1472        if(paw_dmft%dmft_solv>=6.and..not.leg_measure.and.paw_dmft%dmft_solv<=7) then
1473          write(message,'(2a,i3,13x,a)') ch10,'   ===  Direct Fourier Transform t->w of Weiss Field'
1474          call wrtout(std_out,message,'COLL')
1475          call fourier_green(cryst_struc,green,paw_dmft,&
1476 &         pawang,opt_ksloc=2,opt_tw=1)
1477        end if
1478      endif
1479 
1480 
1481    end if
1482 
1483  end do ! iatom
1484 ! =========================================================================================
1485 !  End big loop over atoms to compute hybridization and do the CTQMC
1486 ! =========================================================================================
1487 
1488  if(paw_dmft%dmft_prgn==1) then
1489    call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=2)
1490    call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
1491  end if
1492  !write(message,'(i3,4x,2e21.14)') 6,weiss_for_rot%oper(1)%matlu(1)%mat(1,1,1,1,1)
1493  !call wrtout(std_out,message,'COLL')  ! debug
1494 ! =================================================================
1495 ! Inverse Weiss, then
1496 ! Copy Weiss_for_rot into weiss and rotate back weiss to the original basis
1497 ! =================================================================
1498 
1499 ! ABI_MALLOC(shift,(natom))
1500 ! do ifreq=1,paw_dmft%dmft_nwlo
1501 !  ! First weiss_for_rot contains -G_0^-1+iw_n
1502 !  ! -------------------------------------------
1503 !  ! Compute G_0^-1-iw_n
1504 !  ! --------------------
1505 !       write(6,*) "1"
1506 !  if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone)
1507 !
1508 !
1509 !       write(6,*) "2"
1510 !  ! Compute G_0^-1
1511 !  ! --------------------
1512 !  shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
1513 !  if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift,signe=1)
1514 !
1515 !       write(6,*) "3"
1516 !  ! Compute G_0
1517 !  ! --------------------
1518 !   call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1)
1519 !   ! No need to copy if weiss_for_rot is a pointer to weiss ...
1520 !!   if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0)
1521 !!   if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
1522 !
1523 !  ! Compute G_0 in the original basis
1524 !  ! --------------------
1525 !   call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
1526 ! end do
1527 ! ABI_FREE(shift)
1528 
1529 ! =================================================================
1530 ! Here compute Self energy from Dyson and print it
1531 ! Warning : Weiss_for_rot is inversed inside dyson
1532 ! =================================================================
1533 ! call initialize_self(self,paw_dmft)
1534 ! call dyson(green,paw_dmft,self,weiss_for_rot,opt_weissself=2)
1535 ! call rw_self(self,mpi_enreg,paw_dmft,prtopt=2,opt_rw=2,opt_char="diag")
1536 ! call destroy_self(self)
1537  !write(message,'(i3,4x,2e21.14)') 7,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1)
1538  !call wrtout(std_out,message,'COLL')  ! debug
1539 
1540 ! =================================================================
1541 ! Rotate back green function to original basis (non-diagonal)
1542 !  (and Weiss for further use: might be useful if an back Fourier
1543 !     transformation is done).
1544 ! =================================================================
1545  if(pawprtvol>=3) then
1546    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
1547 &  " == Print green function for small tau after CTQMC"  ! debug
1548    call wrtout(std_out,message,'COLL')  ! debug
1549    call print_matlu(green%oper_tau(1)%matlu,natom,1)  ! debug
1550    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
1551 &  " == Print green function for small freq after CTQMC"  ! debug
1552    call wrtout(std_out,message,'COLL')  ! debug
1553    call print_matlu(green%oper(1)%matlu,natom,1)  ! debug
1554  end if
1555 
1556 
1557 !  === Compute rotated Occupations in green%occup_tau
1558  call occup_green_tau(green)
1559 
1560  if(pawprtvol>=3) then
1561 !  === Compute non rotated Occupations in green%occup_tau
1562    write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the ctqmc basis"
1563    call wrtout(std_out,message,'COLL')
1564    call print_matlu(green%occup_tau%matlu,natom,1)
1565  end if
1566 
1567  write(message,'(a,2x,a,f13.5)') ch10,&
1568 & " == Rotate Green function to original basis "
1569  call wrtout(std_out,message,'COLL')
1570  !write(message,'(i3,4x,2e21.14)') 8,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1)
1571  !call wrtout(std_out,message,'COLL')  ! debug
1572 
1573  ! Rotate oper_tau into Ylm basis and then Slm basis
1574  !-------------------------------------------------------------
1575  do itau=1,paw_dmft%dmftqmc_l
1576    if(opt_diag/=0) call rotate_matlu(green%oper_tau(itau)%matlu,eigvectmatlu,natom,3,0)
1577    if(useylm==1) call slm2ylm_matlu(green%oper_tau(itau)%matlu,natom,2,0)
1578  end do
1579 
1580  ! Rotate occup_tau into Ylm basis and then Slm basis
1581 
1582  ! Rotate occup_tau into Ylm basis and then Slm basis
1583  !-------------------------------------------------------------
1584  if(opt_diag/=0) call rotate_matlu(green%occup_tau%matlu,eigvectmatlu,natom,3,0)
1585  if(useylm==1) then
1586    write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Ylm basis"
1587    call wrtout(std_out,message,'COLL')
1588    call print_matlu(green%occup_tau%matlu,natom,1)
1589  end if
1590 
1591  if(useylm==1) call slm2ylm_matlu(green%occup_tau%matlu,natom,2,0)
1592  write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Slm basis"
1593  call wrtout(std_out,message,'COLL')
1594  call print_matlu(green%occup_tau%matlu,natom,1)
1595 
1596  ! Put Weiss off diagonal terms to zero because Green function will not have any offdiag terms
1597  !------------------------------------------------------------------------------
1598  !   (if opt_nondiag=0 ie dmft_solv=5)
1599  do ifreq=1,paw_dmft%dmft_nwlo
1600    if(opt_nondiag==0) call zero_matlu(weiss%oper(ifreq)%matlu,natom,onlynondiag=1)
1601  end do
1602  !    ( if opt_nondiag=0, then:
1603  !       As Green's function is diagonal, one suppress off diag  terms in Weiss, if any.
1604  !      (If off diag are non zero in the density matrix and thus in the Green's function,
1605  !       there is a warning in checkreal_matlu above).)
1606 
1607  ! Rotate Green's and Weiss functions into Ylm basis and then Slm basis
1608  !-------------------------------------------------------------
1609  do ifreq=1,paw_dmft%dmft_nwlo
1610    if(opt_diag/=0) call rotate_matlu(green%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
1611    if(useylm==1) call slm2ylm_matlu(green%oper(ifreq)%matlu,natom,2,0)
1612    if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
1613    if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0)
1614  end do
1615  !write(message,'(i3,4x,2e21.14)') 10,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1)
1616  !call wrtout(std_out,message,'COLL')  ! debug
1617 
1618  if(pawprtvol>=3) then
1619    write(message,'(a,2x,a,f13.5)') ch10,&                  ! debug
1620 &  " == Print green function for small time after rotation (in the original basis)" ! debug
1621    call wrtout(std_out,message,'COLL')                  ! debug
1622    call print_matlu(green%oper_tau(1)%matlu,natom,1)  ! debug
1623    write(message,'(a,2x,a,f13.5)') ch10,&                  ! debug
1624 &  " == Print green function for small freq after rotation (in the original basis)" ! debug
1625    call wrtout(std_out,message,'COLL')                  ! debug
1626    call print_matlu(green%oper(1)%matlu,natom,1)  ! debug
1627    !< HACK >
1628    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
1629 &  " == Print diagonalized weiss_for_rot function after rotation for small freq in the ctqmc basis"  ! debug
1630    call wrtout(std_out,message,'COLL')  ! debug
1631    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1)  ! debug
1632    !</ HACK >
1633    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
1634 &  " == Print weiss function for small freq in the original basis"  ! debug
1635    call wrtout(std_out,message,'COLL')  ! debug
1636    call print_matlu(weiss%oper(1)%matlu,natom,1)  ! debug
1637 
1638    do ifreq=1,paw_dmft%dmft_nwlo
1639      call sym_matlu(cryst_struc,weiss%oper(ifreq)%matlu,pawang,paw_dmft)
1640    end do
1641    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
1642 &  " == Print symetrized weiss function for small freq in the original basis"  ! debug
1643    call wrtout(std_out,message,'COLL')  ! debug
1644    call print_matlu(weiss%oper(1)%matlu,natom,1)  ! debug
1645  end if
1646 
1647 
1648  ABI_MALLOC(matlu1,(natom))
1649  call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
1650  call copy_matlu(green%occup_tau%matlu,matlu1,natom)
1651  call sym_matlu(cryst_struc,matlu1,pawang,paw_dmft)
1652 
1653  write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the original basis"
1654  call wrtout(std_out,message,'COLL')
1655  call print_matlu(green%occup_tau%matlu,natom,1)
1656 
1657  write(message,'(a,2x,a,f13.5)') ch10," == Symetrized occupations"
1658  call wrtout(std_out,message,'COLL')
1659  call print_matlu(matlu1,natom,1)
1660 
1661  call diff_matlu("CTQMC Occup","CTQMC Occup symetrized",green%occup_tau%matlu,matlu1,natom,0,tol4,ierr)
1662  call destroy_matlu(matlu1,natom)
1663  ABI_FREE(matlu1)
1664 
1665 ! =================================================================
1666 ! Symetrise green function G(tau) and G(ifreq) to recover symetry
1667 ! artificially broken by QMC
1668 ! =================================================================
1669  write(message,'(a,2x,a,f13.5)') ch10,&
1670 & " == Symetrise green function after QMC "
1671  call wrtout(std_out,message,'COLL')
1672  do itau=1,paw_dmft%dmftqmc_l
1673    call sym_matlu(cryst_struc,green%oper_tau(itau)%matlu,pawang,paw_dmft)
1674  end do
1675  do ifreq=1,paw_dmft%dmft_nwlo
1676    call sym_matlu(cryst_struc,green%oper(ifreq)%matlu,pawang,paw_dmft)
1677  end do
1678  if(pawprtvol>=3) then
1679    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
1680 &  " == Print green function for small time after symetrisation"  !  debug
1681    call wrtout(std_out,message,'COLL')  ! debug
1682    call print_matlu(green%oper_tau(1)%matlu,natom,1)  ! debug
1683    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
1684 &  " == Print green function for small freq after symetrisation"  !  debug
1685    call wrtout(std_out,message,'COLL')  ! debug
1686    call print_matlu(green%oper(1)%matlu,natom,1)  ! debug
1687  end if
1688  if(paw_dmft%dmft_prgn==1) then
1689    call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=2)
1690    call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
1691  end if
1692 
1693 !  === Compute Occupations  (Symetrized from oper_tau)
1694  call occup_green_tau(green)
1695 
1696 
1697 !  === Print occupations
1698  call printocc_green(green,6,paw_dmft,3)
1699 
1700  call destroy_oper(energy_level)
1701  call destroy_matlu(dmat_diag,natom)
1702  ABI_FREE(dmat_diag)
1703  call destroy_matlu(identity,natom)
1704  ABI_FREE(identity)
1705  do iatom=1,cryst_struc%natom
1706    lpawu=paw_dmft%lpawu(iatom)
1707    if(lpawu/=-1) then
1708      do isppol=1,nsppol
1709        ABI_FREE(eigvectmatlu(iatom,isppol)%value)
1710        !ABI_FREE(udens_atoms(iatom))
1711      end do
1712      ABI_FREE(udens_atoms(iatom)%value)
1713    end if
1714  end do
1715  ABI_FREE(udens_atoms)
1716  ABI_FREE(eigvectmatlu)
1717  call destroy_green(weiss_for_rot)
1718 ! call destroy_green(gw_loc)
1719 ! call destroy_green(greendft)
1720 
1721 !  destroy limit of hybridization
1722  call destroy_matlu(hybri_coeff,paw_dmft%natom)
1723  ABI_FREE(hybri_coeff)
1724 
1725 end subroutine qmc_prep_ctqmc

m_forctqmc/testcode_ctqmc [ Functions ]

[ Top ] [ m_forctqmc ] [ Functions ]

NAME

 testcode_ctqmc

FUNCTION

 Setup ultra simple hybridization to test CTQMC in simple situations.

INPUTS

 gtmp_nd
 gw_tmp_nd
 temp = temperature
 dmftqmc_l = number of times slices
 nflavor = number of flavor
 testrot = 0/1 if rotation of hybridization is tested or not
 testcode = 1 if tests are activated.
 opt = 1/2 if pre or postprocessing of CTQMC data.

OUTPUT

 fw1_nd = non diagonal hybridization
 fw1 = hybridization
 umod = value of U

SIDE EFFECTS

  gtmp_nd
  gw_tmp_nd

NOTES

SOURCE

1856 subroutine testcode_ctqmc(dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,levels_ctqmc,hybri_limit,&
1857 &   nflavor,opt,temp,testrot,testcode,umod)
1858 
1859 
1860 !Arguments ------------------------------------
1861 !scalars
1862  integer, intent(in) :: dmftqmc_l,nflavor,testrot,testcode,opt
1863  real(dp), intent(in) :: temp
1864  real(dp), intent(out) :: umod(2,2)
1865  complex(dpc), intent(inout) :: gw_tmp_nd(:,:,:)
1866  real(dp),  intent(inout) :: gtmp_nd(:,:,:)
1867  complex(dpc), intent(out) :: fw1(:,:)
1868  complex(dpc), intent(out) :: fw1_nd(:,:,:)
1869  real(dp),  intent(inout) :: levels_ctqmc(:)
1870  complex(dpc),  intent(inout) :: hybri_limit(:,:)
1871 
1872 !Local variables ------------------------------
1873  character(len=500) :: message
1874  integer :: ifreq, itau,realrot,simplehyb
1875  real(dp) :: omega
1876  real(dp) :: tbi1,tbi2,e2,tbi3,tbi4,e3,e4,tbi21,tbi12,e3b,e4b,tbi21b,tbi12b
1877  complex(dpc) :: e1
1878 ! arrays
1879  complex(dpc) :: RR(2,2)
1880  complex(dpc) :: RR1(2,2)
1881  complex(dpc) :: RRi(2,2)
1882  complex(dpc) :: RRt(2,2)
1883 ! ************************************************************************
1884  if (testcode==0) return
1885  if (nflavor/=2) then
1886    write(message,'(2a)') ch10,' testcode nflavor.ne.2'
1887    ABI_ERROR(message)
1888  end if
1889 
1890  simplehyb=2
1891  simplehyb=1
1892  simplehyb=3
1893  !=========================
1894  ! Built rotation matrix
1895  !=========================
1896  realrot=0
1897  realrot=2
1898  if (realrot==1) then
1899    ! Real rotation
1900    !=========================
1901    RR(1,1)  =  SQRT(3.d0)/2.d0
1902    RR(1,2)  = -1.d0/2.d0
1903    RR(2,1)  =  1.d0/2.d0
1904    RR(2,2)  =  SQRT(3.d0)/2.d0
1905  else if (realrot==2) then
1906    ! Real rotation
1907    !=========================
1908    RR(1,1)  =  SQRT(1.d0/2.d0)
1909    RR(1,2)  = -SQRT(1.d0/2.d0)
1910    RR(2,1)  =  SQRT(1.d0/2.d0)
1911    RR(2,2)  =  SQRT(1.d0/2.d0)
1912  else
1913    ! Complex rotation
1914    !=========================
1915    RR(1,1)  =  CMPLX(one,two)
1916    RR(1,2)  =  CMPLX(one,one)
1917    RR(2,1)  =  CMPLX(one,-one)
1918    RR(2,2)  =  CMPLX(-one,two)
1919    RR=RR/sqrt(seven)
1920  end if
1921  ! Check rotation is unitary
1922  !==========================
1923  RRi(1,1) =  conjg(RR(1,1))
1924  RRi(1,2) =  conjg(RR(2,1))
1925  RRi(2,1) =  conjg(RR(1,2))
1926  RRi(2,2) =  conjg(RR(2,2))
1927  RR1(:,:)  = MATMUL ( RR(:,:) , RRi(:,:)          )
1928  !write(6,*) "RR1",RR1
1929  if(abs(RR1(1,1)-one).gt.tol7.or.abs(RR1(1,2)).gt.tol7.or.abs(RR1(2,2)-one).gt.tol7.or.abs(RR1(2,1)).gt.tol7) then
1930    write(message,'(2a)') ch10,' testcode error in rotation matrix'
1931    ABI_ERROR(message)
1932  end if
1933 
1934 
1935  !=================================
1936  ! Built hybridization  for CTQMC
1937  !=================================
1938  if (opt==1) then
1939 
1940  !  Parameters: tight-binding + U
1941  !  firt test of the code try umod=0, and (tbi1,tbi2,e1,e2)=(2,1,0.5,0.0) testrot=1
1942  !  second test of the code try umod=four, and (tbi1,tbi2,e1,e2)=(2,1,0.0,0.0) testrot=1
1943  !=======================================================================================
1944    fw1_nd(:,:,:)= czero
1945    tbi1=2.0_dp
1946    tbi2=1.0_dp
1947    tbi3=1.0_dp
1948    tbi4=1.0_dp
1949    tbi12=2.5_dp
1950    tbi12b=2.5_dp
1951    tbi21=2.5_dp
1952    tbi21b=2.5_dp
1953    e1=cmplx(0.0,0.0,8)
1954    e2=zero
1955    e3=0.2
1956    e4=0.3
1957    e3b=0.3
1958    e4b=-0.2
1959    umod(:,:)=0.d0
1960 
1961    if(testrot==1.and.(abs(tbi1-tbi2)<tol6)) then
1962      write(message,'(3a)') ch10,' testrot=1 with tbi1=tbi2 is equivalent' &
1963      ,'to testrot=0: change testrot'
1964      ABI_WARNING(message)
1965    end if
1966    ! Built fw1_nd
1967    !==============
1968    do ifreq=1,dmftqmc_l
1969 
1970      omega=pi*temp*(two*float(ifreq)-1)
1971 
1972      if(simplehyb==1) then
1973        fw1_nd(ifreq,1,1) =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
1974        fw1_nd(ifreq,2,2) =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
1975        fw1(ifreq,1)      =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
1976        fw1(ifreq,2)      =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
1977        hybri_limit(1,1)=tbi1**2
1978        hybri_limit(2,2)=tbi2**2
1979        hybri_limit(1,2)=0.d0
1980        hybri_limit(2,1)=0.d0
1981      else if(simplehyb==2) then
1982        fw1_nd(ifreq,1,1) =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)+tbi3**2/(dcmplx(0.d0,omega)-e3)
1983        fw1_nd(ifreq,2,2) =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)+tbi4**2/(dcmplx(0.d0,omega)-e4)
1984        fw1(ifreq,1)      =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
1985        fw1(ifreq,2)      =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
1986      else if(simplehyb==3) then
1987        fw1_nd(ifreq,1,1) =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
1988        fw1_nd(ifreq,2,2) =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
1989        fw1_nd(ifreq,1,2) =  tbi12**2/(dcmplx(0.d0,omega)-e3)+tbi12b**2/(dcmplx(0.d0,omega)-e3b)
1990        fw1_nd(ifreq,2,1) =  tbi21**2/(dcmplx(0.d0,omega)-e4)+tbi21b**2/(dcmplx(0.d0,omega)-e4b)
1991        fw1(ifreq,1)      =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
1992        fw1(ifreq,2)      =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
1993        hybri_limit(1,1)=tbi1**2
1994        hybri_limit(2,2)=tbi2**2
1995        hybri_limit(1,2)=tbi12**2+tbi12b**2
1996        hybri_limit(2,1)=tbi21**2+tbi21b**2
1997      end if
1998      write(132,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1))
1999      write(133,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2))
2000      write(134,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1))
2001      write(135,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2))
2002      write(1234,*) omega, real(fw1(ifreq,1)),aimag(fw1(ifreq,1))
2003    end do
2004    ! Built level and limit of hybridization
2005    !=======================================
2006    levels_ctqmc(1:nflavor)=-umod(1,1)/two
2007 
2008    write(std_out,*) "fw1_nd"
2009    write(std_out,*) fw1_nd(1,1,1), fw1_nd(1,1,2)
2010    write(std_out,*) fw1_nd(1,2,1), fw1_nd(1,2,2)
2011    write(std_out,*) "fw1"
2012    write(std_out,*) fw1(1,1), fw1(1,2)
2013    write(std_out,*) fw1(2,1), fw1(2,2)
2014 
2015  ! Rotate hybridization if testrot=1
2016  !==================================
2017    if(testrot==1) then
2018 
2019      do ifreq=1,dmftqmc_l
2020        RRt(:,:)  = MATMUL ( RR(:,:)  , fw1_nd(ifreq,:,:) )
2021    !write(6,*) "RRt"
2022    !write(6,*) RRt(1,1), RRt(1,2)
2023    !write(6,*) RRt(2,1), RRt(2,2)
2024        RR1(:,:)  = MATMUL ( RRt(:,:) , RRi(:,:)          )
2025    !write(6,*) "RR1"
2026    !write(6,*) RR1(1,1), RR1(1,2)
2027    !write(6,*) RR1(2,1), RR1(2,2)
2028        fw1_nd(ifreq,:,:)=RR1(:,:)
2029        omega=pi*temp*(two*float(ifreq)+1)
2030        write(3322,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1))
2031        write(232,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1))
2032        write(233,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2))
2033        write(234,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1))
2034        write(235,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2))
2035      end do
2036 
2037      ! Rotate limit of hybridization
2038      !=======================================
2039      RRt(:,:)  = MATMUL ( RR(:,:)  , hybri_limit(:,:)  )
2040      RR1(:,:)  = MATMUL ( RRt(:,:) , RRi(:,:)          )
2041      hybri_limit(:,:)=RR1(:,:)
2042 
2043    end if
2044    ! rajouter test real(fw1_nd(1,:,:)) doit etre diagonale
2045 
2046  !======================================
2047  ! Rotate Green's function from CTQMC
2048  !======================================
2049  else if(opt==2) then
2050 
2051    write(std_out,*) "gw_tmp_nd"
2052    write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2)
2053    write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2)
2054    ! Rotate Green's function back
2055    !==============================
2056    if(testrot==1) then
2057      do ifreq=1,dmftqmc_l
2058        RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gw_tmp_nd(ifreq,1:nflavor,1:nflavor) )
2059        RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) )
2060        gw_tmp_nd(ifreq,1:nflavor,1:nflavor)=RR1(1:nflavor,1:nflavor)
2061      end do
2062 
2063      write(std_out,*) "gw_tmp_nd after rotation"
2064      write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2)
2065      write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2)
2066 
2067      do itau=1,dmftqmc_l
2068        RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) )
2069        RR1(1:nflavor,1:nflavor)  = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) )
2070        gtmp_nd(itau,1:nflavor,1:nflavor)=real(RR1(1:nflavor,1:nflavor))
2071      end do
2072 
2073    ! Rotate Green's function for comparison with testrot=1
2074    !======================================================
2075    else if (testrot==0) then ! produce rotated green's function to compare to testrot=1 case
2076 
2077      do itau=1,dmftqmc_l
2078        RRt(1:nflavor,1:nflavor) = MATMUL ( RR(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) )
2079        RR1(1:nflavor,1:nflavor)  = MATMUL ( RRt(1:nflavor,1:nflavor),RRi(1:nflavor,1:nflavor) )
2080        write(444,*) real(itau-1)/(temp*real(dmftqmc_l)),real(RR1(1,1)),real(RR1(2,2)),real(RR1(1,2)),real(RR1(2,1))
2081      end do
2082 
2083    end if
2084 
2085    ! Print out rotated Green's function
2086    !=====================================
2087    do itau=1,dmftqmc_l
2088      write(555,'(e14.5,4(2e14.5,3x))') real(itau-1)/(temp*real(dmftqmc_l)),gtmp_nd(itau,1,1),&
2089 &     gtmp_nd(itau,2,2),gtmp_nd(itau,1,2),gtmp_nd(itau,2,1)
2090    end do
2091 
2092    write(message,'(2a)') ch10,' testcode end of test calculation'
2093    ABI_ERROR(message)
2094 
2095  end if
2096  close(444)
2097  close(555)
2098 
2099 end subroutine testcode_ctqmc

m_forctqmc/testcode_ctqmc_b [ Functions ]

[ Top ] [ m_forctqmc ] [ Functions ]

NAME

 testcode_ctqmc_b

FUNCTION

 Setup ultra simple hybridization to test CTQMC in simple situations.

INPUTS

 temp = temperature
 dmftqmc_l = number of times slices
 levels_ctqmc_nd=level matrix

OUTPUT

 fw1_nd=hybridization matrix
 umod = value of U
 hybri_limit= limit of F
 weiss_for_rot= weiss function
 hybri_coeff

SIDE EFFECTS

NOTES

SOURCE

1753 subroutine testcode_ctqmc_b(energy_level,hybri_coeff,weiss_for_rot,dmftqmc_l,fw1_nd,levels_ctqmc,&
1754 &   levels_ctqmc_nd,hybri_limit,temp,umod,opt_diag,opt_fk)
1755 
1756 !Arguments ------------------------------------
1757 !scalars
1758  integer, intent(in) :: dmftqmc_l,opt_diag,opt_fk
1759  real(dp), intent(in) :: temp
1760  real(dp), intent(out) :: umod(2,2)
1761  real(dp), intent(inout) :: levels_ctqmc(:)
1762  complex(dpc), intent(out) :: fw1_nd(:,:,:)
1763  complex(dpc),  intent(inout) :: levels_ctqmc_nd(:,:)
1764  complex(dpc),  intent(inout) :: hybri_limit(:,:)
1765  type(oper_type)  :: energy_level
1766  type(matlu_type), allocatable  :: hybri_coeff(:)
1767  type(green_type)  :: weiss_for_rot
1768 
1769 !Local variables ------------------------------
1770  integer :: ifreq,iatom,ima,imb,ispa,ispb
1771  real(dp) :: omega
1772  real(dp) :: facnd, facd
1773  character(len=30) :: tmpfil
1774 ! ************************************************************************
1775  facnd=0.8d0
1776  facd=1.0d0
1777  !write(6,*) "fac",facnd,facd
1778  levels_ctqmc_nd(2,2)   = energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb)
1779  levels_ctqmc_nd(1,1)   = energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa)
1780  levels_ctqmc(2)   = real(energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb),kind=dp)
1781  levels_ctqmc(1)   = real(energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa),kind=dp)
1782  if(opt_diag/=1) then
1783    levels_ctqmc_nd(1,2)   = energy_level%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1784    levels_ctqmc_nd(2,1)   = energy_level%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1785  end if
1786  hybri_limit(1,1)  = facd*hybri_coeff(iatom)%mat(ima,ima,1,ispa,ispa)
1787  hybri_limit(2,2)  = facd*hybri_coeff(iatom)%mat(imb,imb,1,ispb,ispb)
1788  hybri_limit(1,2)  = facnd*hybri_coeff(iatom)%mat(ima,imb,1,ispa,ispb)
1789  hybri_limit(2,1)  = facnd*hybri_coeff(iatom)%mat(imb,ima,1,ispb,ispa)
1790  !write(6,*) "hybri_limit",hybri_limit
1791  !write(6,*) "levels_ctqmc",levels_ctqmc
1792  umod=zero
1793 
1794  tmpfil = 'fw1_nd_re'
1795  !if (open_file(newunit=unt,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then
1796  !  ABI_ERROR(message)
1797  !end if
1798  tmpfil = 'fw1_nd_im'
1799  !if (open_file(newunit=unt2,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then
1800  !  ABI_ERROR(message)
1801  !end if
1802  write(std_out,*) "testcode==2",ispa,ispb,ima,imb
1803  write(std_out,*) "opt_fk==",opt_fk
1804  do ifreq=1,dmftqmc_l
1805    if (opt_fk==1) then
1806      fw1_nd(ifreq,1,1) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa)
1807      fw1_nd(ifreq,2,2) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb)
1808      !fw1_nd(ifreq,1,2) =  weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1809      !fw1_nd(ifreq,2,1) =  weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1810      fw1_nd(ifreq,1,2) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1811      fw1_nd(ifreq,2,1) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1812      omega=pi*temp*(two*float(ifreq)-1)
1813    else if (opt_fk==0) then
1814      fw1_nd(ifreq,1,1) =  facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa)
1815      fw1_nd(ifreq,2,2) =  facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb)
1816      fw1_nd(ifreq,1,2) =  facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1817      fw1_nd(ifreq,2,1) =  facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1818      call xginv(fw1_nd(ifreq,:,:),2)
1819    end if
1820  end do
1821 end subroutine testcode_ctqmc_b