TABLE OF CONTENTS


ABINIT/m_forctqmc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_forctqmc

FUNCTION

 Prepare CTQMC and call CTQMC

COPYRIGHT

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

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

25 #if defined HAVE_CONFIG_H
26 #include "config.h"
27 #endif
28 
29 #include "abi_common.h"
30 
31 MODULE m_forctqmc
32 
33  use defs_basis
34 
35  implicit none
36 
37  private
38 
39  public :: qmc_prep_ctqmc
40  public :: testcode_ctqmc

m_forctqmc/qmc_prep_ctqmc [ Functions ]

[ Top ] [ m_forctqmc ] [ Functions ]

NAME

 qmc_prep_ctqmc

FUNCTION

 Prepare and call the qmc subroutines

COPYRIGHT

 Copyright (C) 1999-2018 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

PARENTS

      impurity_solve

CHILDREN

      add_matlu,checkreal_matlu,compute_levels,copy_green,copy_matlu
      ctqmc_triqs_run,ctqmcinterface_finalize,ctqmcinterface_init
      ctqmcinterface_run,ctqmcinterface_setopts,data4entropydmft_setdocc
      destroy_green,destroy_matlu,destroy_oper,diag_matlu,diff_matlu
      fac_matlu,flush_unit,fourier_green,hybridization_asymptotic_coefficient
      identity_matlu,init_green,init_matlu,init_oper,int_fct,inverse_oper
      jbessel,occup_green_tau,print_green,print_matlu,printocc_green
      printplot_matlu,prod_matlu,rotate_matlu,rotatevee_hu,sbf8,shift_matlu
      slm2ylm_matlu,sym_matlu,testcode_ctqmc,vee_ndim2tndim_hu_r,wrtout,xginv
      xmpi_barrier,xmpi_bcast

SOURCE

  88 subroutine qmc_prep_ctqmc(cryst_struc,green,self,hu,paw_dmft,pawang,pawprtvol,weiss)
  89 
  90  use defs_basis
  91  use defs_datatypes
  92  use defs_abitypes
  93  use m_errors
  94  use m_xmpi
  95 
  96  use m_special_funcs, only : sbf8
  97  use m_pawang, only : pawang_type
  98  use m_crystal, only : crystal_t
  99  use m_green, only : green_type,occup_green_tau,print_green,printocc_green,spline_fct,copy_green,init_green,destroy_green,&
 100 & int_fct,greenldacompute_green,fourier_green
 101  use m_paw_dmft, only : paw_dmft_type
 102  use m_hide_lapack,         only : xginv
 103  use m_oper, only : oper_type,destroy_oper,init_oper,inverse_oper
 104  use m_self, only : self_type
 105  use m_matlu, only : matlu_type,sym_matlu, print_matlu, &
 106 & diag_matlu,init_matlu,destroy_matlu,rotate_matlu,checkdiag_matlu,checkreal_matlu, &
 107 & copy_matlu, diff_matlu, slm2ylm_matlu, shift_matlu, prod_matlu,fac_matlu,&
 108 & add_matlu,printplot_matlu,identity_matlu
 109  use m_hu, only : hu_type,rotatevee_hu,vee_ndim2tndim_hu_r
 110  use m_Ctqmc
 111  use m_CtqmcInterface
 112  use m_Ctqmcoffdiag
 113  use m_CtqmcoffdiagInterface
 114  use m_GreenHyb
 115  use m_data4entropyDMFT
 116  !use m_self, only : self_type,initialize_self,destroy_self,print_self,rw_self
 117  use m_io_tools, only : flush_unit, open_file
 118  use m_paw_numeric, only : jbessel=>paw_jbessel
 119  use m_datafordmft, only : hybridization_asymptotic_coefficient,compute_levels
 120 
 121 #if defined HAVE_TRIQS_v2_0 || defined HAVE_TRIQS_v1_4
 122  use TRIQS_CTQMC !Triqs module
 123 #endif
 124  use ISO_C_BINDING
 125 
 126 !This section has been created automatically by the script Abilint (TD).
 127 !Do not modify the following lines by hand.
 128 #undef ABI_FUNC
 129 #define ABI_FUNC 'qmc_prep_ctqmc'
 130 !End of the abilint section
 131 
 132  implicit none
 133 
 134 !Arguments ------------------------------------
 135 !scalars
 136 ! type(pawang_type), intent(in) :: pawang
 137  type(crystal_t),intent(in) :: cryst_struc
 138  type(green_type), intent(inout) :: green  ! MGNAG: This fix the problem with v7[27:29] on nag@petrus
 139  type(hu_type), intent(in) :: hu(cryst_struc%ntypat)
 140  type(paw_dmft_type), intent(inout)  :: paw_dmft
 141  type(pawang_type), intent(in) :: pawang
 142  integer, intent(in) :: pawprtvol
 143  type(green_type), intent(inout) :: weiss
 144  type(self_type), intent(in) :: self
 145 
 146 !Local variables ------------------------------
 147  character(len=500) :: message
 148  character(len=2) :: gtau_iter,iatomnb
 149  integer :: iatom,ierr,if1,if2,iflavor,iflavor1,iflavor2,iflavor3,ifreq,im,im1,ispinor,ispinor1,isppol,itau,itypat,im2,ispinor2
 150  integer :: lpawu,master,mbandc,natom,nflavor,nkpt,nspinor,nsppol,nsppol_imp,tndim,ispa,ispb,ima,imb
 151  integer :: nproc,opt_diag,opt_nondiag,testcode,testrot,dmft_nwlo,opt_fk,useylm,nomega,opt_rot
 152  integer :: ier,rot_type_vee
 153  complex(dpc) :: omega_current,integral(2,2)
 154  real(dp) :: doccsum,noise,omega
 155  real(dp) :: facd,facnd
 156  logical :: nondiaglevels
 157 ! arrays
 158  real(dp), allocatable :: docc(:,:)
 159  real(dp), allocatable, target :: gtmp(:,:), levels_ctqmc(:) !modif
 160  complex(dpc), allocatable :: levels_ctqmc_nd(:,:)
 161  complex(dpc), allocatable :: hybri_limit(:,:)
 162  real(dp), allocatable, target :: gtmp_nd(:,:,:)
 163  real(dp) :: umod(2,2)
 164  character(len=30) :: tmpfil
 165  complex(dpc), allocatable :: fw1(:,:),gw_tmp(:,:)
 166  complex(dpc), allocatable, target :: gw_tmp_nd(:,:,:) !modif
 167  complex(dpc), allocatable, target :: fw1_nd(:,:,:) !modif
 168  complex(dpc), allocatable :: gw1_nd(:,:,:)
 169  complex(dpc), allocatable :: shift(:)
 170  integer,parameter :: optdb=0
 171  type(coeff2_type), allocatable :: udens_atoms(:)
 172 ! Type    -----------------------------------------
 173  type(coeff2c_type), allocatable :: eigvectmatlu(:,:)
 174  type(green_type)  :: weiss_for_rot
 175  type(matlu_type), allocatable :: dmat_diag(:)
 176  type(matlu_type), allocatable :: matlu1(:)
 177  type(matlu_type), allocatable :: matlu2(:)
 178  type(matlu_type), allocatable :: matlu3(:)
 179  type(matlu_type), allocatable :: matlu4(:)
 180  type(matlu_type), allocatable :: identity(:)
 181  type(matlu_type), allocatable :: level_diag(:)
 182  type(oper_type)  :: energy_level
 183  !type(self_type) :: self
 184 ! type(green_type) :: gw_loc
 185  type(CtqmcInterface) :: hybrid   !!! WARNING THIS IS A BACKUP PLAN
 186  type(CtqmcoffdiagInterface) :: hybridoffdiag   !!! WARNING THIS IS A BACKUP PLAN
 187  type(green_type) :: greenlda
 188  type(matlu_type), allocatable  :: hybri_coeff(:)
 189 ! Var added to the code for TRIQS_CTQMC test and default value -----------------------------------------------------------
 190  logical(kind=1) :: rot_inv = .false.
 191  logical(kind=1) :: leg_measure = .true.
 192 #if defined HAVE_TRIQS_v2_0 || defined HAVE_TRIQS_v1_4
 193  logical(kind=1) :: hist = .false.
 194  logical(kind=1) :: wrt_files = .true.
 195  logical(kind=1) :: tot_not = .true.
 196 #endif
 197 
 198  integer :: nfreq,unt,unt2
 199  integer :: ntau ! >= 2*nfreq + 1
 200  integer :: nleg
 201  integer :: ileg
 202  integer :: verbosity_solver ! min 0 -> max 3
 203 
 204  real(dp) :: beta,besp,bespp,xx
 205  complex(dpc) :: u_nl
 206 
 207  complex(dpc), allocatable, target ::fw1_nd_tmp(:,:,:)
 208  complex(dpc), allocatable, target :: g_iw(:,:,:)
 209  real(dp), allocatable, target :: u_mat_ij(:,:)
 210  real(dp), allocatable, target :: u_mat_ijkl(:,:,:,:)
 211  real(dp), allocatable, target :: u_mat_ijkl_tmp(:,:,:,:)
 212 !  real(dp), allocatable, target :: gtau(:,:,:)
 213  real(dp), allocatable, target :: gl_nd(:,:,:)
 214  real(dp), allocatable :: jbes(:)
 215 
 216  type(c_ptr) :: levels_ptr, fw1_nd_ptr, u_mat_ij_ptr, u_mat_ijkl_ptr, g_iw_ptr, gtau_ptr, gl_ptr
 217 
 218 ! ************************************************************************
 219  mbandc=paw_dmft%mbandc
 220  nkpt=paw_dmft%nkpt
 221  nsppol=paw_dmft%nsppol
 222  natom=paw_dmft%natom
 223  nspinor=paw_dmft%nspinor
 224  greenlda%whichgreen="LDA"
 225 
 226  call init_green(weiss_for_rot,paw_dmft,opt_oper_ksloc=2)
 227 ! weiss_for_rot=>weiss
 228 ! call init_green(gw_loc,paw_dmft)
 229  call copy_green(weiss,weiss_for_rot,opt_tw=2)
 230 !=======================================================================
 231 !== Use one QMC solver   ===============================================
 232 !=======================================================================
 233  write(message,'(2a)') ch10,'  ===  CT-QMC solver === '
 234  call wrtout(std_out,message,'COLL')
 235 
 236 ! Initialise for compiler
 237  omega_current=czero
 238 
 239 ! Initialise nproc
 240  nproc=paw_dmft%nproc
 241 
 242 ! ======================================
 243 ! Allocations: diagonalization and eigenvectors
 244 ! ======================================
 245  ABI_DATATYPE_ALLOCATE(udens_atoms,(natom))
 246  ABI_DATATYPE_ALLOCATE(eigvectmatlu,(natom,nsppol))
 247  ABI_DATATYPE_ALLOCATE(dmat_diag,(natom))
 248  ABI_DATATYPE_ALLOCATE(identity,(natom))
 249  call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,dmat_diag)
 250  call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,identity)
 251  call identity_matlu(identity,natom)
 252  do iatom=1,cryst_struc%natom
 253    lpawu=paw_dmft%lpawu(iatom)
 254    if(lpawu/=-1) then
 255      tndim=nspinor*(2*lpawu+1)
 256      do isppol=1,nsppol
 257        ABI_ALLOCATE(eigvectmatlu(iatom,isppol)%value,(tndim,tndim))
 258      end do
 259      ABI_ALLOCATE(udens_atoms(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1)))
 260      dmat_diag(iatom)%mat=czero
 261    end if
 262  end do
 263 
 264 ! ___________________________________________________________________________________
 265 !
 266 !  FIRST PART: DIAGONALISATION AND ROTATIONS.
 267 ! ___________________________________________________________________________________
 268 !
 269 
 270 ! =================================================================
 271 ! Choose to diagonalize and how to do it
 272 ! =================================================================
 273 
 274 ! =================================================================
 275 ! Impose diago of density matrix
 276 ! =================================================================
 277 
 278 ! =================================================================
 279 ! Impose diago of levels and Ylm basis if opt_nondiag=1
 280 ! =================================================================
 281 ! opt_diag=1 ! 1: diago the levels (The best choice).
 282 ! opt_diag=2 ! 2: diago density matrix (can be used for historical reasons)
 283 
 284 !  Need in the general case of two input variable for opt_diag and
 285 !  opt_nondiag!
 286 !  The default value of opt_diag should be 2 for historical reasons (or
 287 !  we decide to change the automatic tests)
 288 !  opt_nondiag should be 0 by default
 289  opt_diag    = 1
 290  if(paw_dmft%dmft_solv>=6)  then
 291    opt_nondiag = 1 ! Use cthyb in triqs or ctqmc in abinit with offdiag terms in F
 292  else
 293    opt_nondiag = 0 ! use fast ctqmc in ABINIT without off diagonal terms in F
 294  end if
 295 
 296  useylm=0
 297  if(nspinor==2) then
 298    useylm=1      ! to avoid complex G(tau)
 299  end if
 300 
 301  !write(6,*) "nspinor,useylm",nspinor,useylm
 302  if(useylm==0) then
 303    write(std_out,*) " Slm basis is used (before rotation)"
 304    rot_type_vee=1 ! for rotatevee_hu
 305  else if(useylm==1) then
 306    write(std_out,*) " Ylm basis is used (before rotation)"
 307    rot_type_vee=4 ! for rotatevee_hu
 308  end if
 309 
 310 
 311 ! if(useylm==1.and.opt_diag/=1) MSG_ERROR("useylm==1 and opt_diag/=0 is not possible")
 312  if(hu(1)%jpawu_zero.and.nsppol==2) nsppol_imp=2 ! J=0 and nsppol=2
 313  if(.not.hu(1)%jpawu_zero.or.nsppol/=2) nsppol_imp=1  ! J/=0 ou nsppol=1
 314 ! =================================================================
 315 ! Compute LDA Green's function to compare to weiss_for_rot (check)
 316 ! =================================================================
 317 ! call init_green(greenlda,paw_dmft,opt_oper_ksloc=3)
 318 ! call greenldacompute_green(cryst_struc,greenlda,pawang,paw_dmft)
 319 !! call copy_green(greenlda,weiss_for_rot,2)
 320 
 321 ! =================================================================
 322 ! Compute atomic levels
 323 ! =================================================================
 324  call init_oper(paw_dmft,energy_level,opt_ksloc=3)
 325 
 326  ! Compute atomic levels in Slm basis
 327  ! ----------------------------------
 328  call compute_levels(cryst_struc,energy_level,self%hdc,pawang,paw_dmft,nondiag=nondiaglevels)
 329 
 330  ! If levels are not diagonal, then diagonalize it (according to
 331  ! dmftctqmc_basis)
 332  ! ------------------------------------------------
 333  if(paw_dmft%dmftctqmc_basis==1) then
 334    if(nondiaglevels.or.useylm==1) then
 335      opt_diag=1
 336      write(message,'(3a)') ch10, "   == Hamiltonian in local basis is non diagonal: diagonalise it",ch10
 337    else
 338      opt_diag=0
 339      write(message,'(5a)') ch10, "   == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 &
 340 &     ,"      CTQMC will use this basis",ch10
 341    end if
 342  else if (paw_dmft%dmftctqmc_basis==2) then
 343    if(nondiaglevels.or.useylm==1) then
 344      write(message,'(7a)') ch10, "   == Hamiltonian in local basis is non diagonal",ch10, &
 345 &     "   == According to dmftctqmc_basis: diagonalise density matrix",ch10, &
 346 &     "   == Warning : Check that the Hamiltonian is diagonal !",ch10
 347      opt_diag=2
 348    else
 349      write(message,'(5a)') ch10, "   == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 &
 350 &     ,"      CTQMC will use this basis",ch10
 351      opt_diag=0
 352    end if
 353  else if (paw_dmft%dmftctqmc_basis==0) then
 354    if(nondiaglevels) then
 355      write(message,'(4a)') ch10, "   == Hamiltonian in local basis is non diagonal",ch10, &
 356 &     "   == According to dmftctqmc_basis: keep this non diagonal basis for the calculation"
 357    else
 358      write(message,'(5a)') ch10, "   == Hamiltonian in local basis is diagonal in the Slm basis ",ch10 &
 359 &     ,"      CTQMC will use this basis",ch10
 360    end if
 361    opt_diag=0
 362  end if
 363  call wrtout(std_out,message,'COLL')
 364  if(opt_diag==1) then
 365    write(std_out,*) "  ==  The atomic levels are diagonalized"
 366  else if(opt_diag==2) then
 367    write(std_out,*) "  ==  The correlated occupation matrix is diagonalized"
 368  end if
 369 
 370 ! =================================================================
 371 ! Now, check if diagonalisation is necessary
 372 ! =================================================================
 373 
 374 
 375 ! =================================================================
 376 ! First rotate to Ylm basis the atomic levels
 377 ! =================================================================
 378 
 379  if(useylm==1) then
 380 
 381    ! Rotate from Slm to Ylm the atomic levels
 382    ! ----------------------------------------
 383    call slm2ylm_matlu(energy_level%matlu,natom,1,pawprtvol)
 384 
 385    ! Print atomic energy levels in Ylm basis
 386    ! --------------------------------
 387    if(pawprtvol>=3) then
 388      write(message,'(a,a)') ch10, " == Print Energy levels in Ylm basis"
 389      call wrtout(std_out,message,'COLL')
 390      call print_matlu(energy_level%matlu,natom,1)
 391    end if
 392 
 393  end if ! useylm
 394 
 395 ! ===========================================================================================
 396 ! Start for diagonalization of levels/density matrix according to opt_diag
 397 ! ===========================================================================================
 398  !opt_rot=2 ! do it one time before CTQMC
 399  opt_rot=1 ! do all the rotations successively on all different quantities.
 400  if(opt_diag==1.or.opt_diag==0) then
 401 
 402 
 403    if(opt_diag==1) then
 404 ! =================================================================
 405 ! Diagonalize atomic levels
 406 ! =================================================================
 407      ABI_DATATYPE_ALLOCATE(level_diag,(natom))
 408      call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,level_diag)
 409 
 410      ! Diagonalise atomic levels (opt_real is necessary, because
 411      ! rotation must be real in order for the occupations and Green's
 412      ! function to be real)
 413      ! ---------------------------------------------------------------
 414      call diag_matlu(energy_level%matlu,level_diag,natom,&
 415 &     prtopt=pawprtvol,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,optreal=1,&
 416 &     test=paw_dmft%dmft_solv)  ! temporary: test should be extended to all cases.
 417 
 418 !     call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1)
 419 !       write(message,'(a,2x,a,f13.5)') ch10,&
 420 !&       " == Print first Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie
 421 !       call wrtout(std_out,message,'COLL')
 422 !       call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1)
 423 
 424      if(opt_rot==1) call copy_matlu(level_diag,energy_level%matlu,natom)
 425 
 426 
 427      call destroy_matlu(level_diag,natom)
 428      ABI_DATATYPE_DEALLOCATE(level_diag)
 429 
 430      ! Print diagonalized levels
 431      ! --------------------------
 432      if(pawprtvol>=3) then
 433        write(message,'(a,2x,a,f13.5)') ch10,&
 434 &       " == Print Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie
 435        call wrtout(std_out,message,'COLL')
 436        call print_matlu(energy_level%matlu,natom,1,compl=1,opt_exp=1)
 437      else
 438        write(message,'(a,2x,a,f13.5)') ch10,&
 439 &       " == Energy levels Diagonalized for Fermi Level=",paw_dmft%fermie
 440        call wrtout(std_out,message,'COLL')
 441      end if
 442      call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee)
 443    else if (opt_diag==0) then
 444      do iatom=1,cryst_struc%natom
 445        lpawu=paw_dmft%lpawu(iatom)
 446        itypat=cryst_struc%typat(iatom)
 447        if(lpawu/=-1) then
 448        !  write(6,*) size(udens_atoms(iatom)%value)
 449        !  write(6,*) size(hu(itypat)%udens)
 450        !  write(6,*) udens_atoms(iatom)%value
 451        !  write(6,*) hu(itypat)%udens
 452          udens_atoms(iatom)%value=hu(itypat)%udens
 453        end if
 454      end do
 455    end if
 456   ! call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms)
 457 
 458  else if(opt_diag==2) then
 459 ! =================================================================
 460 ! Diagonalizes density matrix and keep eigenvectors in eigvectmatlu
 461 ! =================================================================
 462 
 463    ! Print density matrix before diagonalization
 464    ! -------------------------------------------
 465    if(pawprtvol>=3) then
 466      write(message,'(a,2x,a)') ch10,        " == Density Matrix before diagonalisation ="
 467      call wrtout(std_out,message,'COLL')
 468      !MGNAG: This call is wrong if green has intent(out), now we use intent(inout)
 469      call print_matlu(green%occup%matlu,natom,1)
 470    end if
 471 
 472 !!  checkstop: we can have two different diagonalisation basis for the up and dn
 473 !!  but one use the same basis, unless the error is really to large(>0.1)
 474 
 475    ! Diagonalize density matrix
 476    ! ---------------------------
 477    call diag_matlu(green%occup%matlu,dmat_diag,natom,&
 478 &   prtopt=4,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,checkstop=.false.)
 479 
 480    ! Print diagonalized density matrix
 481    ! ----------------------------------
 482    if(pawprtvol>=3) then
 483      write(message,'(a,2x,a)') ch10,&
 484 &     " == Diagonalized Density Matrix in the basis used for QMC ="
 485      call wrtout(std_out,message,'COLL')
 486      call print_matlu(dmat_diag,natom,1)
 487 
 488      !write(message,'(2a,i3,13x,a)') ch10,'    ==  Rotation of interaction matrix =='
 489      !call wrtout(std_out,message,'COLL')
 490    end if
 491 
 492    !if (.not.hu(1)%jpawu_zero) &
 493    !MSG_WARNING("In qmc_prep_ctqmc J/=0 and rotation matrix not rotated")
 494 !  Rotate interaction.
 495 !   call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms)
 496    call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,rot_type_vee)
 497 
 498  end if
 499 ! ===========================================================================================
 500 ! END Of diagonalization
 501 ! ===========================================================================================
 502 
 503  call flush_unit(std_out)
 504 
 505 ! ===========================================================================================
 506 ! Broadcast matrix of rotation from processor 0 to the other
 507 ! In case of degenerate levels, severals rotations are possible. Here we
 508 ! choose the rotation of proc 0. It is arbitrary.
 509 ! ===========================================================================================
 510  do iatom=1,cryst_struc%natom
 511    lpawu=paw_dmft%lpawu(iatom)
 512    if(lpawu/=-1) then
 513      tndim=nspinor*(2*lpawu+1)
 514      do isppol=1,nsppol
 515        call xmpi_bcast(eigvectmatlu(iatom,isppol)%value,0,paw_dmft%spacecomm,ier)
 516      end do
 517    end if
 518  end do
 519 
 520 
 521      !unitnb=300000+paw_dmft%myproc
 522      !call int2char4(paw_dmft%myproc,tag_proc)
 523      !tmpfil = 'eigvectmatluaftermpi'//tag_proc
 524      !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
 525      !do iflavor1=1,14
 526      !  do iflavor2=1,14
 527      !    write(unitnb,*) iflavor1,iflavor2,eigvectmatlu(1,1)%value(iflavor1,iflavor2)
 528      !  enddo
 529      !enddo
 530 
 531 ! ===========================================================================================
 532 ! Now rotate various quantities in the new basis
 533 ! ===========================================================================================
 534 
 535 !=======================================================
 536 ! Allocate, Compute, and Rotate atomic levels for CTQMC
 537 !=======================================================
 538 
 539    ! If levels not rotated, rotate them
 540    ! -----------------------------------
 541  if(opt_diag==2.and.opt_rot==1) call rotate_matlu(energy_level%matlu,eigvectmatlu,natom,3,1)
 542 
 543    ! Print atomic levels
 544    ! -------------------
 545  if(pawprtvol>=3) then
 546    write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels after rotation"
 547    call wrtout(std_out,message,'COLL')
 548    call print_matlu(energy_level%matlu,natom,1)
 549  else
 550    write(message,'(a,2x,a,f13.5)') ch10," == CT-QMC Energy levels rotated"
 551    call wrtout(std_out,message,'COLL')
 552  end if
 553 
 554 !====================================================================
 555 ! If levels were diagonalized before, then rotate density matrix for
 556 ! information.
 557 !====================================================================
 558  if(opt_diag==1) then
 559    ABI_DATATYPE_ALLOCATE(matlu1,(natom))
 560    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
 561    call copy_matlu(green%occup%matlu,matlu1,natom)
 562    if(pawprtvol>=3) then
 563      write(message,'(a,2x,a)') ch10,&
 564 &     " == Occupations before rotations"
 565      call wrtout(std_out,message,'COLL')
 566      call print_matlu(green%occup%matlu,natom,1)
 567    end if
 568 
 569    ! 1) rotate density matrix to Ylm basis
 570    ! --------------------------------------
 571    if(useylm==1) then
 572      call slm2ylm_matlu(matlu1,natom,1,pawprtvol)
 573      if(pawprtvol>=3) then
 574        write(message,'(a,a)') ch10, " == Print occupations in Ylm basis"
 575        call wrtout(std_out,message,'COLL')
 576        call print_matlu(matlu1,natom,1)
 577      end if
 578    end if
 579 
 580    ! 2) rotate density matrix to rotated basis
 581    ! -------------------------------------------
 582    if(opt_rot==1.or.opt_rot==2) call rotate_matlu(matlu1,eigvectmatlu,natom,3,1)
 583    write(message,'(a,2x,a,f13.5)') ch10," == Rotated occupations (for information)"
 584    call wrtout(std_out,message,'COLL')
 585    call print_matlu(matlu1,natom,1,compl=1)
 586    call checkreal_matlu(matlu1,natom,tol10)
 587    call destroy_matlu(matlu1,natom)
 588    ABI_DATATYPE_DEALLOCATE(matlu1)
 589 
 590  end if
 591 
 592  call flush_unit(std_out)
 593 
 594 
 595 ! =================================================================
 596 ! Rotate weiss function according to eigenvectors.
 597 ! =================================================================
 598 !!!stop
 599   ! Rotate Weiss function first in Ylm basis
 600   ! -------------------------------------------------------------------
 601  if(useylm==1) then
 602    write(message,'(a,2x,a)') ch10, " == Rotation of weiss and greenlda in the Ylm Basis="
 603    call wrtout(std_out,message,'COLL')
 604    do ifreq=1,paw_dmft%dmft_nwlo
 605      call slm2ylm_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,1,0)
 606      call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,1,0)
 607      ! call slm2ylm_matlu(greenlda%oper(ifreq)%matlu,natom,1,0)
 608    end do
 609  end if
 610 
 611  if(pawprtvol>=3) then
 612    !   write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 613    !   " == Print weiss for small freq 1 before rot" ! debug
 614    !   call wrtout(std_out,message,'COLL') ! debug
 615    !   call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1) !  debug
 616 
 617     ! Print Weiss function
 618     ! --------------------
 619    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 620 &  " == Print weiss for 1st freq before rot" ! debug
 621    call wrtout(std_out,message,'COLL') ! debug
 622    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) !  debug
 623    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 624 &  " == Print weiss for last freq before rot" ! debug
 625    call wrtout(std_out,message,'COLL') ! debug
 626    call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) !  debug
 627 !    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 628 !&   " == Print LDA G for 1st freq before rot" ! debug
 629 !    call wrtout(std_out,message,'COLL') ! debug
 630 !    call print_matlu(greenlda%oper(1)%matlu,natom,1,compl=1,opt_exp=2) !  debug
 631 !    write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 632 !&   " == Print LDA G for last freq before rot" ! debug
 633 !    call wrtout(std_out,message,'COLL') ! debug
 634 !    call print_matlu(greenlda%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) !  debug
 635  end if
 636 
 637  if(opt_diag/=0) then
 638    ! Rotate Weiss function first in Ylm basis then in the rotated basis.
 639    ! -------------------------------------------------------------------
 640    write(message,'(a,2x,a)') ch10, " == Rotation of weiss ="
 641    call wrtout(std_out,message,'COLL')
 642    do ifreq=1,paw_dmft%dmft_nwlo
 643      if(opt_rot==1) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 644      if(opt_rot==1) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 645 !    call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6)
 646    end do
 647 
 648    if(paw_dmft%myproc .eq. mod(nproc+1,nproc)) then
 649      if (open_file(trim(paw_dmft%filapp)//"_atom__G0w_.dat", message, newunit=unt) /= 0) then
 650        MSG_ERROR(message)
 651      end if
 652      do ifreq=1,paw_dmft%dmft_nwlo
 653        write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),&
 654 &       (((weiss_for_rot%oper(ifreq)%matlu(1)%mat(im1,im1,isppol,ispinor,ispinor),&
 655 &       im1=1,3),ispinor=1,nspinor),isppol=1,nsppol)
 656      end do
 657      close(unt)
 658    end if
 659 
 660    call flush_unit(std_out)
 661    if(pawprtvol>=3) then
 662      write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 663 &    " == Print weiss for small freq 1 after rot" ! debug
 664      call wrtout(std_out,message,'COLL') ! debug
 665      call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1) !  debug
 666      write(message,'(a,2x,a,f13.5)') ch10,&   ! debug
 667 &    " == Print weiss for last freq after rot"   ! debug
 668      call wrtout(std_out,message,'COLL')   ! debug
 669      call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug
 670    end if
 671 
 672 !   ! Rotate LDA Green's function first in Ylm basis then in the rotated basis and compare to weiss_for_rot
 673 !   ! -----------------------------------------------------------------------------------------------------
 674 !   write(message,'(a,2x,a)') ch10, " == Rotation of greenlda ="
 675 !   call wrtout(std_out,message,'COLL')
 676 !   do ifreq=1,paw_dmft%dmft_nwlo
 677 !     if(opt_rot==1) call rotate_matlu(greenlda%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 678 !     call diff_matlu("Weiss_for_rot","greenlda",weiss_for_rot%oper(ifreq)%matlu,greenlda%oper(ifreq)%matlu,natom,1,tol14)
 679 !!    call checkdiag_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,tol6)
 680 !   end do
 681 !   if(pawprtvol>=3) then
 682 !     write(message,'(a,2x,a,f13.5)') ch10,& ! debug
 683 !&    " == Print greenlda for small freq 1 after rot" ! debug
 684 !     call wrtout(std_out,message,'COLL') ! debug
 685 !     call print_matlu(greenlda%oper(1)%matlu,natom,1,compl=1,opt_exp=2) !  debug
 686 !     write(message,'(a,2x,a,f13.5)') ch10,&   ! debug
 687 !&    " == Print greenlda for last freq after rot"   ! debug
 688 !     call wrtout(std_out,message,'COLL')   ! debug
 689 !     call print_matlu(greenlda%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1,opt_exp=2) ! debug
 690 !   end if
 691 !   call flush_unit(std_out)
 692  end if
 693 
 694 ! =================================================================
 695 ! Compute analytic limit of hybridization and rotate it
 696 ! =================================================================
 697  ABI_DATATYPE_ALLOCATE(hybri_coeff,(paw_dmft%natom))
 698  call init_matlu(paw_dmft%natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,hybri_coeff)
 699  !write(6,*)"hybri1",hybri_coeff(1)%mat(1,1,1,1,1),paw_dmft%natom,cryst_struc%natom
 700 
 701  ! Compute analytical C_ij such that F_ij -> C_ij/iw_n
 702  ! ---------------------------------------
 703  call hybridization_asymptotic_coefficient(cryst_struc,paw_dmft,pawang,hybri_coeff)
 704  write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n for large frequency"
 705  call wrtout(std_out,message,'COLL')
 706 
 707  ! Print analytical C_ij (not rotated)
 708  ! ---------------------------------------
 709  call print_matlu(hybri_coeff,natom,1)
 710 
 711  ! Rotate analytical C_ij in Ylm basis
 712  ! ---------------------------------------
 713  if(useylm==1) call slm2ylm_matlu(hybri_coeff,natom,1,pawprtvol)
 714  if(opt_diag/=0)  then
 715 
 716  ! Rotate analytical C_ij in rotated basis
 717  ! ---------------------------------------
 718    if(opt_rot==1.or.opt_rot==2) call rotate_matlu(hybri_coeff,eigvectmatlu,natom,3,1)
 719 
 720  ! Print analytical C_ij (rotated)
 721  ! ---------------------------------------
 722    write(message,'(a,2x,a)') ch10," == Coeff analytical C_ij such that F -> C_ij/iw_n after rotation"
 723    call wrtout(std_out,message,'COLL')
 724    call print_matlu(hybri_coeff,natom,1,compl=1,opt_exp=1)
 725  end if
 726 
 727 ! =================================================================
 728 ! Check if rotation is properly done.
 729 ! =================================================================
 730  if(3==4) then
 731    write(message,'(a,2x,a)') ch10,&
 732 &   " == Print  dmat before rot"
 733    call wrtout(std_out,message,'COLL')
 734    call print_matlu(green%occup%matlu,natom,1)
 735    if(useylm==1) call slm2ylm_matlu(green%occup%matlu,natom,1,pawprtvol)
 736    if(opt_rot==1) call rotate_matlu(green%occup%matlu,eigvectmatlu,natom,3,1)
 737    write(message,'(a,2x,a)') ch10,&
 738 &   " == Print  dmat after rot"
 739    call wrtout(std_out,message,'COLL')
 740    call print_matlu(green%occup%matlu,natom,1)
 741 
 742    write(message,'(2a)') ch10,' QMC STOP: DEBUG'
 743    call wrtout(std_out,message,'COLL')
 744    MSG_ERROR(message)
 745  end if
 746 ! =================================================================
 747 ! Check
 748 ! =================================================================
 749 
 750 ! write(message,'(a,2x,a,f13.5)') ch10,&
 751 !&   " == Print weiss for small tau"
 752 ! call wrtout(std_out,message,'COLL')
 753 ! call print_matlu(weiss%oper(1)%matlu,natom,1)
 754 ! write(message,'(a,2x,a,f13.5)') ch10,&
 755 !&   " == Print weiss for large tau"
 756 ! call wrtout(std_out,message,'COLL')
 757 ! call print_matlu(weiss%oper(paw_dmft%dmft_nwlo)%matlu,natom,1)
 758 ! call flush_unit(std_out)
 759 ! write(message,'(2a)') ch10,' Check weiss_for_rot(last freq)'
 760 ! call wrtout(std_out,message,'COLL')
 761 ! call checkdiag_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,tol6,opt=nspinor)
 762 ! call flush_unit(std_out)
 763 ! write(message,'(2a)') ch10,' Check weiss_for_rot(ifreq=1)'
 764 ! call wrtout(std_out,message,'COLL')
 765 ! call checkdiag_matlu(weiss_for_rot%oper(1)%matlu,natom,tol6,opt=nspinor)
 766 ! call flush_unit(std_out)
 767 
 768  master=0
 769 
 770 ! =================================================================
 771 ! Print out
 772 ! =================================================================
 773 
 774 ! Print Weiss
 775 ! -------------
 776  if(paw_dmft%dmft_prgn==1) then
 777    call print_green('Weiss_diag',weiss_for_rot,1,paw_dmft,pawprtvol=1,opt_wt=1,opt_decim=1)
 778  end if
 779 
 780  write(message,'(a,2x,a,f13.5)') ch10,&
 781 & " == Preparing data for CTQMC"
 782  call wrtout(std_out,message,'COLL')
 783 
 784 ! Print Rotate Weiss for 1st and last frequencies
 785 ! ------------------------------------------------
 786  if (pawprtvol>=3) then
 787    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 788 &  " == Print rotated weiss function for small freq in the rotated basis"  ! debug
 789    call wrtout(std_out,message,'COLL')  ! debug
 790    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1,compl=1)  ! debug
 791    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 792 &  " == Print rotated weiss function for largest freq in the rotated basis"  ! debug
 793    call wrtout(std_out,message,'COLL')  ! debug
 794    call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1)  ! debug
 795  end if
 796 
 797 ! =================================================================
 798 !  VARIABLES FOR CTQMC TESTS
 799  testcode = 0
 800  testrot  = 0
 801  opt_fk=0 ! for developpers to check Fourier transform and computes G0(tau)
 802  opt_fk=1 ! usual case: for real calculations
 803 ! =================================================================
 804 
 805 ! ___________________________________________________________________________________
 806 !
 807 !  SECOND PART : BUILT HYBRIDIZATION FROM G0
 808 ! ___________________________________________________________________________________
 809 !
 810 ! ===========================================================================================
 811 ! Compute inverse of weiss  and compute hybridization
 812 ! ===========================================================================================
 813 
 814 ! Compute inverse of weiss  for each Frequency
 815 ! ----------------------------------------------
 816  do ifreq=1,paw_dmft%dmft_nwlo
 817    ABI_DATATYPE_ALLOCATE(matlu1,(natom))
 818    ABI_DATATYPE_ALLOCATE(matlu2,(natom))
 819    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
 820    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2)
 821 
 822    call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,natom)
 823 
 824    ! Print G_0(iw_n)
 825    ! ----------------
 826    if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"go",60000,imre=1)
 827 
 828    ! Compute G_0^-1
 829    ! -------------------------------------------
 830    ! if opt_fk=1 or testcode/=0  Do the inversion
 831    ! if opt_fk=0                 Do not inverse.
 832    ! If testcode=2 and opt_fk=0  Do the inversion
 833    ! If testcode=1 and opt_fk=0  Do the inversion but no effect, because it will nevertheless be erased
 834    ! If opt_fk=1                 Do the inversion
 835    ! -------------------------------------------
 836    if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"weiss",12000,imre=1)
 837    if(opt_fk==1.or.testcode/=0) call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1)
 838 
 839    ! Print G_0^-1(iw_n)
 840    ! ----------------
 841    if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"goinv",70000,imre=1)
 842 
 843    if(pawprtvol>=4.or.ifreq==paw_dmft%dmft_nwlo) then
 844      if(opt_fk==1.or.testcode/=0) then
 845       ! Check inversion : do the product
 846       ! ----------------------------------------------
 847        call prod_matlu(weiss_for_rot%oper(ifreq)%matlu,matlu1,matlu2,natom)
 848        write(message,'(a,2x,a,i7)') ch10,&  ! debug
 849 &      " == Print product of  weiss times invers for freq",ifreq
 850        call wrtout(std_out,message,'COLL')  ! debug
 851        call print_matlu(matlu2,natom,1)  ! debug
 852      end if
 853    end if
 854 
 855    call destroy_matlu(matlu1,natom)
 856    call destroy_matlu(matlu2,natom)
 857    ABI_DATATYPE_DEALLOCATE(matlu1)
 858    ABI_DATATYPE_DEALLOCATE(matlu2)
 859  end do
 860 
 861  ! Copy weiss_for_rot into weiss
 862  ! -------------------------------
 863  !call copy_matlu(weiss_for_rot%oper(ifreq)%matlu,weiss%oper(ifreq)%matlu,natom)
 864 
 865 
 866  ! Print G_0^-1 for 1st and last frequencies.
 867  ! -----------------------------------------
 868  if(pawprtvol>=3) then
 869    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 870 &  " == Print G_0^-1 for small freq in the rotated basis"  ! debug
 871    call wrtout(std_out,message,'COLL')  ! debug
 872    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1)  ! debug
 873    write(message,'(a,2x,a,e18.10,a)') ch10,&   ! debug
 874 &  " == Print G_0^-1 for last freq in the rotated basis (last freq=", paw_dmft%omega_lo(paw_dmft%dmft_nwlo),")"  ! 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 ! Substract frequency from diagonal part
 880 ! ======================================
 881 
 882  ABI_ALLOCATE(shift,(natom))
 883  do ifreq=1,paw_dmft%dmft_nwlo
 884    shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
 885 
 886 !  write(5555,'(400e17.4)') paw_dmft%omega_lo(ifreq),((((((weiss_for_rot%oper(ifreq)%matlu(1)%mat&
 887 !  & (im,im1,isppol,ispinor,ispinor1)-cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)),im=1,2*3+1),&
 888 !&      im1=1,2*3+1),isppol=1,nsppol),ispinor=1,nspinor),ispinor1=1,nspinor)
 889 
 890   ! Compute G_0^-1-iw_n
 891   ! --------------------
 892    if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift)
 893 
 894   ! Compute -G_0^-1+iw_n
 895   ! --------------------
 896    if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone)
 897 
 898   ! Print -G_0^-1+iw_n
 899   ! --------------------
 900    if(optdb==1) then
 901      call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"G0inv_minus_omega",20000,imre=1)
 902    end if
 903  end do
 904 
 905  ! Print -G_0^+1-iw_n=(F-levels) for last freq in the rotated basis"
 906  ! ------------------------------------------------------------------
 907  ABI_DEALLOCATE(shift)
 908  if(pawprtvol>=3) then
 909    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
 910 &  " == Print G_0^-1-iw_n=-(F-levels) for last freq in the rotated basis"  ! debug
 911    call wrtout(std_out,message,'COLL')   ! debug
 912    call print_matlu(weiss_for_rot%oper(paw_dmft%dmft_nwlo)%matlu,natom,1,compl=1) ! debug
 913  end if
 914 
 915 ! Check numerical limit of F(i_wn)*iw_n (can be used also to compute F )
 916 ! ======================================
 917 
 918  if(opt_nondiag==1) then
 919    ABI_DATATYPE_ALLOCATE(matlu1,(natom))
 920    ABI_DATATYPE_ALLOCATE(matlu2,(natom))
 921    ABI_DATATYPE_ALLOCATE(matlu3,(natom))
 922    ABI_DATATYPE_ALLOCATE(matlu4,(natom))
 923    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
 924    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu2)
 925    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu3)
 926    call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu4)
 927 
 928    write(message,'(a,2x,a,f13.5)') ch10,  " == energy_levels"
 929    call wrtout(std_out,message,'COLL')
 930    call print_matlu(energy_level%matlu,natom,1,opt_exp=2,compl=1)
 931 
 932    do ifreq=paw_dmft%dmft_nwlo,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency.
 933    !do ifreq=paw_dmft%dmftqmc_l,1,-1 ! necessary to have matlu4 computed for the max frequency and available for all frequency.
 934       ! Compute F (substract levels) for max frequency
 935       ! -----------------------------------------------
 936      call add_matlu(weiss_for_rot%oper(ifreq)%matlu,energy_level%matlu,matlu1,natom,-1)
 937 
 938       ! Print F(iw_n)=-(G_0^-1-iw_n+levels)  for last frequency.
 939       ! --------------------------------------------------------
 940      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
 941        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
 942 &       " == Print F(iw_n)=-(G_0^-1-iw_n+levels) for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
 943        call wrtout(std_out,message,'COLL')
 944        call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 945      end if
 946      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"Hybridization",10000,imre=1)
 947 
 948       ! Put F in weiss_for_rot -> CTQMC
 949       ! -------------------------------
 950      if(opt_rot==2) call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,1)
 951 !   The following line will produce directly the weiss function for the CTQMC code
 952      if(opt_fk==1) call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom)
 953 
 954 
 955       ! Multiply F by frequency
 956       ! ------------------------
 957      call copy_matlu(matlu1,matlu2,natom)
 958      call fac_matlu(matlu1,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp))
 959      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
 960        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
 961 &       " == Print numerical C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
 962        call wrtout(std_out,message,'COLL')
 963        call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 964      end if
 965      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij",72800,imre=1)
 966      !call rotate_matlu(matlu1,eigvectmatlu,natom,3,1)
 967 
 968      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
 969        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
 970 &       " == Print numerical after back rotation C_ij = F(iw_n)*iw_n for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
 971        call wrtout(std_out,message,'COLL')
 972        call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 973      end if
 974      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_rotated",72900,imre=1)
 975 
 976       ! Built C_ij/iw_n
 977       ! ------------------------
 978      call copy_matlu(hybri_coeff,matlu1,natom)
 979      call fac_matlu(matlu1,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp))
 980      if(optdb==1) call printplot_matlu(matlu1,natom,paw_dmft%omega_lo(ifreq),"cij_over_omega",72000)
 981     ! if(ifreq==paw_dmft%dmft_nwlo) then
 982     !   write(message,'(a,2x,a,f13.5)') ch10,  " == Print numerical C_ij/iw_n for frequency",paw_dmft%omega_lo(ifreq)
 983     !   call wrtout(std_out,message,'COLL')
 984     !   call print_matlu(matlu1,natom,1,opt_exp=1,compl=1)
 985     ! endif
 986 
 987       ! For test: put C_ij/i_wn into weiss_for_rot
 988       ! --------------------------------------------
 989      !call copy_matlu(matlu1,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1)
 990 
 991       ! Compute Hybri - C_ij/iw_n
 992       ! ------------------------
 993      call add_matlu(matlu2,matlu1,matlu3,natom,-1)
 994 
 995       ! Print Hybri - C_ij/iw_n
 996       ! ------------------------
 997      if(optdb==1) call printplot_matlu(matlu3,natom,paw_dmft%omega_lo(ifreq),"hybri_minus_asymp",74000,imre=1)
 998 
 999       ! 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.
1000       ! ------------------------------------------------------------------------------------------------------------------
1001      call copy_matlu(matlu3,matlu2,natom)
1002      call fac_matlu(matlu2,natom,cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2)
1003      if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"fminuscijtimesw2",75000,imre=1)
1004      if(ifreq==paw_dmft%dmft_nwlo.or.ifreq==paw_dmft%dmftqmc_l) then
1005        call copy_matlu(matlu2,matlu4,natom)
1006        write(message,'(a,2x,a,i4,a,f13.5,a)') ch10, &
1007 &       " == Print numerical (F(iw_n)-C_ij/iw_n)%iw_n^2 for freq nb",ifreq," (=",paw_dmft%omega_lo(ifreq),")"
1008        call wrtout(std_out,message,'COLL')
1009        call print_matlu(matlu4,natom,1)
1010      end if
1011 
1012       ! Built C_ij/iw_n+D_ij/(iw_n)^2
1013       ! ------------------------
1014      call copy_matlu(matlu4,matlu3,natom,opt_re=1)
1015      call fac_matlu(matlu3,natom,1.d0/cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)**2)
1016      call add_matlu(matlu1,matlu3,matlu2,natom,1)
1017      if(optdb==1) call printplot_matlu(matlu2,natom,paw_dmft%omega_lo(ifreq),"cij_w_plus_dij_w2",72700,imre=1)
1018       ! For test: put C_ij/i_wn +D_ij/(iw_n)^2 into weiss_for_rot
1019       ! --------------------------------------------
1020      !call copy_matlu(matlu2,weiss_for_rot%oper(ifreq)%matlu,natom,opt_non_diag=1)
1021 
1022 
1023    end do
1024 
1025    call destroy_matlu(matlu1,natom)
1026    call destroy_matlu(matlu2,natom)
1027    call destroy_matlu(matlu3,natom)
1028    call destroy_matlu(matlu4,natom)
1029    ABI_DATATYPE_DEALLOCATE(matlu1)
1030    ABI_DATATYPE_DEALLOCATE(matlu2)
1031    ABI_DATATYPE_DEALLOCATE(matlu3)
1032    ABI_DATATYPE_DEALLOCATE(matlu4)
1033  end if
1034 
1035 ! =========================================================================================
1036 ! Start big loop to compute hybridization
1037 ! =========================================================================================
1038  do iatom=1,cryst_struc%natom
1039    green%ecorr_qmc(iatom)=zero
1040    itypat=cryst_struc%typat(iatom)
1041    lpawu=paw_dmft%lpawu(iatom)
1042    tndim=2*lpawu+1
1043    if(lpawu/=-1) then
1044 
1045      nflavor=2*(tndim)
1046      if(testcode>=1) then
1047        nflavor=2
1048        if(testcode==2) then
1049          ispa=1
1050          ispb=2
1051          if(nspinor==1) ispb=1
1052          ima=1
1053          imb=1
1054          if(tndim>4) then
1055            ima=5 ! row
1056            imb=4 ! column
1057          end if
1058        end if
1059      end if
1060 
1061      !allocate(correl_loc(nflavor,nflavor))
1062      !ABI_ALLOCATE(f_with_k,(MIN(paw_dmft%dmft_nwli,paw_dmft%dmftqmc_l),nflavor))
1063      ABI_ALLOCATE(fw1,(paw_dmft%dmft_nwlo,nflavor))
1064      ABI_ALLOCATE(fw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor))
1065      ABI_ALLOCATE(levels_ctqmc,(nflavor))
1066      ABI_ALLOCATE(levels_ctqmc_nd,(nflavor,nflavor))
1067      levels_ctqmc_nd=czero
1068      ABI_ALLOCATE(hybri_limit,(nflavor,nflavor))
1069      hybri_limit=czero
1070      fw1_nd=czero
1071      fw1=czero
1072      !allocate(fw2(paw_dmft%dmft_nwli))
1073 ! =================================================================
1074 ! Compute Hybridization
1075 ! =================================================================
1076 
1077      if (testcode==0) then
1078        iflavor1=0
1079        iflavor2=0
1080        do isppol=1,nsppol
1081          do ispinor1=1,nspinor
1082            do ispinor2=1,nspinor
1083              do im1=1,tndim
1084                do im2=1,tndim
1085                  ! first diagonal terms whatever opt_nondiag
1086                  iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
1087                  iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
1088                  !write(6,*) isppol,ispinor1,ispinor2,im1,im2
1089                  !write(6,*) iflavor1,iflavor2
1090 
1091                  if ( iflavor1==iflavor2 ) then
1092 
1093 !              Do spline of weiss function for im and isppol
1094 !              Construction of fw1
1095 !             =================================================================
1096 
1097                    do ifreq=1,paw_dmft%dmft_nwlo
1098                      if(opt_fk==1) then
1099                        fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1)
1100                      else if (opt_fk==0) then
1101                        fw1(ifreq,iflavor1)= weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1)
1102                      end if
1103                    end do
1104                    fw1_nd(:,iflavor1,iflavor1)=fw1(:,iflavor1)
1105 
1106                    levels_ctqmc(iflavor1)=real(energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1),kind=dp)
1107                    hybri_limit(iflavor1,iflavor1)=hybri_coeff(iatom)%mat(im1,im1,isppol,ispinor1,ispinor1)
1108 
1109 
1110                    if(nsppol==1.and.nspinor==1) then
1111                      !f_with_k(:,iflavor+tndim)=f_with_k(:,iflavor)
1112                      fw1(:,iflavor1+tndim)=fw1(:,iflavor1)
1113                      fw1_nd(:,iflavor1+tndim,iflavor1+tndim)=fw1(:,iflavor1)
1114                      levels_ctqmc(iflavor1+tndim)=levels_ctqmc(iflavor1)
1115                      hybri_limit(iflavor1+tndim,iflavor1+tndim)=hybri_limit(iflavor1,iflavor1)
1116                    end if
1117 
1118                  else
1119 
1120                    do ifreq=1,paw_dmft%dmft_nwlo
1121                      if(opt_fk==1) then
1122                        fw1_nd(ifreq,iflavor1,iflavor2)= &
1123 &                       weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1124                      else if (opt_fk==0) then
1125                        fw1_nd(ifreq,iflavor1,iflavor2)= &
1126 &                       weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1127                      end if
1128 
1129         ! omega=pi*paw_dmft%temp*(two*float(ifreq)-1)
1130          !write(3333,*) omega,weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1131          !write(4444,*) omega,fw1_nd(ifreq,iflavor1,iflavor2),"#",iflavor1,iflavor2
1132         ! if(iflavor1/=iflavor2)write(5555,*) omega,imag(fw1_nd(ifreq,iflavor1,iflavor2)),"#",iflavor1,iflavor2
1133 
1134                    end do
1135                    hybri_limit(iflavor1,iflavor2)=hybri_coeff(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1136 
1137         ! write(3333,*)
1138         ! write(4444,*)
1139         ! write(5555,*)
1140 
1141                    if(nsppol==1.and.nspinor==1) then
1142                      fw1_nd(:,iflavor1+tndim,iflavor2+tndim) = fw1_nd(:,iflavor1,iflavor2)
1143                      hybri_limit(iflavor1+tndim,iflavor2+tndim)=hybri_limit(iflavor1,iflavor2)
1144                    end if
1145 
1146                  end if
1147 
1148 ! <  / HACK >
1149                end do !im2
1150              end do !im1
1151            end do  !ispinor2
1152          end do  !ispinor1
1153        end do  !isppol
1154 ! < HACK >
1155        ! JB. On 1000 cpus this can not work since all CPU try to open/write the files
1156        ! Action : Don't print it or check only one cpu does it.
1157        if(pawprtvol>=10000000) then
1158          write(message,'(a,2x,a)') ch10,  " == Hybri for all flavors for CTQMC "
1159          call wrtout(std_out,message,'COLL')
1160          do iflavor1=1,nflavor
1161            write(message,'(4x,14(2e14.5,2x))') (hybri_limit(iflavor1,iflavor2),iflavor2=1,nflavor)
1162            call wrtout(std_out,message,'COLL')
1163          end do
1164 
1165          if (open_file('Hybri_cijoveromega',message, newunit=unt, status='unknown', form='formatted') /= 0) then
1166            MSG_ERROR(message)
1167          end if
1168          if (open_file('Hybri',message,newunit=unt2,status='unknown',form='formatted') /= 0) then
1169            MSG_ERROR(message)
1170          end if
1171          do ifreq=1,paw_dmft%dmft_nwlo
1172 !              weiss_for_rot is G_0^-1-iw_n=-(F-levels)
1173            if(optdb==1) call printplot_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,paw_dmft%omega_lo(ifreq),"weissbefore112",30000)
1174          end do
1175          do iflavor1=1,nflavor
1176            do iflavor2=1,nflavor
1177              do ifreq=1,paw_dmft%dmft_nwlo
1178                omega=pi*paw_dmft%temp*(two*float(ifreq)-1)
1179 !              fw1_nd is -G_0^+1-iw_n=(F-levels)
1180                write(unt,'(300e16.5)') paw_dmft%omega_lo(ifreq)&
1181 &               ,fw1_nd(ifreq,iflavor1,iflavor2)-hybri_limit(iflavor1,iflavor2)/cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)
1182                write(unt2,'(300e16.5)') paw_dmft%omega_lo(ifreq),fw1_nd(ifreq,iflavor1,iflavor2)
1183                !write(1111,*) omega,real(fw1_nd(ifreq,iflavor1,iflavor2))
1184                !write(1112,*) omega,imag(fw1_nd(ifreq,iflavor1,iflavor2))
1185              end do
1186              write(unt,*)
1187              write(unt2,*)
1188             ! write(1111,*)
1189             ! write(1112,*)
1190            end do
1191          end do
1192          close(unt)
1193          close(unt2)
1194        end if
1195      end if ! testcode
1196    ! </ HACK >
1197 
1198 ! ====================================================================================
1199 !  TEST
1200 !  For testing purpose, built ultra simple hybridization (constant in
1201 !  imaginary time or very simple) or extract some part of the calculated hybridization
1202 ! ====================================================================================
1203 
1204      if(testcode>=1) then
1205        dmft_nwlo=paw_dmft%dmft_nwlo
1206        paw_dmft%dmft_nwlo=paw_dmft%dmftqmc_l
1207        ABI_ALLOCATE(gw1_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor))
1208        gw1_nd=czero
1209 
1210        !  Call testcode_ctqmc: built simple hybridization
1211        !--------------------------------------------------
1212        if (testcode==1) then
1213          call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,&
1214 &         levels_ctqmc,hybri_limit,nflavor,1,paw_dmft%temp,testrot,testcode,umod)
1215        !  Select 2x2 hybridization matrix from the current larger matrix
1216        !  ima and imb are defined above.
1217        !----------------------------------------------------------------
1218        else if (testcode==2) then
1219          facnd=0.8d0
1220          facd=1.0d0
1221          !write(6,*) "fac",facnd,facd
1222          levels_ctqmc_nd(2,2)   = energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb)
1223          levels_ctqmc_nd(1,1)   = energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa)
1224          levels_ctqmc(2)   = real(energy_level%matlu(iatom)%mat(imb,imb,1,ispb,ispb),kind=dp)
1225          levels_ctqmc(1)   = real(energy_level%matlu(iatom)%mat(ima,ima,1,ispa,ispa),kind=dp)
1226          if(opt_diag/=1) then
1227            levels_ctqmc_nd(1,2)   = energy_level%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1228            levels_ctqmc_nd(2,1)   = energy_level%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1229          end if
1230          hybri_limit(1,1)  = facd*hybri_coeff(iatom)%mat(ima,ima,1,ispa,ispa)
1231          hybri_limit(2,2)  = facd*hybri_coeff(iatom)%mat(imb,imb,1,ispb,ispb)
1232          hybri_limit(1,2)  = facnd*hybri_coeff(iatom)%mat(ima,imb,1,ispa,ispb)
1233          hybri_limit(2,1)  = facnd*hybri_coeff(iatom)%mat(imb,ima,1,ispb,ispa)
1234          !write(6,*) "hybri_limit",hybri_limit
1235          !write(6,*) "levels_ctqmc",levels_ctqmc
1236          umod=zero
1237 
1238          tmpfil = 'fw1_nd_re'
1239          !if (open_file(newunit=unt,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then
1240          !  MSG_ERROR(message)
1241          !end if
1242          tmpfil = 'fw1_nd_im'
1243          !if (open_file(newunit=unt2,message,file=trim(tmpfil),status='unknown',form='formatted')/=0) then
1244          !  MSG_ERROR(message)
1245          !end if
1246          write(std_out,*) "testcode==2",ispa,ispb,ima,imb
1247          write(std_out,*) "opt_fk==",opt_fk
1248          do ifreq=1,paw_dmft%dmftqmc_l
1249            if (opt_fk==1) then
1250              fw1_nd(ifreq,1,1) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa)
1251              fw1_nd(ifreq,2,2) = facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb)
1252              !fw1_nd(ifreq,1,2) =  weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1253              !fw1_nd(ifreq,2,1) =  weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1254              fw1_nd(ifreq,1,2) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1255              fw1_nd(ifreq,2,1) = facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1256              omega=pi*paw_dmft%temp*(two*float(ifreq)-1)
1257            else if (opt_fk==0) then
1258              fw1_nd(ifreq,1,1) =  facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,ima,1,ispa,ispa)
1259              fw1_nd(ifreq,2,2) =  facd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,imb,1,ispb,ispb)
1260              fw1_nd(ifreq,1,2) =  facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(ima,imb,1,ispa,ispb)
1261              fw1_nd(ifreq,2,1) =  facnd*weiss_for_rot%oper(ifreq)%matlu(iatom)%mat(imb,ima,1,ispb,ispa)
1262              call xginv(fw1_nd(ifreq,:,:),2)
1263            end if
1264          end do
1265          !close(unt)
1266          !close(unt2)
1267        end if
1268 
1269        ! Calculation of Inverse Green's function from hybridization
1270        !-------------------------------------------------------------
1271        do if1=1,2
1272          do if2=1,2
1273            do ifreq=1,paw_dmft%dmftqmc_l
1274              omega=pi*paw_dmft%temp*(two*float(ifreq)-1)
1275              if(if1==if2) then
1276                gw1_nd(ifreq,if1,if2) =  (cmplx(0.d0,omega,kind=dp)-fw1_nd(ifreq,if1,if2))
1277              else
1278                gw1_nd(ifreq,if1,if2) =  (-fw1_nd(ifreq,if1,if2))
1279              end if
1280            end do
1281          end do
1282        end do
1283        ! Calculation of Green's function (call inverse)
1284        !-------------------------------------------------------------
1285        do ifreq=1,paw_dmft%dmftqmc_l
1286          call xginv(gw1_nd(ifreq,:,:),2)
1287        end do
1288        write(std_out,*) " testctqmc high frequency limit of hybridization",fw1_nd(paw_dmft%dmftqmc_l,:,:)
1289 
1290        ! Integrate Green's function
1291        !-------------------------------------------------------------
1292        do if1=1,2
1293          do if2=1,2
1294            do ifreq=1,paw_dmft%dmftqmc_l
1295              omega=pi*paw_dmft%temp*(two*float(ifreq)-1)
1296              !write(999,*) omega,gw1_nd(ifreq,if1,if2)
1297            end do
1298            !write(999,*)
1299            call int_fct(gw1_nd(:,if1,if2),(if1==if2),2,paw_dmft,integral(if1,if2))  ! test_1
1300           ! write(std_out,*) "testctqmc occupations of input Green's function",(integral(if1,if2))
1301          end do
1302        end do
1303        ! Write Occupations
1304        !-------------------------------------------------------------
1305        do if1=1,2
1306          do if2=1,2
1307           ! write(std_out,*) "testctqmc occupations of input Green's function",(integral(if1,if2)+conjg(integral(if2,if1)))/two
1308          end do
1309        end do
1310        write(std_out,*) "Occupation of model in matrix form"
1311        do if1=1,2
1312          write(std_out,'(2(2f13.5,3x))') ((integral(if1,if2)+conjg(integral(if2,if1)))/two,if2=1,2)
1313        end do
1314        write(std_out,*) "Limit of hybridization "
1315        do if1=1,2
1316          write(std_out,'(2(2f13.5,3x))') (hybri_limit(if1,if2),if2=1,2)
1317        end do
1318 
1319        ! If opt_fk=0, give Green's function to CTQMC code instead of
1320        ! hybridization
1321        !-------------------------------------------------------------
1322        if(opt_fk==0) then
1323          fw1_nd=gw1_nd
1324        end if
1325 
1326        ABI_DEALLOCATE(gw1_nd)
1327        paw_dmft%dmft_nwlo=dmft_nwlo
1328      end if
1329 
1330      call flush_unit(std_out)
1331 ! =================================================================
1332 
1333 ! ___________________________________________________________________________________
1334 !
1335 !  THIRD PART : CALL CTQMC
1336 ! ___________________________________________________________________________________
1337 
1338 ! =================================================================
1339 !    Main calls to CTQMC code in ABINIT
1340 ! =================================================================
1341      if(paw_dmft%dmft_solv==5.or.paw_dmft%dmft_solv==8) then
1342 
1343        write(message,'(a,2x,a)') ch10,&
1344 &       " == Initializing CTQMC"
1345        call wrtout(std_out,message,'COLL')
1346 
1347 !    Initialisation
1348 ! =================================================================
1349      if(paw_dmft%dmft_solv==5) then
1350        nomega=paw_dmft%dmftqmc_l
1351        call CtqmcInterface_init(hybrid,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, &
1352 &       paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,&
1353 &       std_out,paw_dmft%spacecomm)
1354 !    options
1355 ! =================================================================
1356        call CtqmcInterface_setOpts(hybrid,&
1357        opt_Fk      =opt_fk,&
1358 &       opt_order   =paw_dmft%dmftctqmc_order ,&
1359 &       opt_movie   =paw_dmft%dmftctqmc_mov   ,&
1360 &       opt_analysis=paw_dmft%dmftctqmc_correl,&
1361 &       opt_check   =paw_dmft%dmftctqmc_check ,&
1362 &       opt_noise   =paw_dmft%dmftctqmc_grnns ,&
1363 &       opt_spectra =paw_dmft%dmftctqmc_mrka  ,&
1364 &       opt_gmove   =paw_dmft%dmftctqmc_gmove )
1365      endif
1366 
1367      if(paw_dmft%dmft_solv==8) then
1368        nomega=paw_dmft%dmftqmc_l
1369        call CtqmcoffdiagInterface_init(hybridoffdiag,paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, &
1370 &        paw_dmft%dmftqmc_therm, paw_dmft%dmftctqmc_meas,nflavor,&
1371 &        paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,&
1372 &        std_out,paw_dmft%spacecomm,opt_nondiag)
1373 !    options
1374 ! =================================================================
1375        call CtqmcoffdiagInterface_setOpts(hybridoffdiag,&
1376        opt_Fk      =opt_fk,&
1377 &       opt_order   =paw_dmft%dmftctqmc_order ,&
1378 &       opt_movie   =paw_dmft%dmftctqmc_mov   ,&
1379 &       opt_analysis=paw_dmft%dmftctqmc_correl,&
1380 &       opt_check   =paw_dmft%dmftctqmc_check ,&
1381 &       opt_noise   =paw_dmft%dmftctqmc_grnns ,&
1382 &       opt_spectra =paw_dmft%dmftctqmc_mrka  ,&
1383 &       opt_gmove   =paw_dmft%dmftctqmc_gmove )
1384      endif
1385 
1386        write(message,'(a,2x,2a)') ch10,&
1387 &       " == Initialization CTQMC done", ch10
1388        call wrtout(std_out,message,'COLL')
1389      end if
1390 
1391      if(paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7) then
1392        ABI_ALLOCATE(gw_tmp_nd,(paw_dmft%dmft_nwli,nflavor,nflavor)) !because size allocation problem with TRIQS paw_dmft%dmft_nwlo must be >= paw_dmft%dmft_nwli
1393          open(unit=505,file=trim(paw_dmft%filapp)//"_Legendre_coefficients.dat", status='unknown',form='formatted')
1394      else
1395        if(paw_dmft%dmft_solv==5) ABI_ALLOCATE(gw_tmp,(paw_dmft%dmft_nwlo,nflavor+1))
1396        ABI_ALLOCATE(gw_tmp_nd,(paw_dmft%dmft_nwlo,nflavor,nflavor+1))
1397 !     use  gw_tmp to put freq
1398        do ifreq=1,paw_dmft%dmft_nwlo
1399           if(paw_dmft%dmft_solv==5) gw_tmp(ifreq,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
1400          gw_tmp_nd(ifreq,nflavor,nflavor+1)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
1401        end do
1402      end if
1403 
1404      ABI_ALLOCATE(gtmp,(paw_dmft%dmftqmc_l,nflavor))
1405      ! THIS IS A BACKUP PLAN. USING paw_dmft%hybrid makes a segfault on TIKAL
1406      ! PSC with MPI only (and max2_open64). paw_dmf%hybrid is corrupted
1407      ! somewhere but I could not find the place in all DMFT routines
1408      ABI_ALLOCATE(gtmp_nd,(paw_dmft%dmftqmc_l,nflavor,nflavor))
1409      call flush_unit(std_out)
1410 
1411      if(testcode==0) then
1412       !unitnb=100000+paw_dmft%myproc
1413       !call int2char4(paw_dmft%myproc,tag_proc)
1414       !tmpfil = 'hybrilimit'//tag_proc
1415       !open (unit=unitnb,file=trim(tmpfil),status='unknown',form='formatted')
1416       !do iflavor1=1,nflavor
1417       !  do iflavor2=1,nflavor
1418       !    write(unitnb,*) iflavor1,iflavor2,hybri_limit(iflavor1,iflavor2)
1419       !  enddo
1420       !enddo
1421 
1422 ! =================================================================
1423 !    CTQMC run Abinit
1424 ! =================================================================
1425        if(paw_dmft%dmft_solv==5) then
1426 
1427          ABI_ALLOCATE(docc,(1:nflavor,1:nflavor))
1428          docc(:,:) = zero
1429          call CtqmcInterface_run(hybrid,fw1(1:paw_dmft%dmftqmc_l,:),Gtau=gtmp,&
1430 &         Gw=gw_tmp,D=docc(:,:),E=green%ecorr_qmc(iatom),&
1431 !&       matU=hu(itypat)%udens,opt_levels=levels_ctqmc)
1432 &         matU=udens_atoms(iatom)%value,opt_levels=levels_ctqmc)
1433          call data4entropyDMFT_setDocc(paw_dmft%forentropyDMFT,iatom,docc)
1434          ABI_DEALLOCATE(docc)
1435 
1436 ! =================================================================
1437 !    CTQMC run TRIQS
1438 ! =================================================================
1439        else if (paw_dmft%dmft_solv==8) then
1440          ABI_ALLOCATE(docc,(1:nflavor,1:nflavor))
1441          docc(:,:) = zero
1442          call CtqmcoffdiagInterface_run(hybridoffdiag,fw1_nd(1:paw_dmft%dmftqmc_l,:,:),Gtau=gtmp_nd,&
1443 &        Gw=gw_tmp_nd,D=doccsum,E=green%ecorr_qmc(iatom),&
1444 &        Noise=noise,matU=udens_atoms(iatom)%value,Docc=docc,opt_levels=levels_ctqmc,hybri_limit=hybri_limit)
1445          ABI_DEALLOCATE(docc)
1446        ! TODO: Handle de luj0 case for entropy
1447 
1448        else if (paw_dmft%dmft_solv>=6.and.paw_dmft%dmft_solv<=7) then
1449 
1450        ! fw1_nd: Hybridation
1451        ! levels_ctqmc: niveaux
1452        ! hu(itypat)%udens(:,:) : U_ij
1453        ! hu(itypat)%u(:,:,:,:) : uijkl
1454        ! temperature : paw_dmft%temp
1455        ! paw_dmft%dmftqmc_l: nombre de points en temps -1
1456        ! paw_dmft%dmftqmc_n: nombre de cycles
1457        ! ?? Quelles sorties: Les fonctions de Green
1458        ! frequence/temps/Legendre.
1459        ! Double occupations ?? <n_i n_j>
1460        ! test n_tau > 2*nfreq => ntau = 2*nfreq + 1
1461 !        for non diagonal code:
1462 !        call CtqmcInterface_run(hybrid,fw1_nd(1:paw_dmft%dmftqmc_l,:,:),Gtau=gtmp_nd,&
1463 !&       Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),&
1464 !&       Noise=Noise,matU=hu(itypat)%udens,opt_levels=levels_ctqmc,hybri_limit=hybri_limit)
1465        !Check choice of user to fix model bool var for the solver
1466          if (paw_dmft%dmft_solv==6) then
1467            rot_inv = .false.
1468          else !obviously paw_dmft%dmft_solv==7 with rot invariant terms
1469            rot_inv = .true.
1470          end if
1471 
1472          nfreq = paw_dmft%dmft_nwli
1473        !paw_dmft%dmft_nwlo = paw_dmft%dmft_nwli !transparent for user
1474          ntau  = paw_dmft%dmftqmc_l !(2*paw_dmft%dmftqmc_l)+1 !nfreq=paw_dmft%dmft_nwli
1475          nleg  = paw_dmft%dmftctqmc_triqs_nleg
1476 
1477          if ( ntau >= (2*nfreq)+1 ) then
1478 
1479            verbosity_solver = paw_dmft%prtvol
1480            beta = 1.0/(paw_dmft%temp*Ha_eV)
1481 
1482          !Allocation in/output array phase:
1483            ABI_ALLOCATE(fw1_nd_tmp,(1:nflavor,1:nflavor,1:nfreq)) !column major
1484            ABI_ALLOCATE(g_iw,(1:nflavor,1:nflavor,1:nfreq)) !column major
1485 !          ABI_ALLOCATE(gtau,(1:nflavor,1:nflavor,1:ntau)) !column major
1486            ABI_ALLOCATE(u_mat_ij,(1:nflavor,1:nflavor)) !column major
1487            ABI_ALLOCATE(u_mat_ijkl,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major
1488            ABI_ALLOCATE(u_mat_ijkl_tmp,(1:nflavor,1:nflavor,1:nflavor,1:nflavor)) !column major
1489 
1490            if ( leg_measure ) then !only if functionality is enabled
1491              ABI_ALLOCATE(gl_nd,(1:nleg,1:nflavor,1:nflavor)) !column major !nl = 30 by default
1492            end if
1493 
1494        !Conversion datas Ha -> eV (some duplications for test...)
1495        !fw1_nd_tmp = fw1_nd(1:paw_dmft%dmftqmc_l,:,:) * Ha_eV !fw1_nd = fw1_nd * Ha_eV !Ok?
1496 
1497            do iflavor=1,nflavor
1498              do iflavor1=1,nflavor
1499                do ifreq=1,nfreq
1500                  fw1_nd_tmp(iflavor,iflavor1,ifreq) = fw1_nd(ifreq,iflavor,iflavor1) * Ha_eV
1501 !              WRITE(500,*) "[IN Fortran] F[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",fw1_nd(ifreq,iflavor,iflavor1)
1502                end do
1503              end do
1504            end do
1505 
1506          !Report test
1507 !         WRITE(502,*) hu(itypat)%udens
1508 !         do ifreq=1,paw_dmft%dmftqmc_l
1509 !           write(501,*) ((fw1_nd(ifreq,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor)
1510 !         enddo
1511          !write(866,*)paw_dmft%dmft_nwlo,paw_dmft%dmftqmc_l
1512          !write(866,*) u_mat_ij
1513 !        do iflavor=1,nflavor+1
1514 !          do iflavor1=1,nflavor+1
1515 !            WRITE(502,*) "[OUT Fortran] U(i,j)[ l= ",iflavor," l_= ",iflavor1,"] = ",hu(itypat)%udens(iflavor,iflavor1)
1516 !          enddo
1517 !        enddo
1518 
1519 !          if(paw_dmft%myproc==0) then
1520 !          do iflavor=1,nflavor
1521 !            do iflavor1=1,nflavor
1522 !               do iflavor2=1,nflavor
1523 !                  do iflavor3=1,nflavor
1524 !                    write(490,*), hu(itypat)%vee(iflavor,iflavor1,iflavor2,iflavor3)
1525 !                  enddo
1526 !                 enddo
1527 !                enddo
1528 !              enddo
1529 !          endif
1530 
1531 !          if(paw_dmft%myproc==0) then
1532 !          do iflavor=1,nflavor
1533 !            do iflavor1=1,nflavor
1534 !            write(491,*), hu(itypat)%udens(iflavor,iflavor1) !(1,1,1,1)
1535 !            enddo
1536 !          enddo
1537 !          endif
1538 
1539 !          do iflavor=1,nflavor
1540 !            do iflavor1=1,nflavor
1541 !          do iflavor2=1,nflavor
1542 !            do iflavor3=1,nflavor
1543                  ! WRITE(552,*), hu(itypat)%vee!(iflavor,iflavor1,iflavor2,iflavor3)
1544 !            enddo
1545 !          enddo
1546 !            enddo
1547 !          enddo
1548 
1549            call vee_ndim2tndim_hu_r(lpawu,hu(itypat)%vee,u_mat_ijkl_tmp,1)
1550            do iflavor=1,nflavor
1551              do iflavor1=1,nflavor
1552                do iflavor2=1,nflavor
1553                  do iflavor3=1,nflavor
1554                    u_mat_ijkl(iflavor,iflavor1,iflavor2,iflavor3)   =  Ha_eV * u_mat_ijkl_tmp(iflavor,iflavor1,iflavor2,iflavor3)
1555                  end do
1556                end do
1557              end do
1558            end do
1559 
1560          !u_mat_ijkl   =  Ha_eV * reshape( u_mat_ijkl , [nflavor,nflavor,nflavor,nflavor] )  !column -> row major + conversion
1561            u_mat_ij     = transpose( hu(itypat)%udens ) * Ha_eV !column -> row major + conversion
1562            levels_ctqmc = levels_ctqmc * Ha_eV
1563 
1564          !Location array in memory for C++ pointer args to pass
1565            g_iw_ptr       = C_LOC( gw_tmp_nd ) !C_LOC( g_iw )
1566            gtau_ptr       = C_LOC( gtmp_nd ) !C_LOC( gtau )
1567            gl_ptr         = C_LOC( gl_nd )
1568            fw1_nd_ptr     = C_LOC( fw1_nd_tmp )
1569            u_mat_ij_ptr   = C_LOC( u_mat_ij )
1570            u_mat_ijkl_ptr = C_LOC( u_mat_ijkl )
1571            levels_ptr     = C_LOC( levels_ctqmc )
1572 
1573          !Calling interfaced TRIQS solver subroutine from src/01_triqs_ext package
1574 #if defined HAVE_TRIQS_v2_0 || defined HAVE_TRIQS_v1_4
1575            call Ctqmc_triqs_run (     rot_inv, leg_measure, hist, wrt_files, tot_not,                            &
1576 &           nflavor, nfreq, ntau , nleg, int(paw_dmft%dmftqmc_n/paw_dmft%nproc),       &
1577 &           paw_dmft%dmftctqmc_meas*2*2*nflavor, paw_dmft%dmftqmc_therm,               &
1578 &           verbosity_solver, paw_dmft%dmftqmc_seed,beta,                              &
1579 &           levels_ptr,  u_mat_ij_ptr, u_mat_ijkl_ptr, fw1_nd_ptr,                     &
1580 &           g_iw_ptr, gtau_ptr, gl_ptr, paw_dmft%spacecomm                             )
1581 #endif
1582 
1583          !WRITE(*,*) "Hello Debug"
1584          !call xmpi_barrier(paw_dmft%spacecomm) !Resynch all processus after calling Impurity solver from TRIQS
1585 
1586          !Report output datas from TRIQS to Abinit
1587          !Interacting G(iw)
1588            do ifreq=1,nfreq
1589              do iflavor1=1,nflavor
1590                do iflavor=1,nflavor
1591 !                gw_tmp_nd(ifreq,iflavor,iflavor1) = g_iw(iflavor,iflavor1,ifreq) !* Ha_eV !because 1/ G0(eV)
1592             !  WRITE(503,*) "[OUT Fortran] G(iw)[ w= ",ifreq," l= ",iflavor," l_= ",iflavor1,"] = ",gw_tmp_nd(ifreq,iflavor,iflavor1)!g_iw(iflavor,iflavor1,ifreq)
1593                end do
1594              end do
1595            end do
1596 
1597 !        Convert in Ha
1598            gw_tmp_nd = gw_tmp_nd*Ha_eV
1599 
1600 !            do iflavor1=1,nflavor
1601 !              do iflavor=1,nflavor
1602 !
1603 !               WRITE(510,*) "[OUT Fortran] U[ l= ",iflavor," l_= ",iflavor1,"] = ",u_mat_ij(iflavor,iflavor1)
1604 !              enddo
1605 !            enddo
1606 
1607 !        if(paw_dmft%myproc==0) write(6,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1)
1608 !        if(paw_dmft%myproc==1) write(6,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1)
1609 !        if(paw_dmft%myproc==0) write(621,*) "essai",paw_dmft%myproc, gw_tmp_nd(2,1,1)
1610 !        if(paw_dmft%myproc==1) write(622,*) "essai",paw_dmft%myproc,gw_tmp_nd(2,1,1)
1611 !        call flush_unit(621)
1612 !        call flush_unit(622)
1613 !        write(message,*) ch10, "essai",paw_dmft%myproc, paw_dmft%myproc,paw_dmft%dmftqmc_seed!gw_tmp_nd(2,1,1)
1614 !        call wrtout(555,message,'PERS',.true.)
1615 !        if(paw_dmft%myproc==0) write(499,*) "essai",paw_dmft%myproc, paw_dmft%dmftqmc_seed
1616 !        if(paw_dmft%myproc==1) write(498,*) "essai",paw_dmft%myproc,paw_dmft%dmftqmc_seed
1617 
1618          !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)...
1619 !          do iflavor=1,nflavor
1620 !            do iflavor1=1,nflavor
1621 !              do itau=1,ntau
1622 !                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
1623                  !gtmp_nd(itau,iflavor,iflavor1) = gtau(iflavor,iflavor1,itau)
1624 !                endif
1625 !                if(paw_dmft%myproc==0) then
1626 !                  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)
1627 !                endif
1628 !              enddo
1629 !            enddo
1630 !          enddo
1631 
1632          !Legendre Polynoms G(L) for extrapolation of Interacting G(iw) by FT, only if leg_measure == TRUE
1633            if (leg_measure) then
1634              do ileg=1,nleg
1635                WRITE(505,*) ileg,((gl_nd(ileg,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor)
1636              end do
1637              close(505)
1638            end if
1639 !        if(paw_dmft%myproc==0) then
1640 !          do itau=1,paw_dmft%dmftqmc_l
1641 !            write(490,*) ((gtmp_nd(itau,iflavor,iflavor1),iflavor=1,nflavor),iflavor1=1,nflavor)
1642 !          enddo
1643 !        endif
1644            ABI_DEALLOCATE( fw1_nd_tmp )
1645            ABI_DEALLOCATE( g_iw )
1646            ABI_DEALLOCATE( u_mat_ijkl )
1647            ABI_DEALLOCATE( u_mat_ijkl_tmp )
1648            ABI_DEALLOCATE( u_mat_ij )
1649 ! =========================================================================
1650 !   Compute Green's function in imaginary freq using Legendre coefficients
1651 ! =========================================================================
1652            if (leg_measure) then
1653              call xmpi_barrier(paw_dmft%spacecomm)
1654              call flush_unit(std_out)
1655              write(message,'(2a)') ch10,"    ==  Compute G(iw_n) from Legendre coefficients"
1656              call wrtout(std_out,message,'COLL')
1657              ABI_ALLOCATE( jbes, (nleg))
1658              gw_tmp_nd=czero
1659 
1660         !   write(77,*) " TEST OF BESSEL S ROUTINES 0 0"
1661 
1662         !   xx=0_dp
1663         !   ileg=0
1664         !   call sbf8(ileg+1,xx,jbes)
1665         !   write(77,*) "T0 A",jbes(ileg+1)
1666         !   call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx)
1667         !   write(77,*) "T0 B",jbes(ileg+1)
1668         !   write(77,*) "T0 C",bessel_jn(ileg,xx)
1669 
1670         !   write(77,*) " TEST OF BESSEL S ROUTINES 1.5 0"
1671 
1672         !   xx=1.5_dp
1673         !   ileg=0
1674         !   call sbf8(ileg+1,xx,jbes)
1675         !   write(77,*) "T1 A",jbes(ileg+1)
1676         !   call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx)
1677         !   write(77,*) "T1 B",jbes(ileg+1)
1678         !   write(77,*) "T1 C",bessel_jn(ileg,xx)
1679 
1680         !   write(77,*) " TEST OF BESSEL S ROUTINES 1.5 1"
1681 
1682         !   xx=1.5_dp
1683         !   ileg=1
1684         !   call sbf8(ileg+1,xx,jbes)
1685         !   write(77,*) "T2 A",jbes(ileg+1)
1686         !   call jbessel(jbes(ileg+1),besp,bespp,ileg,1,xx)
1687         !   write(77,*) "T2 B",jbes(ileg+1)
1688         !   write(77,*) "T2 C",bessel_jn(ileg,xx)
1689 
1690 
1691              do ifreq=1,paw_dmft%dmft_nwli
1692                xx=real(2*ifreq-1,kind=dp)*pi/two
1693                if(xx<=100_dp) call sbf8(nleg,xx,jbes)
1694                do ileg=1,nleg
1695               ! write(77,*) "A",ifreq,jbes(ileg),xx
1696 
1697                  if(xx>=99) call jbessel(jbes(ileg),besp,bespp,ileg-1,1,xx)
1698               ! write(77,*) "B",ifreq,jbes(ileg),xx
1699 
1700                !write(77,*) "C",ifreq,jbes(ileg),xx
1701 
1702                  u_nl=sqrt(float(2*ileg-1))*(-1)**(ifreq-1)*cmplx(0_dp,one)**(ileg)*jbes(ileg)
1703            !    write(77,*) "----------",ileg,jbes(ileg), u_nl,gl_nd(ileg,1,1)
1704 
1705                  do iflavor=1,nflavor
1706                    do iflavor1=1,nflavor
1707                      gw_tmp_nd(ifreq,iflavor,iflavor1)= gw_tmp_nd(ifreq,iflavor,iflavor1) + &
1708 &                     u_nl*gl_nd(ileg,iflavor,iflavor1)
1709                    end do
1710                  end do
1711 
1712            !    write(77,*) "------------------", gw_tmp_nd(ifreq,1,1)
1713 
1714                end do
1715            !  write(77,*) "------------------ sum ", gw_tmp_nd(ifreq,1,1)
1716              end do
1717              ABI_DEALLOCATE( jbes )
1718              call xmpi_barrier(paw_dmft%spacecomm)
1719              call flush_unit(std_out)
1720            end if
1721            gw_tmp_nd = gw_tmp_nd*Ha_eV
1722 
1723 
1724 
1725 ! =========================================================================
1726 
1727            if ( leg_measure ) then !only if functionality is enabled
1728              ABI_DEALLOCATE(gl_nd)
1729            end if
1730 
1731          else
1732            write(message,'(2a)') ch10," Can't launch TRIQS CTHYB solver because dmftqmc_l must be >= 2*dmft_nwli + 1"
1733            MSG_ERROR(message)
1734          end if
1735 
1736        end if
1737 
1738 ! =================================================================
1739 !    END CTQMC run TRIQS
1740 ! =================================================================
1741 
1742      else if (testcode>=1) then
1743 
1744 
1745 ! =================================================================
1746 !    CTQMC run for tests
1747 ! =================================================================
1748        write(std_out,*) "nomega,dmftqmc_l",nomega,paw_dmft%dmftqmc_l
1749        call CtqmcInterface_run(hybrid,fw1(1:nomega,:),Gtau=gtmp,&
1750 &       Gw=gw_tmp,E=green%ecorr_qmc(iatom),&
1751 &       matU=umod,opt_levels=levels_ctqmc)
1752 
1753 ! for non diagonal code
1754 !       call CtqmcInterface_run(hybrid,fw1_nd(1:nomega,:,:),Gtau=gtmp_nd,&
1755 !&       Gw=gw_tmp_nd,D=Doccsum,E=green%ecorr_qmc(iatom),&
1756 !&       Noise=Noise,matU=umod,opt_levels=levels_ctqmc,hybri_limit=hybri_limit)
1757 
1758      end if
1759 
1760 ! =================================================================
1761 !  TEST
1762 !  If test of the code is activated, and testrot =1 rotate back green's function
1763 !  and stop the code.
1764 ! =================================================================
1765      if(testcode==1) then
1766 
1767        call testcode_ctqmc(paw_dmft%dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,&
1768 &       levels_ctqmc,hybri_limit,nflavor,2,paw_dmft%temp,testrot,testcode,umod)
1769 
1770        write(message,'(2a)') ch10,' testcode end of test calculation'
1771        MSG_ERROR(message)
1772      end if
1773      if(testcode==2) then
1774        write(message,'(2a)') ch10,' testcode 2 end of test calculation'
1775        MSG_ERROR(message)
1776      end if
1777 
1778 ! =================================================================
1779      ABI_DEALLOCATE(fw1)
1780      ABI_DEALLOCATE(fw1_nd)
1781 
1782      !----------------------------------------
1783      ! <DEBUG>
1784      !----------------------------------------
1785      ! Construct UNIT
1786      if(paw_dmft%idmftloop < 10) then
1787        write(gtau_iter,'("0",i1)') paw_dmft%idmftloop
1788      elseif(paw_dmft%idmftloop >= 10 .and. paw_dmft%idmftloop < 100) then
1789        write(gtau_iter,'(i2)') paw_dmft%idmftloop
1790      else
1791        gtau_iter="xx"
1792      end if
1793      if(iatom < 10) then
1794        write(iatomnb,'("0",i1)') iatom
1795      elseif(iatom >= 10 .and. iatom < 100) then
1796        write(iatomnb,'(i2)') iatom
1797      else
1798        iatomnb='xx'
1799      end if
1800 
1801      if(paw_dmft%myproc .eq. mod(nproc+1,nproc)) then
1802 ! < HACK >
1803        if(paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7) then
1804          if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /=0) then
1805            MSG_ERROR(message)
1806          end if
1807          do ifreq=1,paw_dmft%dmft_nwli
1808            write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq),((gw_tmp_nd(ifreq,iflavor,iflavor)), iflavor=1, nflavor)
1809          end do
1810          close(unt)
1811        else
1812          if(paw_dmft%dmft_solv==5) then
1813            if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_"//gtau_iter//".dat", message, newunit=unt) /= 0) then
1814              MSG_ERROR(message)
1815            end if
1816            do itau=1,paw_dmft%dmftqmc_l
1817              write(unt,'(29f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,&
1818              (gtmp(itau,iflavor), iflavor=1, nflavor)
1819            end do
1820            write(unt,'(29f21.14)') 1/paw_dmft%temp, (-1_dp-gtmp(1,iflavor), iflavor=1, nflavor)
1821            close(unt)
1822          endif
1823          if(paw_dmft%dmft_solv==8) then
1824            if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_offdiag_unsym_"//gtau_iter//".dat",&
1825 &           message, newunit=unt) /= 0) then
1826              MSG_ERROR(message)
1827            end if
1828            do itau=1,paw_dmft%dmftqmc_l
1829              write(unt,'(196f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,&
1830              ((gtmp_nd(itau,iflavor,iflavor1), iflavor=1, nflavor),iflavor1=1, nflavor)
1831            end do
1832            close(unt)
1833            ABI_DATATYPE_ALLOCATE(matlu1,(natom))
1834            call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
1835            do itau=1,paw_dmft%dmftqmc_l
1836              do isppol=1,nsppol
1837                do ispinor1=1,nspinor
1838                  do im1=1,tndim
1839                    iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
1840                    do ispinor2=1,nspinor
1841                      do im2=1,tndim
1842                        iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
1843                        matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
1844 &                        gtmp_nd(itau,iflavor1,iflavor2)
1845                      end do  ! im2
1846                    end do  ! ispinor2
1847                  end do  ! im1
1848                end do  ! ispinor
1849              end do ! isppol
1850              call rotate_matlu(matlu1,eigvectmatlu,natom,3,0)
1851              call slm2ylm_matlu(matlu1,natom,2,0)
1852              call sym_matlu(cryst_struc,matlu1,pawang)
1853              call slm2ylm_matlu(matlu1,natom,1,0)
1854              call rotate_matlu(matlu1,eigvectmatlu,natom,3,1)
1855              do isppol=1,nsppol
1856                do ispinor1=1,nspinor
1857                  do im1=1,tndim
1858                    iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
1859                    do ispinor2=1,nspinor
1860                      do im2=1,tndim
1861                        iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
1862                        gtmp_nd(itau,iflavor1,iflavor2)=&
1863                         matlu1(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)
1864                      end do  ! im2
1865                    end do  ! ispinor2
1866                  end do  ! im1
1867                end do  ! ispinor
1868              end do ! isppol
1869            end do  !itau
1870            call destroy_matlu(matlu1,natom)
1871            ABI_DATATYPE_DEALLOCATE(matlu1)
1872            if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gtau_offdiag_"//gtau_iter//".dat",&
1873 &           message, newunit=unt) /= 0) then
1874              MSG_ERROR(message)
1875            end if
1876            do itau=1,paw_dmft%dmftqmc_l
1877              write(unt,'(196f21.14)') float(itau-1)/float(paw_dmft%dmftqmc_l)/paw_dmft%temp,&
1878              ((gtmp_nd(itau,iflavor,iflavor1), iflavor=1, nflavor),iflavor1=1, nflavor)
1879            end do
1880            close(unt)
1881          endif
1882          !open(unit=4243, file=trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_F_"//gtau_iter//".dat")
1883          !call BathOperator_printF(paw_dmft%hybrid(iatom)%hybrid%bath,4243) !Already comment here
1884          !close(4243)
1885          if(paw_dmft%dmft_solv==5) then
1886            if (open_file(trim(paw_dmft%filapp)//"_atom_"//iatomnb//"_Gw_"//gtau_iter//".dat", message, newunit=unt) /= 0) then
1887              MSG_ERROR(message)
1888            end if
1889            do ifreq=1,paw_dmft%dmft_nwlo
1890              write(unt,'(29f21.14)') paw_dmft%omega_lo(ifreq), &
1891 &             (gw_tmp(ifreq,iflavor), iflavor=1, nflavor)
1892            end do
1893          endif
1894          close(unt)
1895        end if
1896 ! </ HACK >
1897      end if
1898 
1899      !----------------------------------------
1900      ! </DEBUG>
1901      !----------------------------------------
1902      if(paw_dmft%dmft_solv>=6.and.paw_dmft%dmft_solv<=7) then
1903      !Nothing just hybrid var problem
1904      else
1905        write(message,'(a,2x,a)') ch10,&
1906 &       " == Destroy CTQMC"
1907        call wrtout(std_out,message,'COLL')
1908        if(paw_dmft%dmft_solv==5) call CtqmcInterface_finalize(hybrid)
1909        if(paw_dmft%dmft_solv==8) call CtqmcoffdiagInterface_finalize(hybridoffdiag)
1910        write(message,'(a,2x,a)') ch10,&
1911 &       " == Destroy CTQMC done"
1912        call wrtout(std_out,message,'COLL')
1913      end if
1914      !ABI_DEALLOCATE(f_with_k)
1915      ABI_DEALLOCATE(hybri_limit)
1916      ABI_DEALLOCATE(levels_ctqmc_nd)
1917      ABI_DEALLOCATE(levels_ctqmc)
1918 
1919 ! ___________________________________________________________________________________
1920 !
1921 !  FOURTH PART : USE OUTPUT OF CTQMC AND DO BACK ROTATION
1922 ! ___________________________________________________________________________________
1923 !
1924      do itau=1,paw_dmft%dmftqmc_l
1925        green%oper_tau(itau)%matlu(iatom)%mat(:,:,:,:,:)=czero
1926      end do
1927      green%occup_tau%matlu(iatom)%mat(nflavor:,:,:,:,:)=czero
1928 
1929      do ifreq=1,paw_dmft%dmft_nwlo
1930        green%oper(ifreq)%matlu(iatom)%mat(:,:,:,:,:)=czero
1931      end do
1932      green%occup%matlu(iatom)%mat(:,:,:,:,:)=czero
1933 
1934 !   built time and frequency green's function from output of CTQMC
1935 ! =================================================================
1936      if(opt_nondiag==1) then
1937        do isppol=1,nsppol
1938          do ispinor1=1,nspinor
1939            do im1=1,tndim
1940              iflavor1=im1+tndim*(ispinor1-1)+tndim*(isppol-1)
1941              do ispinor2=1,nspinor
1942                do im2=1,tndim
1943                  iflavor2=im2+tndim*(ispinor2-1)+tndim*(isppol-1)
1944                   ! NNtodo here: nsppol=1 above, only, with symetrisation
1945                   ! automatic
1946 !               iflavor1=im+(isppol-1)*tndim
1947                  do itau=1,paw_dmft%dmftqmc_l
1948                    green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
1949 &                   gtmp_nd(itau,iflavor1,iflavor2)
1950                    if(nsppol==1.and.nspinor==1) then
1951                      green%oper_tau(itau)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
1952 &                     (gtmp_nd(itau,iflavor1,iflavor2)+gtmp_nd(itau,iflavor1+tndim,iflavor2+tndim))/two
1953                    end if
1954                   ! NNtodo here: isppol above should be one and symetrized
1955                   ! gtmp
1956                  end do  !itau
1957 !               ifreq2=0
1958                  if(paw_dmft%dmft_solv<6.or.leg_measure) then
1959                    do ifreq=1,paw_dmft%dmft_nwlo
1960 !                   if(paw_dmft%select_log(ifreq)==1) then
1961 !                     ifreq2=ifreq2+1
1962                      green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
1963 &                     gw_tmp_nd(ifreq,iflavor1,iflavor2)
1964                      if(nsppol==1.and.nspinor==1) then
1965                        green%oper(ifreq)%matlu(iatom)%mat(im1,im2,isppol,ispinor1,ispinor2)=&
1966 &                       (gw_tmp_nd(ifreq,iflavor1,iflavor2)+&
1967 &                       gw_tmp_nd(ifreq,iflavor1+tndim,iflavor2+tndim))/two
1968                      end if
1969                     ! NNtodo here: isppol above should be one and symetrized
1970                     ! gw_tmp
1971 !                   endif
1972                    end do ! ifreq
1973                  end if
1974                end do  ! im2
1975              end do  ! ispinor2
1976            end do  ! im1
1977          end do  ! ispinor
1978        end do ! isppol
1979        if(paw_dmft%dmft_solv>=6.and..not.leg_measure) then
1980          write(message,'(2a,i3,13x,a)') ch10,'   ===  Direct Fourier Transform t->w of Weiss Field'
1981          call wrtout(std_out,message,'COLL')
1982          call fourier_green(cryst_struc,green,paw_dmft,&
1983 &         pawang,opt_ksloc=2,opt_tw=1)
1984        end if
1985      else
1986        iflavor=0
1987        do isppol=1,nsppol
1988          do ispinor=1,nspinor
1989            do im=1,tndim
1990               ! NNtodo here: nsppol=1 above, only, with symetrisation
1991               ! automatic
1992 !           iflavor=im+(isppol-1)*tndim
1993              iflavor=iflavor+1
1994              do itau=1,paw_dmft%dmftqmc_l
1995                green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gtmp(itau,iflavor)
1996                if(nsppol==1.and.nspinor==1) then
1997                  green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=&
1998 &                 (gtmp(itau,iflavor)+gtmp(itau,iflavor+tndim))/two
1999                end if
2000               ! NNtodo here: isppol above should be one and symetrized
2001               ! gtmp
2002              end do
2003 !           ifreq2=0
2004              do ifreq=1,paw_dmft%dmft_nwlo
2005 !             if(paw_dmft%select_log(ifreq)==1) then
2006 !               ifreq2=ifreq2+1
2007                green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=gw_tmp(ifreq,iflavor)
2008                if(nsppol==1.and.nspinor==1) then
2009                  green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=&
2010 &                 (gw_tmp(ifreq,iflavor)+gw_tmp(ifreq,iflavor+tndim))/two
2011                end if
2012                ! if opt_nondiag=0, then:
2013                ! As Green's function is diagonal, one suppress off diag  terms in Weiss, if any.
2014                ! (If off diag are non zero in the density matrix and thus in the Green's function,
2015                !  there is a warning in checkreal_matlu above).
2016                do im1=1,tndim
2017                  do ispinor1=1,nspinor
2018                    if(im/=im1.or.ispinor/=ispinor1) then
2019                      weiss%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=czero
2020                    end if
2021                  end do
2022                end do
2023               ! NNtodo here: isppol above should be one and symetrized
2024               ! gw_tmp
2025 !             endif
2026              end do
2027            end do
2028          end do
2029        end do
2030      end if
2031      if(paw_dmft%dmft_solv<6) then
2032       ! write(6,*) "dmft_solv",paw_dmft%dmft_solv
2033        ABI_DEALLOCATE(gw_tmp)
2034      endif
2035      ABI_DEALLOCATE(gw_tmp_nd)
2036      ABI_DEALLOCATE(gtmp)
2037      ABI_DEALLOCATE(gtmp_nd)
2038      if(nsppol==1.and.nspinor==1) then
2039        write(message,'(a,2x,a,f13.5)') ch10,&
2040 &       " == nsppol==1 and nspden==1: Green functions from CTQMC have been symetrized over spin"
2041        call wrtout(std_out,message,'COLL')
2042      end if
2043      !write(message,'(i3,4x,2e21.14)') 5,weiss_for_rot%oper(1)%matlu(1)%mat(1,1,1,1,1)
2044      !call wrtout(std_out,message,'COLL')  ! debug
2045 !     do im=1,tndim
2046 !       do itau=1,paw_dmft%dmftqmc_l
2047 !         gtt=(green%oper_tau(itau)%matlu(iatom)%mat(im,im,1,1,1)+&
2048 !&         green%oper_tau(itau)%matlu(iatom)%mat(im,im,2,1,1))/two
2049 !         green%oper_tau(itau)%matlu(iatom)%mat(im,im,1,1,1)=gtt
2050 !         green%oper_tau(itau)%matlu(iatom)%mat(im,im,2,1,1)=gtt
2051 !       enddo
2052 !     enddo
2053 !     write(6,*)" SYMETRISATION OVER ISPPOL"
2054      !deallocate(correl_loc,f_with_k,gtmp,fw1,fw2)
2055    end if
2056  end do ! iatom
2057 
2058 
2059  if(paw_dmft%dmft_prgn==1) then
2060    call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=2)
2061    call print_green('QMC_diag_notsym',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
2062  end if
2063  !write(message,'(i3,4x,2e21.14)') 6,weiss_for_rot%oper(1)%matlu(1)%mat(1,1,1,1,1)
2064  !call wrtout(std_out,message,'COLL')  ! debug
2065 ! =================================================================
2066 ! Inverse Weiss, then
2067 ! Copy Weiss_for_rot into weiss and rotate back weiss to the original basis
2068 ! =================================================================
2069 
2070 ! ABI_ALLOCATE(shift,(natom))
2071 ! do ifreq=1,paw_dmft%dmft_nwlo
2072 !  ! First weiss_for_rot contains -G_0^-1+iw_n
2073 !  ! -------------------------------------------
2074 !  ! Compute G_0^-1-iw_n
2075 !  ! --------------------
2076 !       write(6,*) "1"
2077 !  if(opt_fk==1) call fac_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,-cone)
2078 !
2079 !
2080 !       write(6,*) "2"
2081 !  ! Compute G_0^-1
2082 !  ! --------------------
2083 !  shift(:)=cmplx(zero,paw_dmft%omega_lo(ifreq),kind=dp)
2084 !  if(opt_fk==1) call shift_matlu(weiss_for_rot%oper(ifreq)%matlu,natom,shift,signe=1)
2085 !
2086 !       write(6,*) "3"
2087 !  ! Compute G_0
2088 !  ! --------------------
2089 !   call inverse_oper(weiss_for_rot%oper(ifreq),option=1,prtopt=1)
2090 !   ! No need to copy if weiss_for_rot is a pointer to weiss ...
2091 !!   if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0)
2092 !!   if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
2093 !
2094 !  ! Compute G_0 in the original basis
2095 !  ! --------------------
2096 !   call rotate_matlu(weiss_for_rot%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
2097 ! end do
2098 ! ABI_DEALLOCATE(shift)
2099 
2100 ! =================================================================
2101 ! Here compute Self energy from Dyson and print it
2102 ! Warning : Weiss_for_rot is inversed inside dyson
2103 ! =================================================================
2104 ! call initialize_self(self,paw_dmft)
2105 ! call dyson(green,paw_dmft,self,weiss_for_rot,opt_weissself=2)
2106 ! call rw_self(self,mpi_enreg,paw_dmft,prtopt=2,opt_rw=2,opt_char="diag")
2107 ! call destroy_self(self)
2108  !write(message,'(i3,4x,2e21.14)') 7,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1)
2109  !call wrtout(std_out,message,'COLL')  ! debug
2110 
2111 ! =================================================================
2112 ! Rotate back green function to original basis (non-diagonal)
2113 !  (and Weiss for further use: might be useful if an back Fourier
2114 !     transformation is done).
2115 ! =================================================================
2116  if(pawprtvol>=3) then
2117    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
2118 &  " == Print green function for small tau after CTQMC"  ! debug
2119    call wrtout(std_out,message,'COLL')  ! debug
2120    call print_matlu(green%oper_tau(1)%matlu,natom,1)  ! debug
2121    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
2122 &  " == Print green function for small freq after CTQMC"  ! debug
2123    call wrtout(std_out,message,'COLL')  ! debug
2124    call print_matlu(green%oper(1)%matlu,natom,1)  ! debug
2125  end if
2126 
2127 
2128 !  === Compute rotated Occupations in green%occup_tau
2129  call occup_green_tau(green)
2130 
2131  if(pawprtvol>=3) then
2132 !  === Compute non rotated Occupations in green%occup_tau
2133    write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the ctqmc basis"
2134    call wrtout(std_out,message,'COLL')
2135    call print_matlu(green%occup_tau%matlu,natom,1)
2136  end if
2137 
2138  write(message,'(a,2x,a,f13.5)') ch10,&
2139 & " == Rotate Green function to original basis "
2140  call wrtout(std_out,message,'COLL')
2141  !write(message,'(i3,4x,2e21.14)') 8,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1)
2142  !call wrtout(std_out,message,'COLL')  ! debug
2143 
2144  ! Rotate oper_tau into Ylm basis and then Slm basis
2145  !-------------------------------------------------------------
2146  do itau=1,paw_dmft%dmftqmc_l
2147    if(opt_diag/=0) call rotate_matlu(green%oper_tau(itau)%matlu,eigvectmatlu,natom,3,0)
2148    if(useylm==1) call slm2ylm_matlu(green%oper_tau(itau)%matlu,natom,2,0)
2149  end do
2150 
2151  ! Rotate occup_tau into Ylm basis and then Slm basis
2152 
2153  ! Rotate occup_tau into Ylm basis and then Slm basis
2154  !-------------------------------------------------------------
2155  if(opt_diag/=0) call rotate_matlu(green%occup_tau%matlu,eigvectmatlu,natom,3,0)
2156  if(useylm==1) then
2157    write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Ylm basis"
2158    call wrtout(std_out,message,'COLL')
2159    call print_matlu(green%occup_tau%matlu,natom,1)
2160  end if
2161 
2162  if(useylm==1) call slm2ylm_matlu(green%occup_tau%matlu,natom,2,0)
2163  write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the Slm basis"
2164  call wrtout(std_out,message,'COLL')
2165  call print_matlu(green%occup_tau%matlu,natom,1)
2166 
2167  ! Rotate Green's and Weiss functions into Ylm basis and then Slm basis
2168  !-------------------------------------------------------------
2169  do ifreq=1,paw_dmft%dmft_nwlo
2170    if(opt_diag/=0) call rotate_matlu(green%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
2171    if(useylm==1) call slm2ylm_matlu(green%oper(ifreq)%matlu,natom,2,0)
2172    if(opt_diag/=0) call rotate_matlu(weiss%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
2173    if(useylm==1) call slm2ylm_matlu(weiss%oper(ifreq)%matlu,natom,2,0)
2174  end do
2175  !write(message,'(i3,4x,2e21.14)') 10,weiss%oper(1)%matlu(1)%mat(1,1,1,1,1)
2176  !call wrtout(std_out,message,'COLL')  ! debug
2177 
2178  if(pawprtvol>=3) then
2179    write(message,'(a,2x,a,f13.5)') ch10,&                  ! debug
2180 &  " == Print green function for small time after rotation (in the original basis)" ! debug
2181    call wrtout(std_out,message,'COLL')                  ! debug
2182    call print_matlu(green%oper_tau(1)%matlu,natom,1)  ! debug
2183    write(message,'(a,2x,a,f13.5)') ch10,&                  ! debug
2184 &  " == Print green function for small freq after rotation (in the original basis)" ! debug
2185    call wrtout(std_out,message,'COLL')                  ! debug
2186    call print_matlu(green%oper(1)%matlu,natom,1)  ! debug
2187    !< HACK >
2188    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
2189 &  " == Print diagonalized weiss_for_rot function after rotation for small freq in the ctqmc basis"  ! debug
2190    call wrtout(std_out,message,'COLL')  ! debug
2191    call print_matlu(weiss_for_rot%oper(1)%matlu,natom,1)  ! debug
2192    !</ HACK >
2193    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
2194 &  " == Print weiss function for small freq in the original basis"  ! debug
2195    call wrtout(std_out,message,'COLL')  ! debug
2196    call print_matlu(weiss%oper(1)%matlu,natom,1)  ! debug
2197 
2198    do ifreq=1,paw_dmft%dmft_nwlo
2199      call sym_matlu(cryst_struc,weiss%oper(ifreq)%matlu,pawang)
2200    end do
2201    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
2202 &  " == Print symetrized weiss function for small freq in the original basis"  ! debug
2203    call wrtout(std_out,message,'COLL')  ! debug
2204    call print_matlu(weiss%oper(1)%matlu,natom,1)  ! debug
2205  end if
2206 
2207 
2208  ABI_DATATYPE_ALLOCATE(matlu1,(natom))
2209  call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,matlu1)
2210  call copy_matlu(green%occup_tau%matlu,matlu1,natom)
2211  call sym_matlu(cryst_struc,matlu1,pawang)
2212 
2213  write(message,'(a,2x,a,f13.5)') ch10," == Occupation from G(tau) in the original basis"
2214  call wrtout(std_out,message,'COLL')
2215  call print_matlu(green%occup_tau%matlu,natom,1)
2216 
2217  write(message,'(a,2x,a,f13.5)') ch10," == Symetrized occupations"
2218  call wrtout(std_out,message,'COLL')
2219  call print_matlu(matlu1,natom,1)
2220 
2221  call diff_matlu("CTQMC Occup","CTQMC Occup symetrized",green%occup_tau%matlu,matlu1,natom,0,tol4,ierr)
2222  call destroy_matlu(matlu1,natom)
2223  ABI_DATATYPE_DEALLOCATE(matlu1)
2224 
2225 ! =================================================================
2226 ! Symetrise green function G(tau) and G(ifreq) to recover symetry
2227 ! artificially broken by QMC
2228 ! =================================================================
2229  write(message,'(a,2x,a,f13.5)') ch10,&
2230 & " == Symetrise green function after QMC "
2231  call wrtout(std_out,message,'COLL')
2232  do itau=1,paw_dmft%dmftqmc_l
2233    call sym_matlu(cryst_struc,green%oper_tau(itau)%matlu,pawang)
2234  end do
2235  do ifreq=1,paw_dmft%dmft_nwlo
2236    call sym_matlu(cryst_struc,green%oper(ifreq)%matlu,pawang)
2237  end do
2238  if(pawprtvol>=3) then
2239    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
2240 &  " == Print green function for small time after symetrisation"  !  debug
2241    call wrtout(std_out,message,'COLL')  ! debug
2242    call print_matlu(green%oper_tau(1)%matlu,natom,1)  ! debug
2243    write(message,'(a,2x,a,f13.5)') ch10,&  ! debug
2244 &  " == Print green function for small freq after symetrisation"  !  debug
2245    call wrtout(std_out,message,'COLL')  ! debug
2246    call print_matlu(green%oper(1)%matlu,natom,1)  ! debug
2247  end if
2248  if(paw_dmft%dmft_prgn==1) then
2249    call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=2)
2250    call print_green('QMC_sym',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
2251  end if
2252 
2253 !  === Compute Occupations  (Symetrized from oper_tau)
2254  call occup_green_tau(green)
2255 
2256 
2257 !  === Print occupations
2258  call printocc_green(green,6,paw_dmft,3)
2259 
2260  call destroy_oper(energy_level)
2261  call destroy_matlu(dmat_diag,natom)
2262  ABI_DATATYPE_DEALLOCATE(dmat_diag)
2263  call destroy_matlu(identity,natom)
2264  ABI_DATATYPE_DEALLOCATE(identity)
2265  do iatom=1,cryst_struc%natom
2266    lpawu=paw_dmft%lpawu(iatom)
2267    if(lpawu/=-1) then
2268      do isppol=1,nsppol
2269        ABI_DEALLOCATE(eigvectmatlu(iatom,isppol)%value)
2270        !ABI_DEALLOCATE(udens_atoms(iatom))
2271      end do
2272      ABI_DEALLOCATE(udens_atoms(iatom)%value)
2273    end if
2274  end do
2275  ABI_DATATYPE_DEALLOCATE(udens_atoms)
2276  ABI_DATATYPE_DEALLOCATE(eigvectmatlu)
2277  call destroy_green(weiss_for_rot)
2278 ! call destroy_green(gw_loc)
2279 ! call destroy_green(greenlda)
2280 
2281 !  destroy limit of hybridization
2282  call destroy_matlu(hybri_coeff,paw_dmft%natom)
2283  ABI_DATATYPE_DEALLOCATE(hybri_coeff)
2284 
2285 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.

COPYRIGHT

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

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

PARENTS

      qmc_prep_ctqmc

CHILDREN

SOURCE

2332 subroutine testcode_ctqmc(dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,levels_ctqmc,hybri_limit,&
2333 &   nflavor,opt,temp,testrot,testcode,umod)
2334 
2335  use defs_basis
2336  use defs_datatypes
2337  use defs_abitypes
2338  use m_errors
2339  use m_ctqmc
2340  use m_CtqmcInterface
2341  use m_greenhyb
2342  !use m_self, only : self_type,initialize_self,destroy_self,print_self,rw_self
2343  use m_io_tools, only : flush_unit
2344 
2345 !This section has been created automatically by the script Abilint (TD).
2346 !Do not modify the following lines by hand.
2347 #undef ABI_FUNC
2348 #define ABI_FUNC 'testcode_ctqmc'
2349 !End of the abilint section
2350 
2351  implicit none
2352 
2353 !Arguments ------------------------------------
2354 !scalars
2355  integer, intent(in) :: dmftqmc_l,nflavor,testrot,testcode,opt
2356  real(dp), intent(in) :: temp
2357  real(dp), intent(out) :: umod(2,2)
2358  complex(dpc), intent(inout) :: gw_tmp_nd(:,:,:)
2359  real(dp),  intent(inout) :: gtmp_nd(:,:,:)
2360  complex(dpc), intent(out) :: fw1(:,:)
2361  complex(dpc), intent(out) :: fw1_nd(:,:,:)
2362  real(dp),  intent(inout) :: levels_ctqmc(:)
2363  complex(dpc),  intent(inout) :: hybri_limit(:,:)
2364 
2365 !Local variables ------------------------------
2366  character(len=500) :: message
2367  integer :: ifreq, itau,realrot,simplehyb
2368  real(dp) :: omega
2369  real(dp) :: tbi1,tbi2,e2,tbi3,tbi4,e3,e4,tbi21,tbi12,e3b,e4b,tbi21b,tbi12b
2370  complex(dpc) :: e1
2371 ! arrays
2372  complex(dpc) :: RR(2,2)
2373  complex(dpc) :: RR1(2,2)
2374  complex(dpc) :: RRi(2,2)
2375  complex(dpc) :: RRt(2,2)
2376 ! ************************************************************************
2377  if (testcode==0) return
2378  if (nflavor/=2) then
2379    write(message,'(2a)') ch10,' testcode nflavor.ne.2'
2380    MSG_ERROR(message)
2381  end if
2382 
2383  simplehyb=2
2384  simplehyb=1
2385  simplehyb=3
2386  !=========================
2387  ! Built rotation matrix
2388  !=========================
2389  realrot=0
2390  realrot=2
2391  if (realrot==1) then
2392    ! Real rotation
2393    !=========================
2394    RR(1,1)  =  SQRT(3.d0)/2.d0
2395    RR(1,2)  = -1.d0/2.d0
2396    RR(2,1)  =  1.d0/2.d0
2397    RR(2,2)  =  SQRT(3.d0)/2.d0
2398  else if (realrot==2) then
2399    ! Real rotation
2400    !=========================
2401    RR(1,1)  =  SQRT(1.d0/2.d0)
2402    RR(1,2)  = -SQRT(1.d0/2.d0)
2403    RR(2,1)  =  SQRT(1.d0/2.d0)
2404    RR(2,2)  =  SQRT(1.d0/2.d0)
2405  else
2406    ! Complex rotation
2407    !=========================
2408    RR(1,1)  =  CMPLX(one,two)
2409    RR(1,2)  =  CMPLX(one,one)
2410    RR(2,1)  =  CMPLX(one,-one)
2411    RR(2,2)  =  CMPLX(-one,two)
2412    RR=RR/sqrt(seven)
2413  end if
2414  ! Check rotation is unitary
2415  !==========================
2416  RRi(1,1) =  conjg(RR(1,1))
2417  RRi(1,2) =  conjg(RR(2,1))
2418  RRi(2,1) =  conjg(RR(1,2))
2419  RRi(2,2) =  conjg(RR(2,2))
2420  RR1(:,:)  = MATMUL ( RR(:,:) , RRi(:,:)          )
2421  !write(6,*) "RR1",RR1
2422  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
2423    write(message,'(2a)') ch10,' testcode error in rotation matrix'
2424    MSG_ERROR(message)
2425  end if
2426 
2427 
2428  !=================================
2429  ! Built hybridization  for CTQMC
2430  !=================================
2431  if (opt==1) then
2432 
2433  !  Parameters: tight-binding + U
2434  !  firt test of the code try umod=0, and (tbi1,tbi2,e1,e2)=(2,1,0.5,0.0) testrot=1
2435  !  second test of the code try umod=four, and (tbi1,tbi2,e1,e2)=(2,1,0.0,0.0) testrot=1
2436  !=======================================================================================
2437    fw1_nd(:,:,:)= czero
2438    tbi1=2.0_dp
2439    tbi2=1.0_dp
2440    tbi3=1.0_dp
2441    tbi4=1.0_dp
2442    tbi12=2.5_dp
2443    tbi12b=2.5_dp
2444    tbi21=2.5_dp
2445    tbi21b=2.5_dp
2446    e1=cmplx(0.0,0.0,8)
2447    e2=zero
2448    e3=0.2
2449    e4=0.3
2450    e3b=0.3
2451    e4b=-0.2
2452    umod(:,:)=0.d0
2453 
2454    if(testrot==1.and.(abs(tbi1-tbi2)<tol6)) then
2455      write(message,'(3a)') ch10,' testrot=1 with tbi1=tbi2 is equivalent' &
2456      ,'to testrot=0: change testrot'
2457      MSG_WARNING(message)
2458    end if
2459    ! Built fw1_nd
2460    !==============
2461    do ifreq=1,dmftqmc_l
2462 
2463      omega=pi*temp*(two*float(ifreq)-1)
2464 
2465      if(simplehyb==1) then
2466        fw1_nd(ifreq,1,1) =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
2467        fw1_nd(ifreq,2,2) =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
2468        fw1(ifreq,1)      =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
2469        fw1(ifreq,2)      =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
2470        hybri_limit(1,1)=tbi1**2
2471        hybri_limit(2,2)=tbi2**2
2472        hybri_limit(1,2)=0.d0
2473        hybri_limit(2,1)=0.d0
2474      else if(simplehyb==2) then
2475        fw1_nd(ifreq,1,1) =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)+tbi3**2/(dcmplx(0.d0,omega)-e3)
2476        fw1_nd(ifreq,2,2) =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)+tbi4**2/(dcmplx(0.d0,omega)-e4)
2477        fw1(ifreq,1)      =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
2478        fw1(ifreq,2)      =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
2479      else if(simplehyb==3) then
2480        fw1_nd(ifreq,1,1) =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
2481        fw1_nd(ifreq,2,2) =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
2482        fw1_nd(ifreq,1,2) =  tbi12**2/(dcmplx(0.d0,omega)-e3)+tbi12b**2/(dcmplx(0.d0,omega)-e3b)
2483        fw1_nd(ifreq,2,1) =  tbi21**2/(dcmplx(0.d0,omega)-e4)+tbi21b**2/(dcmplx(0.d0,omega)-e4b)
2484        fw1(ifreq,1)      =  -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)
2485        fw1(ifreq,2)      =  -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)
2486        hybri_limit(1,1)=tbi1**2
2487        hybri_limit(2,2)=tbi2**2
2488        hybri_limit(1,2)=tbi12**2+tbi12b**2
2489        hybri_limit(2,1)=tbi21**2+tbi21b**2
2490      end if
2491      write(132,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1))
2492      write(133,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2))
2493      write(134,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1))
2494      write(135,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2))
2495      write(1234,*) omega, real(fw1(ifreq,1)),aimag(fw1(ifreq,1))
2496    end do
2497    ! Built level and limit of hybridization
2498    !=======================================
2499    levels_ctqmc(1:nflavor)=-umod(1,1)/two
2500 
2501    write(std_out,*) "fw1_nd"
2502    write(std_out,*) fw1_nd(1,1,1), fw1_nd(1,1,2)
2503    write(std_out,*) fw1_nd(1,2,1), fw1_nd(1,2,2)
2504    write(std_out,*) "fw1"
2505    write(std_out,*) fw1(1,1), fw1(1,2)
2506    write(std_out,*) fw1(2,1), fw1(2,2)
2507 
2508  ! Rotate hybridization if testrot=1
2509  !==================================
2510    if(testrot==1) then
2511 
2512      do ifreq=1,dmftqmc_l
2513        RRt(:,:)  = MATMUL ( RR(:,:)  , fw1_nd(ifreq,:,:) )
2514    !write(6,*) "RRt"
2515    !write(6,*) RRt(1,1), RRt(1,2)
2516    !write(6,*) RRt(2,1), RRt(2,2)
2517        RR1(:,:)  = MATMUL ( RRt(:,:) , RRi(:,:)          )
2518    !write(6,*) "RR1"
2519    !write(6,*) RR1(1,1), RR1(1,2)
2520    !write(6,*) RR1(2,1), RR1(2,2)
2521        fw1_nd(ifreq,:,:)=RR1(:,:)
2522        omega=pi*temp*(two*float(ifreq)+1)
2523        write(3322,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1))
2524        write(232,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1))
2525        write(233,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2))
2526        write(234,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1))
2527        write(235,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2))
2528      end do
2529 
2530      ! Rotate limit of hybridization
2531      !=======================================
2532      RRt(:,:)  = MATMUL ( RR(:,:)  , hybri_limit(:,:)  )
2533      RR1(:,:)  = MATMUL ( RRt(:,:) , RRi(:,:)          )
2534      hybri_limit(:,:)=RR1(:,:)
2535 
2536    end if
2537    ! rajouter test real(fw1_nd(1,:,:)) doit etre diagonale
2538 
2539  !======================================
2540  ! Rotate Green's function from CTQMC
2541  !======================================
2542  else if(opt==2) then
2543 
2544    write(std_out,*) "gw_tmp_nd"
2545    write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2)
2546    write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2)
2547    ! Rotate Green's function back
2548    !==============================
2549    if(testrot==1) then
2550      do ifreq=1,dmftqmc_l
2551        RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gw_tmp_nd(ifreq,1:nflavor,1:nflavor) )
2552        RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) )
2553        gw_tmp_nd(ifreq,1:nflavor,1:nflavor)=RR1(1:nflavor,1:nflavor)
2554      end do
2555 
2556      write(std_out,*) "gw_tmp_nd after rotation"
2557      write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2)
2558      write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2)
2559 
2560      do itau=1,dmftqmc_l
2561        RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) )
2562        RR1(1:nflavor,1:nflavor)  = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) )
2563        gtmp_nd(itau,1:nflavor,1:nflavor)=real(RR1(1:nflavor,1:nflavor))
2564      end do
2565 
2566    ! Rotate Green's function for comparison with testrot=1
2567    !======================================================
2568    else if (testrot==0) then ! produce rotated green's function to compare to testrot=1 case
2569 
2570      do itau=1,dmftqmc_l
2571        RRt(1:nflavor,1:nflavor) = MATMUL ( RR(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) )
2572        RR1(1:nflavor,1:nflavor)  = MATMUL ( RRt(1:nflavor,1:nflavor),RRi(1:nflavor,1:nflavor) )
2573        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))
2574      end do
2575 
2576    end if
2577 
2578    ! Print out rotated Green's function
2579    !=====================================
2580    do itau=1,dmftqmc_l
2581      write(555,'(e14.5,4(2e14.5,3x))') real(itau-1)/(temp*real(dmftqmc_l)),gtmp_nd(itau,1,1),&
2582 &     gtmp_nd(itau,2,2),gtmp_nd(itau,1,2),gtmp_nd(itau,2,1)
2583    end do
2584 
2585    write(message,'(2a)') ch10,' testcode end of test calculation'
2586    MSG_ERROR(message)
2587 
2588  end if
2589  close(444)
2590  close(555)
2591 
2592 end subroutine testcode_ctqmc