TABLE OF CONTENTS


ABINIT/greenldacompute_green [ Functions ]

[ Top ] [ Functions ]

NAME

 greenldacompute_green

FUNCTION

 Compute levels for ctqmc

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

OUTPUT

PARENTS

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

3199  subroutine greenldacompute_green(cryst_struc,green,pawang,paw_dmft)
3200 
3201  use defs_basis
3202  use defs_datatypes
3203  use defs_abitypes
3204  use m_errors
3205  use m_abicore
3206 
3207  use m_pawang, only : pawang_type
3208  use m_crystal, only : crystal_t
3209  use m_paw_dmft, only : paw_dmft_type
3210  use m_oper, only : oper_type,loc_oper
3211  use m_matlu, only : sym_matlu, print_matlu
3212 
3213 !This section has been created automatically by the script Abilint (TD).
3214 !Do not modify the following lines by hand.
3215 #undef ABI_FUNC
3216 #define ABI_FUNC 'greenldacompute_green'
3217 !End of the abilint section
3218 
3219  implicit none
3220 
3221 !Arguments ------------------------------------
3222 !scalars
3223  type(crystal_t),intent(in) :: cryst_struc
3224  type(pawang_type), intent(in) :: pawang
3225  type(paw_dmft_type), intent(in)  :: paw_dmft
3226  type(green_type),intent(inout) :: green
3227 
3228 !Local variables ------------------------------
3229 ! scalars
3230  integer :: ifreq,iband,ikpt,isppol
3231  integer :: mbandc,natom,nspinor,nsppol,nkpt
3232  character(len=500) :: message
3233 ! arrays
3234 !************************************************************************
3235 
3236  mbandc=paw_dmft%mbandc
3237  nkpt=paw_dmft%nkpt
3238  nsppol=paw_dmft%nsppol
3239  nspinor=paw_dmft%nspinor
3240  natom=paw_dmft%natom
3241 
3242      !  write(6,*) green%oper(1)%ks(1,1,1,1)
3243  if(green%oper(1)%has_operks==0) then
3244   MSG_ERROR("greenlda%oper(1)%ks not allocated")
3245  endif
3246 
3247 !======================================
3248 !Get Green's function G=1/(iw_n+mu-e_nk)
3249 !======================================
3250  do ifreq=1,green%nw
3251    do iband=1,mbandc
3252      do ikpt=1,nkpt
3253        do isppol=1,nsppol
3254      !  write(6,*) green%oper(ifreq)%ks(isppol,ikpt,iband,iband)
3255      !  write(6,*) ifreq,iband,ikpt,isppol
3256      !  write(6,*) cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)
3257      !  write(6,*) paw_dmft%fermie
3258      !  write(6,*) paw_dmft%eigen_lda(isppol,ikpt,iband)
3259          green%oper(ifreq)%ks(isppol,ikpt,iband,iband)=&
3260          cone/(cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)+paw_dmft%fermie-paw_dmft%eigen_lda(isppol,ikpt,iband))
3261        end do
3262      end do
3263    end do
3264 
3265 
3266 !======================================================================
3267 ! Compute local Green's function
3268 !======================================================================
3269    call loc_oper(green%oper(ifreq),paw_dmft,1)
3270 
3271 !======================================================================
3272 ! Symetrize
3273 !======================================================================
3274    call sym_matlu(cryst_struc,green%oper(ifreq)%matlu,pawang)
3275 
3276  enddo
3277  write(message,'(a,2x,a,f13.5)') ch10," == Print LDA Green's function for last frequency"
3278  call wrtout(std_out,message,'COLL')
3279  call print_matlu(green%oper(paw_dmft%dmft_nwlo)%matlu,natom,1)
3280 
3281  end subroutine greenldacompute_green

ABINIT/m_green [ Modules ]

[ Top ] [ Modules ]

NAME

  m_green

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 
25 #include "abi_common.h"
26 
27  MODULE m_green
28 
29  use defs_basis
30  use defs_abitypes
31  use m_xmpi
32  use m_abicore
33  use m_errors
34  use m_lib_four
35 
36  use m_io_tools, only : flush_unit, open_file
37  use m_time,     only : timab
38  use m_oper, only : oper_type
39  use m_matlu, only : matlu_type
40 
41  implicit none
42 
43  private
44 
45  public :: init_green
46  public :: destroy_green
47  public :: init_green_tau
48  public :: destroy_green_tau
49  public :: print_green
50  public :: printocc_green
51  public :: compute_green
52  public :: integrate_green
53  public :: icip_green
54  public :: fourier_green
55  public :: check_fourier_green
56  public :: compa_occup_ks
57  public :: copy_green
58  public :: occup_green_tau
59  public :: add_int_fct
60  public :: int_fct
61  public :: fourier_fct
62  public :: spline_fct
63  private :: occupfd
64  private :: distrib_paral
65  public :: greenldacompute_green
66  public :: fermi_green
67  public :: newton
68  public :: local_ks_green

m_green/add_int_fct [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 add_int_fct

FUNCTION

 Do integration in matsubara space

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

  ff= function is frequency space
  ldiag    = option according to diagonal or non-diagonal elements
  option = nspinor
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  proct= for parallelism

SIDE EFFECTS

  * integral = integral of ff over matsubara frequencies (there is an accumulation in the present routine, so intent(inout))
  * ft= function is time space

PARENTS

      m_green

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2494 subroutine add_int_fct(ifreq,ff,ldiag,omega_current,option,integral,temp,wgt_wlo,dmft_nwlo)
2495 
2496  use defs_basis
2497  use m_paw_dmft, only : paw_dmft_type
2498 
2499 !This section has been created automatically by the script Abilint (TD).
2500 !Do not modify the following lines by hand.
2501 #undef ABI_FUNC
2502 #define ABI_FUNC 'add_int_fct'
2503 !End of the abilint section
2504 
2505  implicit none
2506 
2507 !Arguments ------------------------------------
2508 !type
2509  integer,intent(in) :: ifreq
2510  logical,intent(in) :: ldiag
2511  integer,intent(in) :: option,dmft_nwlo
2512  complex(dpc),intent(inout) :: integral
2513  complex(dpc), intent(in) :: ff
2514  complex(dpc), intent(in) :: omega_current
2515  real(dp), intent(in) :: temp, wgt_wlo
2516 
2517 !local variables-------------------------------
2518  real(dp)  :: omega
2519 ! *********************************************************************
2520  omega=aimag(omega_current)
2521  if(ldiag) then
2522 
2523   if(option==1) then ! nspinor==1
2524 !   write(500,*) paw_dmft%omega_lo(ifreq),real(ff(ifreq)),imag(ff(ifreq))
2525     integral=integral+2.d0*temp *                         &
2526 &      real( ff-one / ( j_dpc*omega ) ) *  &
2527 &      wgt_wlo
2528     if(ifreq==dmft_nwlo) integral=integral+half
2529 !    integral=integral+half
2530     ! the if is here, to count only one time this correction
2531   endif
2532 
2533   if(option==2) then ! nspinor==2
2534     integral=integral+2.d0*temp *                         &
2535 &       ( ff-one / ( j_dpc*omega ) ) *  &
2536 &   wgt_wlo
2537     if(ifreq==dmft_nwlo) integral=integral+half
2538 !    integral=integral+half
2539     ! the if is here, to count only one time this correction
2540   endif
2541 
2542 
2543  else   ! ldiag
2544 
2545 ! write(std_out,*) "nondiag"
2546   if(option==1) then
2547     integral=integral+2.d0*temp *   &
2548 &    real( ff ) *                     &
2549 &   wgt_wlo
2550   endif
2551   if(option==2) then
2552     integral=integral+2.d0*temp *   &
2553 &        ff   *                     &
2554 &    wgt_wlo
2555   endif
2556  endif  ! ldiag
2557 
2558 end subroutine add_int_fct

m_green/check_fourier_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 check_fourier_green

FUNCTION

  Check fourier transformations

INPUTS

  cryst_struc <type(crystal_t)>= crystal structure data.
  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawang <type(pawang)>=paw angular mesh and related data

OUTPUT

PARENTS

      m_dmft

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2312 subroutine check_fourier_green(cryst_struc,green,paw_dmft,pawang)
2313 
2314  use m_pawang, only : pawang_type
2315  use m_crystal, only : crystal_t
2316  use m_paw_dmft, only : paw_dmft_type
2317  use m_matlu, only : print_matlu
2318 
2319 !This section has been created automatically by the script Abilint (TD).
2320 !Do not modify the following lines by hand.
2321 #undef ABI_FUNC
2322 #define ABI_FUNC 'check_fourier_green'
2323 !End of the abilint section
2324 
2325  implicit none
2326 
2327 !Arguments ------------------------------------
2328 !type
2329  type(crystal_t),intent(in) :: cryst_struc
2330  type(green_type),intent(inout) :: green
2331  !type(MPI_type), intent(in) :: mpi_enreg
2332  type(pawang_type),intent(in) :: pawang
2333  type(paw_dmft_type), intent(inout) :: paw_dmft
2334 
2335 !local variables-------------------------------
2336  type(green_type) :: green_check
2337  character(len=500) :: message
2338 ! *********************************************************************
2339 ! Only imaginary frequencies here
2340  if(green%w_type=="real") then
2341    message = 'check_fourier_green not implemented for real frequency'
2342    MSG_BUG(message)
2343  endif
2344 
2345  call init_green(green_check,paw_dmft)
2346  call init_green_tau(green_check,paw_dmft)
2347  call copy_green(green,green_check,opt_tw=2)
2348 
2349  write(message,'(2a,i3,13x,a)') ch10,'   ===  Inverse Fourier Transform w->t of Weiss Field'
2350  call wrtout(std_out,message,'COLL')
2351  call fourier_green(cryst_struc,green_check,paw_dmft,&
2352 & pawang,opt_ksloc=2,opt_tw=-1)
2353 
2354  write(message,'(3a)') ch10,' === Print (for check by user) of occupation matrix'&
2355 &   ,' after  fourier transform with respect to initial one'
2356  call wrtout(std_out,message,'COLL')
2357  call printocc_green(green_check,6,paw_dmft,3)
2358 
2359  write(message,'(2a,i3,13x,a)') ch10,'   ===  Direct Fourier Transform t->w of Weiss Field'
2360  call wrtout(std_out,message,'COLL')
2361  call fourier_green(cryst_struc,green_check,paw_dmft,&
2362 & pawang,opt_ksloc=2,opt_tw=1)
2363 ! call print_matlu(green%oper(1)%matlu,paw_dmft%natom,1) ! debug
2364 
2365  call integrate_green(cryst_struc,green_check,paw_dmft,&
2366 & pawang,prtopt=2,opt_ksloc=2)
2367 
2368  write(message,'(3a)') ch10,' === Print (for check by user) of occupation matrix'&
2369 &   ,' after double fourier transform with respect to initial one'
2370  call wrtout(std_out,message,'COLL')
2371  call printocc_green(green_check,5,paw_dmft,3)
2372 
2373  call destroy_green_tau(green_check)
2374  call destroy_green(green_check)
2375 
2376 end subroutine check_fourier_green

m_green/compa_occup_ks [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 compa_occup_ks

FUNCTION

  Compare occupations to Fermi Dirac Occupations

INPUTS

  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data

OUTPUT

PARENTS

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2399 subroutine compa_occup_ks(green,paw_dmft)
2400 
2401  use m_paw_dmft, only : paw_dmft_type
2402 
2403 !This section has been created automatically by the script Abilint (TD).
2404 !Do not modify the following lines by hand.
2405 #undef ABI_FUNC
2406 #define ABI_FUNC 'compa_occup_ks'
2407 !End of the abilint section
2408 
2409  implicit none
2410 
2411 !Arguments ------------------------------------
2412 !type
2413  type(green_type),intent(inout) :: green
2414  type(paw_dmft_type), intent(inout) :: paw_dmft
2415 
2416 !local variables-------------------------------
2417  character(len=500) :: message
2418  integer :: ib,ikpt,isppol
2419  integer :: ib_,ikpt_,isppol_
2420  integer :: ib__,ikpt__,isppol__
2421  real(dp) :: diffmax,occ1,occ2,occ_a,occ_b,diffrel,occ_aa,occ_bb
2422 ! *********************************************************************
2423  diffmax=zero
2424  diffrel=zero
2425  do isppol=1,paw_dmft%nsppol
2426    do ikpt=1,paw_dmft%nkpt
2427      do ib=1,paw_dmft%mbandc
2428        occ1=occupfd(paw_dmft%eigen_lda(isppol,ikpt,ib),paw_dmft%fermie,paw_dmft%temp)
2429        occ2=real(green%occup%ks(isppol,ikpt,ib,ib))
2430        if(abs(occ1-occ2)>diffmax) then
2431          diffmax=abs(occ1-occ2)
2432          occ_a=occ1
2433          occ_b=occ2
2434          ib_=ib;isppol_=isppol;ikpt_=ikpt
2435        endif
2436        if(abs(two*(occ1-occ2)/(occ1+occ2))>diffrel) then
2437          diffrel=abs(two*(occ1-occ2)/(occ1+occ2))
2438          occ_aa=occ1
2439          occ_bb=occ2
2440          ib__=ib;isppol__=isppol;ikpt__=ikpt
2441        endif
2442      enddo
2443    enddo
2444  enddo
2445 
2446 
2447  write(message,'(2a)') ch10,'   ===  Compare green function occupations and Fermi Dirac occupations'
2448  call wrtout(std_out,message,'COLL')
2449  write(message,'(2a,f12.5,2a,f12.5,2a,f12.5,2a,3i5)') ch10,'     =  Max difference is',diffmax,&
2450 &                           ch10,'        Corresponding Occupation from green function is',occ_b,&
2451 &                           ch10,'        Corresponding Occupation Fermi Dirac weight  is',occ_a,&
2452 &                           ch10,'        (For polarization, k-point and band index)     ',isppol_,ikpt_,ib_
2453  call wrtout(std_out,message,'COLL')
2454  write(message,'(2a,f12.5,2a,f12.5,2a,f12.5,2a,3i5)') ch10,'     =  Max relative difference is',diffrel,&
2455 &                           ch10,'        Corresponding Occupation from green function is',occ_bb,&
2456 &                           ch10,'        Corresponding Occupation Fermi Dirac weight  is',occ_aa,&
2457 &                           ch10,'        (For polarization, k-point and band index)     ',isppol__,ikpt__,ib__
2458  call wrtout(std_out,message,'COLL')
2459 
2460 end subroutine compa_occup_ks

m_green/compute_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 compute_green

FUNCTION

 compute green function from LDA and self-energy

INPUTS

  pawang <type(pawang)>=paw angular mesh and related data
  opt_self = optional argument, if =1, upfold self-energy
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  green  <type(green_type)>= green function data
  opt_nonxsum = 0 : do usual xsum after calculation of green(freq)%ks
              = 1 : do not do xsum after calculation of green(freq)%ks: each proc as only a part
                    of the data: this is useful where only the total number of electron will be computed.
  opt_nonxsum2 0 : green(ifreq)%matlu will be broadcasted
               1 : green(ifreq)%matlu will not be broadcasted in compute_green: calc
                   if occupations will not possible.
                   (a keyword: compute_local_green would be in fact equivalent and more clear)

OUTPUT

PARENTS

      m_dmft,fermi_green,m_green,newton,spectral_function

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

1036 subroutine compute_green(cryst_struc,green,paw_dmft,pawang,prtopt,self,opt_self,&
1037 &           opt_nonxsum,opt_nonxsum2)
1038 
1039  use m_pawang, only : pawang_type
1040  use m_crystal, only : crystal_t
1041  use m_matlu, only : sym_matlu,zero_matlu,add_matlu,print_matlu,shift_matlu,init_matlu,destroy_matlu,copy_matlu
1042  use m_oper, only : inverse_oper,loc_oper,print_oper,upfold_oper,init_oper,destroy_oper
1043  use m_paw_dmft, only : paw_dmft_type
1044  use m_self, only : self_type
1045 
1046 !This section has been created automatically by the script Abilint (TD).
1047 !Do not modify the following lines by hand.
1048 #undef ABI_FUNC
1049 #define ABI_FUNC 'compute_green'
1050 !End of the abilint section
1051 
1052  implicit none
1053 
1054 !Arguments ------------------------------------
1055 !type
1056  type(crystal_t),intent(in) :: cryst_struc
1057  type(green_type),intent(inout) :: green
1058  type(paw_dmft_type), intent(in) :: paw_dmft
1059  !type(MPI_type), intent(in) :: mpi_enreg
1060  type(pawang_type),intent(in) :: pawang
1061  type(self_type), intent(inout) :: self
1062  integer, intent(in) :: prtopt
1063  integer, optional, intent(in) :: opt_self
1064  integer, optional, intent(in) :: opt_nonxsum
1065  integer, optional, intent(in) :: opt_nonxsum2
1066 
1067 !local variables-------------------------------
1068  !logical :: lintegrate
1069  integer :: iatom,ib,ib1,ierr,ifreq,ikpt,is
1070  integer :: icomp_chloc
1071  integer :: natom,mband,mbandc,nkpt
1072  integer :: myproc,nproc,spacecomm,nspinor,nsppol,option,optself,optnonxsum
1073  integer :: optnonxsum2
1074  integer, allocatable :: procb_ifreq(:)
1075  character(len=500) :: message
1076  real(dp) :: fermilevel
1077  real(dp), allocatable :: Id(:,:)
1078 ! integer, allocatable :: procb(:,:),proct(:,:)
1079  real(dp) :: tsec(2)
1080  type(oper_type) :: self_minus_hdc_oper
1081  complex(dpc) :: omega_current
1082  type(oper_type) :: green_temp
1083 ! *********************************************************************
1084 
1085  !lintegrate=.true.
1086  !if(lintegrate.and.green%w_type=="real") then
1087  !if(green%w_type=="real") then
1088  !  message = 'integrate_green not implemented for real frequency'
1089  !  MSG_BUG(message)
1090  !endif
1091  call timab(624,1,tsec)
1092  if(present(opt_self)) then
1093    optself=opt_self
1094  else
1095    optself=0
1096  endif
1097  if(present(opt_nonxsum)) then
1098    optnonxsum=opt_nonxsum
1099  else
1100    optnonxsum=0
1101  endif
1102  if(present(opt_nonxsum2)) then
1103    optnonxsum2=opt_nonxsum2
1104  else
1105    optnonxsum2=0
1106  endif
1107 
1108  if(prtopt>0)  then
1109    write(message,'(2a,i3,13x,a)') ch10,'  ===  Compute green function '
1110    call wrtout(std_out,message,'COLL')
1111  endif
1112 
1113  if(self%nw/=green%nw)  then
1114    message = ' BUG: frequencies for green and self not coherent'
1115    MSG_BUG(message)
1116  endif
1117 
1118 ! Initialise spaceComm, myproc, and nproc
1119  spacecomm=paw_dmft%spacecomm
1120  myproc=paw_dmft%myproc
1121  nproc=paw_dmft%nproc
1122 
1123 ! Initialise integers
1124  mband   = paw_dmft%mband
1125  mbandc  = paw_dmft%mbandc
1126  nkpt    = paw_dmft%nkpt
1127  nspinor = paw_dmft%nspinor
1128  nsppol  = paw_dmft%nsppol
1129  natom   = paw_dmft%natom
1130 
1131 ! Initialise for compiler
1132  omega_current=czero
1133  icomp_chloc=0
1134 
1135 ! ==============================
1136 ! Initialise Identity
1137 ! ==============================
1138  ABI_ALLOCATE(Id,(paw_dmft%mbandc,paw_dmft%mbandc))
1139  Id=zero
1140  do ib=1, paw_dmft%mbandc
1141    Id(ib,ib)=one
1142  enddo ! ib
1143 
1144  if(prtopt/=0.and.prtopt>-100)  then
1145    write(message,'(2a)') ch10,'  == Green function is computed:'
1146    call wrtout(std_out,message,'COLL')
1147  endif
1148  option=1
1149  fermilevel=paw_dmft%fermie
1150  if(option==123) then
1151    fermilevel=2.d0
1152    write(message,'(2a,e14.3,a)') ch10,&
1153 &  '  Warning (special case for check: fermi level=',fermilevel,')'
1154    call wrtout(std_out,message,'COLL')
1155  endif
1156 
1157 ! ====================================================
1158 ! Upfold self-energy and double counting  Self_imp -> self(k)
1159 ! ====================================================
1160 ! if(optself==1) then
1161 !   do ifreq=1,green%nw
1162 !     call upfold_oper(self%oper(ifreq),paw_dmft,1)
1163 !   enddo ! ifreq
1164 !   call upfold_oper(self%hdc,paw_dmft,1)
1165 ! endif
1166 
1167 ! =================================================================
1168 ! Initialize green%oper before calculation (important for xmpi_sum)
1169 ! =================================================================
1170  do ifreq=1,green%nw
1171    call zero_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom)
1172  enddo ! ifreq
1173  !call xmpi_barrier(spacecomm)
1174 
1175 ! ================================
1176 ! == Compute Green function G(k)
1177 ! ================================
1178  green%occup%ks=czero
1179  ABI_ALLOCATE(procb_ifreq,(paw_dmft%nkpt))
1180  do ifreq=1,green%nw
1181    !if(present(iii)) write(6,*) ch10,'ifreq  self', ifreq,self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
1182 !   ====================================================
1183 !   First Upfold self-energy and double counting  Self_imp -> self(k)
1184 !   ====================================================
1185 !   if(mod(ifreq-1,nproc)==myproc) then
1186 !    write(6,*) "compute_green ifreq",ifreq, mod(ifreq-1,nproc)==myproc,proct(ifreq,myproc)==1
1187    if(green%proct(ifreq,myproc)==1) then
1188      if(green%w_type=="imag") then
1189        omega_current=cmplx(zero,green%omega(ifreq),kind=dp)
1190      else if(green%w_type=="real") then
1191        omega_current=cmplx(green%omega(ifreq),0.1/27.211_dp,kind=dp)
1192      endif
1193      call init_oper(paw_dmft,self_minus_hdc_oper)
1194      call init_oper(paw_dmft,green_temp)
1195 
1196      call add_matlu(self%oper(ifreq)%matlu,self%hdc%matlu,&
1197 &                   self_minus_hdc_oper%matlu,green%oper(ifreq)%natom,-1)
1198      if(paw_dmft%dmft_solv==4)  then
1199        call shift_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,cmplx(self%qmc_shift,0.d0,kind=dp),-1)
1200        call shift_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,cmplx(self%qmc_xmu,0.d0,kind=dp),-1)
1201      endif
1202      if(optself==0) then
1203        call zero_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom)
1204      end if
1205 
1206      procb_ifreq=green%procb(ifreq,:)
1207      call upfold_oper(self_minus_hdc_oper,paw_dmft,1,procb=procb_ifreq,iproc=myproc)
1208      do ib1 = 1 , paw_dmft%mbandc
1209        do ib = 1 , paw_dmft%mbandc
1210          do ikpt = 1 , paw_dmft%nkpt
1211            do is = 1 , paw_dmft%nsppol
1212              if (green%procb(ifreq,ikpt)==myproc) then
1213 !               green%oper(ifreq)%ks(is,ikpt,ib,ib1)=       &
1214                green_temp%ks(is,ikpt,ib,ib1)=       &
1215 &               ( omega_current     &
1216 &               + fermilevel                               &
1217 &               - paw_dmft%eigen_lda(is,ikpt,ib)) * Id(ib,ib1) &
1218 &               - self_minus_hdc_oper%ks(is,ikpt,ib,ib1)
1219 
1220 
1221 
1222 
1223 !&               - (self%oper(ifreq)%ks(is,ikpt,ib,ib1)-self%hdc%ks(is,ikpt,ib,ib1))
1224 !               if(prtopt>5) then
1225 !               if(ikpt==1.and.(ifreq==1.or.ifreq==3).and.ib==16.and.ib1==16) then
1226 !                write(std_out,*) 'omega_current                         ',omega_current
1227 !                write(std_out,*) 'fermilevel                            ',fermilevel
1228 !                write(std_out,*) ' paw_dmft%eigen_lda(is,ikpt,ib)       ', paw_dmft%eigen_lda(is,ikpt,ib),Id(ib,ib1)
1229 !                write(std_out,*) 'self_minus_hdc_oper%ks(is,ikpt,ib,ib1)',self_minus_hdc_oper%ks(is,ikpt,ib,ib1)
1230 !                write(std_out,*) 'green                                 ',green%oper(ifreq)%ks(is,ikpt,ib,ib1)
1231 !               endif
1232 !               if(ib==1.and.ib1==3) then
1233 !                 write(std_out,*) "ff compute",ikpt,ifreq,is,ikpt,ib,ib1
1234 !                 write(std_out,*) "ff compute",ikpt,ifreq, green_temp%ks(is,ikpt,ib,ib1)
1235 !                 write(std_out,*) "ff details",paw_dmft%eigen_lda(is,ikpt,ib)
1236 !                 write(std_out,*) "ff details2",fermilevel
1237 !                 write(std_out,*) "ff details3",Id(ib,ib1)
1238 !        !        write(std_out,*) "ff details4",self_minus_hdc_oper%ks(is,ikpt,ib,ib1)
1239 !               endif
1240              endif
1241            enddo ! is
1242          enddo ! ikpt
1243        enddo ! ib
1244      enddo ! ib1
1245 !     call print_oper(green%oper(ifreq),9,paw_dmft,3)
1246 !       write(std_out,*) 'after print_oper'
1247 !     if(ifreq==1.or.ifreq==3) then
1248 !         write(std_out,*) 'after print_oper', ifreq
1249 !       write(std_out,*) 'green1  ifreq         %ks(1,1,16,16)',ifreq,green%oper(ifreq)%ks(1,1,16,16)
1250 !     endif
1251 !       write(std_out,*) 'before inverse_oper'
1252     call inverse_oper(green_temp,2,prtopt,procb=procb_ifreq,iproc=myproc)
1253      if (green%oper(1)%has_operks==1) green%oper(ifreq)%ks=green_temp%ks
1254     !if(ifreq==1) then
1255     !  write(std_out,*) "1188",green_temp%ks(1,1,8,8)
1256     !  write(std_out,*) "1189",green_temp%ks(1,1,8,9)
1257     !  write(std_out,*) "1198",green_temp%ks(1,1,9,8)
1258     !  write(std_out,*) "1199",green_temp%ks(1,1,9,9)
1259     !endif
1260 
1261      !if(lintegrate) then
1262 !    accumulate integration
1263      if(green%w_type/="real") then
1264        do is = 1 , paw_dmft%nsppol
1265          do ib = 1 , paw_dmft%mbandc
1266            do ib1 = 1 , paw_dmft%mbandc
1267              do ikpt = 1 , paw_dmft%nkpt
1268                if (green%procb(ifreq,ikpt)==myproc) then
1269                  call add_int_fct(ifreq,green_temp%ks(is,ikpt,ib,ib1),ib==ib1,    &
1270 &                      omega_current,2,green%occup%ks(is,ikpt,ib,ib1),            &
1271 &                      paw_dmft%temp,paw_dmft%wgt_wlo(ifreq),paw_dmft%dmft_nwlo)
1272                endif
1273              enddo ! ikpt
1274            enddo ! ib1
1275          enddo ! ib
1276        enddo ! is
1277      endif
1278 !        write(std_out,*) 'after inverse_oper'
1279 !      if(ikpt==1.and.is==1.and.ib==1.and.ib1==1) then
1280 !        write(6,*) 'occup(is,ikpt,ib,ib1)',ifreq,green%occup%ks(1,1,1,1),green_temp%ks(1,1,1,1)
1281 !      endif
1282 !        write(std_out,*) 'green1afterinversion  %ks(1,1,16,16)',ifreq,green%oper(ifreq)%ks(1,1,16,16)
1283 !      endif
1284 !        write(std_out,*) 'before flush'
1285 !      call flush(std_out)
1286      !endif
1287 ! ================================
1288 ! == Compute Local Green function
1289 ! ================================
1290    !write(message,'(2a)') ch10,' loc'
1291    !call wrtout(std_out,message,'COLL')
1292    !call flush(std_out)
1293      call loc_oper(green_temp,paw_dmft,1,procb=procb_ifreq,iproc=myproc)
1294     !if(ifreq==1) then
1295     !  write(std_out,*) "4411",green_temp%matlu(1)%mat(4,4,1,1,1)
1296     !  write(std_out,*) "4512",green_temp%matlu(1)%mat(4,5,1,1,2)
1297     !  write(std_out,*) "5421",green_temp%matlu(1)%mat(5,4,1,2,1)
1298     !  write(std_out,*) "5522",green_temp%matlu(1)%mat(5,5,1,2,2)
1299     !  write(std_out,*) "(5512)",green_temp%matlu(1)%mat(5,5,1,1,2)
1300     !endif
1301 !     write(std_out,*) ifreq,nproc,'before if after loc_oper'
1302 !     if(ifreq==1.or.ifreq==11) then
1303 !       write(std_out,*) ifreq,nproc,'before sym'
1304 !       call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0)
1305 !       ! ok
1306 !     endif
1307 ! call flush(std_out)
1308      call sym_matlu(cryst_struc,green_temp%matlu,pawang)
1309      call copy_matlu(green_temp%matlu,green%oper(ifreq)%matlu,natom)
1310 !     if(ifreq==1.and.ifreq==11) then
1311 !       write(std_out,*) ifreq,nproc,'after sym'
1312 !       call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0)
1313 !       ! ok
1314 !     endif
1315      call destroy_oper(self_minus_hdc_oper)
1316      call destroy_oper(green_temp)
1317 ! call flush(std_out)
1318    endif ! parallelisation
1319  enddo ! ifreq
1320  ABI_DEALLOCATE(procb_ifreq)
1321 
1322 ! =============================================
1323 ! built total green function (sum over procs).
1324 ! =============================================
1325  !call xmpi_barrier(spacecomm)
1326  do ifreq=1,green%nw
1327 !  call xmpi_sum(green%oper(ifreq)%ks,spacecomm,ierr)
1328 !       print *, "myproc, proct, ifreq ------------------- ", myproc, ifreq
1329 !   do ikpt=1,paw_dmft%nkpt
1330 !     call xmpi_bcast(green%oper(ifreq)%ks(:,ikpt,:,:),procb(ifreq,ikpt),spacecomm,ierr)
1331 !   enddo
1332 !! or
1333    if(optnonxsum==0.and.green%oper(ifreq)%has_operks==1) then
1334      call xmpi_sum(green%oper(ifreq)%ks,spacecomm,ierr)
1335    else if(optnonxsum==0.and.green%oper(ifreq)%has_operks==0) then
1336      message = 'optnonxsum==0 and green%oper(ifreq)%has_operks==0: not compatible'
1337      MSG_BUG(message)
1338    endif
1339 !   endif
1340 ! enddo ! ifreq
1341 !  print *,"myproc", myproc
1342    if(optnonxsum2==0) then
1343       do iatom=1,green%oper(ifreq)%natom
1344        if(green%oper(ifreq)%matlu(iatom)%lpawu.ne.-1) then
1345          call xmpi_sum(green%oper(ifreq)%matlu(iatom)%mat,spacecomm,ierr)
1346        endif
1347      enddo ! iatom
1348      green%has_greenmatlu_xsum=1
1349    else if(optnonxsum2==1) then
1350      green%has_greenmatlu_xsum=0
1351    endif
1352 !     if(ifreq==1.or.ifreq==11) then
1353 !       write(std_out,*) ifreq,nproc,'after xsum'
1354 !       call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0)
1355 !       ! ok
1356 !     endif
1357  enddo ! ifreq
1358  !call xmpi_barrier(spacecomm)
1359  call xmpi_sum(green%occup%ks,spacecomm,ierr)
1360 ! write(std_out,*) 'afterxsum sym     %matlu(1)%mat(2,5,1,1,1) 1',green%oper(1)%matlu(1)%mat(2,5,1,1,1)
1361 
1362  if(prtopt/=0.and.prtopt>-100)  then
1363    write(message,'(2a)') ch10,&
1364 &   '  == Local Green function has been computed and projected on local orbitals'
1365    call wrtout(std_out,message,'COLL')
1366  endif
1367 ! useless test
1368  if(abs(prtopt)>=4.and.prtopt>-100) then
1369    write(message,'(2a)') ch10,' == Green function is now printed for first frequency'
1370    call wrtout(std_out,message,'COLL')
1371    call print_oper(green%oper(1),9,paw_dmft,3)
1372    write(message,'(2a)') ch10,' == Green function is now printed for second frequency'
1373    call wrtout(std_out,message,'COLL')
1374    call print_oper(green%oper(2),9,paw_dmft,3)
1375    if(paw_dmft%dmft_nwlo>=11) then
1376      write(message,'(2a)') ch10,' == Green function is now printed for 11th frequency'
1377      call wrtout(std_out,message,'COLL')
1378      call print_oper(green%oper(11),9,paw_dmft,3)
1379    endif
1380  endif
1381 ! call flush(std_out)
1382 
1383  ABI_DEALLOCATE(Id)
1384  call timab(624,2,tsec)
1385 
1386 end subroutine compute_green

m_green/copy_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 copy_green

FUNCTION

  copy one data structure green1 into green2

INPUTS

  green1  <type(green_type)>= green function data
  green2  <type(green_type)>= green function data
  opt_tw = option to precise which data to copy
          1: copy only green%occup_tau and green%oper_tau data
          2: copy only green%occup and green%oper  data (frequency)

OUTPUT

PARENTS

      dyson,impurity_solve,m_green,qmc_prep_ctqmc,spectral_function

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

549 subroutine copy_green(green1,green2,opt_tw)
550 
551  use m_crystal, only : crystal_t
552  use m_oper, only : copy_oper
553 
554 !This section has been created automatically by the script Abilint (TD).
555 !Do not modify the following lines by hand.
556 #undef ABI_FUNC
557 #define ABI_FUNC 'copy_green'
558 !End of the abilint section
559 
560  implicit none
561 
562 !Arguments ------------------------------------
563 !type
564  type(green_type),intent(in) :: green1
565  type(green_type),intent(inout) :: green2
566  integer, intent(in) :: opt_tw
567 
568 !local variables-------------------------------
569  integer :: ifreq,itau
570 ! *********************************************************************
571 
572  if(opt_tw==2) then
573    call copy_oper(green1%occup,green2%occup)
574    do ifreq=1, green1%nw
575      call copy_oper(green1%oper(ifreq),green2%oper(ifreq))
576      if(green1%has_greenmatlu_xsum==1) then ! indicate to integrate_green  that xsum has been done
577 !  for matlu in compute_green.
578         green2%has_greenmatlu_xsum=1
579      endif
580    enddo
581  else if (opt_tw==1) then
582    call copy_oper(green1%occup_tau,green2%occup_tau)
583    do itau=1,green1%dmftqmc_l
584      call copy_oper(green1%oper_tau(itau),green2%oper_tau(itau))
585    enddo
586  endif
587 
588 end subroutine copy_green

m_green/destroy_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 destroy_green

FUNCTION

  deallocate green

INPUTS

  green  <type(green_type)>= green function data

OUTPUT

PARENTS

      m_dmft,dyson,m_hubbard_one,m_green,qmc_prep_ctqmc,spectral_function

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

408 subroutine destroy_green(green)
409 
410  use m_oper, only : destroy_oper
411 
412 !This section has been created automatically by the script Abilint (TD).
413 !Do not modify the following lines by hand.
414 #undef ABI_FUNC
415 #define ABI_FUNC 'destroy_green'
416 !End of the abilint section
417 
418  implicit none
419 
420 !Arguments ------------------------------------
421 !scalars
422  type(green_type),intent(inout) :: green
423 
424 !local variables-------------------------------
425  integer :: ifreq
426 
427 ! *********************************************************************
428  call destroy_oper(green%occup)
429  if ( allocated(green%oper))       then
430    do ifreq=1,green%nw
431     call destroy_oper(green%oper(ifreq))
432    enddo
433    ABI_DATATYPE_DEALLOCATE(green%oper)
434  end if
435  if ( allocated(green%charge_matlu))       then
436    ABI_DEALLOCATE(green%charge_matlu)
437  end if
438  green%has_charge_matlu=0
439 
440  if ( allocated(green%charge_matlu_prev))       then
441    ABI_DEALLOCATE(green%charge_matlu_prev)
442  end if
443  green%has_charge_matlu_prev=0
444 
445  if ( allocated(green%charge_matlu_solver))       then
446    ABI_DEALLOCATE(green%charge_matlu_solver)
447  end if
448  green%has_charge_matlu_solver=0
449 
450  if ( allocated(green%ecorr_qmc))   then
451     ABI_DEALLOCATE(green%ecorr_qmc)
452  end if
453  if ( allocated(green%procb))   then
454     ABI_DEALLOCATE(green%procb)
455  end if
456  if ( allocated(green%proct))   then
457     ABI_DEALLOCATE(green%proct)
458  end if
459  green%omega => null()
460 
461 end subroutine destroy_green

m_green/destroy_green_tau [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 destroy_green_tau

FUNCTION

  deallocate green

INPUTS

  green  <type(green_type)>= green function data

OUTPUT

PARENTS

      impurity_solve,m_green

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

484 subroutine destroy_green_tau(green)
485 
486  use m_crystal, only : crystal_t
487  use m_oper, only : destroy_oper
488 
489 !This section has been created automatically by the script Abilint (TD).
490 !Do not modify the following lines by hand.
491 #undef ABI_FUNC
492 #define ABI_FUNC 'destroy_green_tau'
493 !End of the abilint section
494 
495  implicit none
496 
497 !Arguments ------------------------------------
498 !scalars
499  type(green_type),intent(inout) :: green
500 ! integer, optional, intent(in) :: opt_ksloc
501 !local variables-------------------------------
502  integer :: itau
503 ! integer :: optksloc
504 ! *********************************************************************
505 ! if(present(opt_ksloc)) then
506 !   optksloc=opt_ksloc
507 ! else
508 !   optksloc=3
509 ! endif
510 
511  call destroy_oper(green%occup_tau)
512  if ( allocated(green%oper_tau))       then
513    do itau=1,green%dmftqmc_l
514     call destroy_oper(green%oper_tau(itau))
515    enddo
516    ABI_DATATYPE_DEALLOCATE(green%oper_tau)
517  end if
518  if ( allocated(green%tau))            then
519    ABI_DEALLOCATE(green%tau)
520  end if
521 
522 end subroutine destroy_green_tau

m_green/distrib_paral [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 distrib_paral

FUNCTION

INPUTS

  nw : number of frequencies
  nkpt : number of k-points
  nproc : number of procs

OUTPUT

  procb(iw,ikpt) :  number of the proc  that is in charge of the combination {iw,ikpt}.
  proct(iw,iproc) : 1 if the frequency "iw" should be computed by the proc "iproc"
                    0 if iw should not   "      "
                    careful: if proct=1, it does not mean that each combination {iw,ikpt} is treated by this proc.

PARENTS

      m_green

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

3071  subroutine distrib_paral(nkpt,nproc,nw,nw_perproc,procb,proct)
3072 
3073  use defs_basis
3074  use defs_datatypes
3075  use defs_abitypes
3076 
3077 !This section has been created automatically by the script Abilint (TD).
3078 !Do not modify the following lines by hand.
3079 #undef ABI_FUNC
3080 #define ABI_FUNC 'distrib_paral'
3081 !End of the abilint section
3082 
3083  implicit none
3084 
3085 !Arguments ------------------------------------
3086 !type
3087 ! Integrate analytic tail 1/(iw-mu)
3088  integer, intent(in) :: nw,nkpt,nproc
3089  integer, intent(out) :: nw_perproc
3090  integer, intent(out):: procb(nw,nkpt),proct(nw,0:nproc-1)
3091 !Local variables-------------------------------
3092  integer :: ikpt,iw,ir,ratio,m,iproc
3093  integer, allocatable :: proca(:,:)
3094 ! *********************************************************************
3095 
3096 
3097  proct(:,:)=0
3098  procb(:,:)=-1
3099 
3100  if(nproc.ge.2*nw) then
3101    ratio=nproc/nw
3102    ABI_ALLOCATE(proca,(nw,nproc/nw))
3103    do ir=1,ratio
3104      do iw=1,nw
3105        proca(iw,ir)=iw+(ir-1)*nw-1
3106        proct(iw,iw+(ir-1)*nw-1)=1
3107      enddo
3108    enddo
3109 !     do iw=1,nw
3110 !        write(6,*) " freq, procs"
3111 !        write(6,*) iw,proca(iw,:)
3112 !     enddo
3113    do iw=1,nw
3114      do ikpt=1,nkpt
3115        if(ikpt.le.ratio) procb(iw,ikpt)=proca(iw,ikpt)
3116        if(ikpt.ge.ratio) then
3117          m=mod(ikpt,ratio)
3118          if(m==0) then
3119             procb(iw,ikpt)=proca(iw,ratio)
3120          else
3121             procb(iw,ikpt)=proca(iw,m)
3122          endif
3123        endif
3124 !       write(6,*) " freq, k-point, proc"
3125 !       write(6,*) iw,ikpt,procb(iw,ikpt)
3126      enddo
3127    enddo
3128    nw_perproc=1
3129    ABI_DEALLOCATE(proca)
3130 
3131  else if (nproc.ge.nw) then
3132    do iw=1,nw
3133      procb(iw,:)= iw
3134      proct(iw,iw)=1
3135    enddo
3136 !  do iw=1,nw
3137 !        write(6,*) " freq, proc"
3138 !        write(6,*) iw, procb(iw,1)
3139 !  enddo
3140    nw_perproc=1
3141 
3142  else if (nproc.le.nw) then
3143    ratio=nw/nproc
3144    do iproc=0,nproc-1
3145      do iw=1,nw
3146        if (mod(iw-1,nproc)==iproc) then
3147          procb(iw,:)=iproc
3148          proct(iw,iproc)=1
3149        endif
3150      enddo
3151    enddo
3152    nw_perproc=ratio+1
3153 !  some procs will compute a number of frequency which is ratio and some
3154 !  other will compute a number of frequency which is ratio+1.
3155 
3156 !  do iw=1,nw
3157 !        write(6,*) " freq, proc"
3158 !        write(6,*) iw, procb(iw,1)
3159 !  enddo
3160   endif
3161 
3162 !     do iw=1,nw
3163 !      do iproc=0,nproc-1
3164 !        write(6,*) " freq, procs,?"
3165 !        write(6,*) iw,iproc,proct(iw,iproc)
3166 !      enddo
3167 !     enddo
3168 
3169 
3170 
3171  end subroutine distrib_paral

m_green/fermi_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 fermi_green

FUNCTION

  Compute Fermi level for DMFT or LDA.

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

 fermie: current value
 f_precision: precision of f required
 ITRLST: =1 if last iteration of DMFT
 opt_noninter   : if one wants the LDA fermi level
 max_iter : max number of iterations.

OUTPUT

 fermie: output value

PARENTS

      m_dmft

CHILDREN

      compute_green,integrate_green,newton,wrtout

SOURCE

3316 subroutine fermi_green(cryst_struc,green,paw_dmft,pawang,self)
3317 
3318  use m_abicore
3319 
3320  use defs_basis
3321  use defs_abitypes
3322  use m_pawang, only : pawang_type
3323  use m_crystal, only : crystal_t
3324  use m_paw_dmft, only: paw_dmft_type
3325  use m_self, only : self_type
3326 
3327 !This section has been created automatically by the script Abilint (TD).
3328 !Do not modify the following lines by hand.
3329 #undef ABI_FUNC
3330 #define ABI_FUNC 'fermi_green'
3331 !End of the abilint section
3332 
3333  implicit none
3334 
3335 !Arguments ------------------------------------
3336 !scalars
3337  type(crystal_t),intent(in) :: cryst_struc
3338  type(green_type),intent(inout) :: green
3339  type(paw_dmft_type), intent(inout) :: paw_dmft
3340  !type(MPI_type), intent(in) :: mpi_enreg
3341  type(pawang_type),intent(in) :: pawang
3342  type(self_type), intent(inout) :: self
3343 
3344 !Local variables-------------------------------
3345  integer :: ierr_hh,opt_noninter,max_iter
3346  real(dp):: x_precision,f_precision,fermi_old
3347 ! real(dp) :: hx
3348  character(len=500) :: message
3349 !************************************************************************
3350 !
3351  write(message,'(a,8x,a)') ch10,"  == Compute Fermi level"
3352  call wrtout(std_out,message,'COLL')
3353 
3354 !=============
3355 !headers
3356 !=============
3357  write(message,'(2a)') ch10, "  |---Newton method to search Fermi level ------------|"
3358  call wrtout(std_out,message,'COLL')
3359  write(message,'(2a,f13.6)') ch10, "  |--- Initial value for Fermi level",paw_dmft%fermie
3360  call wrtout(std_out,message,'COLL')
3361 
3362 !========================================
3363 !Define precision and nb of iterations
3364 !=========================================
3365  fermi_old=paw_dmft%fermie
3366  ierr_hh=0
3367  f_precision=paw_dmft%dmft_chpr
3368  !f_precision=0.01
3369  x_precision=tol5
3370 !if(option==1) then
3371 !f_precision=(erreursurlacharge)/100d0
3372 !else
3373 !f_precision=tol11
3374 !endif
3375  max_iter=1 ! for tests only
3376  !write(6,*) "for tests max_iter=1"
3377  max_iter=50
3378  opt_noninter=4
3379 
3380 !=====================
3381 !Call newton method
3382 !=====================
3383  write(message,'(a,4x,a,e13.6)') ch10," Precision required :",f_precision
3384  call wrtout(std_out,message,'COLL')
3385  if (f_precision<10_dp)  then
3386    call newton(cryst_struc,green,paw_dmft,pawang,self,&
3387 &   paw_dmft%fermie,x_precision,max_iter,f_precision,ierr_hh,opt_noninter)
3388  end if
3389 
3390 !===========================
3391 !Deals with errors signals
3392 !===========================
3393  if(ierr_hh==-314) then
3394    write(message,'(a)') "Warning, check Fermi level"
3395    call wrtout(std_out,message,'COLL')
3396 !  call abi_abort('COLL')
3397    write(message,'(2a,f13.6)') ch10, "  |---  Final  value for Fermi level (check)",paw_dmft%fermie
3398    call wrtout(std_out,message,'COLL')
3399  else if (ierr_hh==-123) then
3400    write(message,'(a,f13.6)') " Fermi level is put to",fermi_old
3401    paw_dmft%fermie=fermi_old
3402    call wrtout(std_out,message,'COLL')
3403 
3404 !  =====================================
3405 !  If fermi level search was successful
3406 !  =====================================
3407  else
3408    write(message,'(a,4x,a,e13.6)') ch10," Precision achieved on Fermi Level :",x_precision
3409    call wrtout(std_out,message,'COLL')
3410    write(message,'(4x,a,e13.6)') " Precision achieved on number of electrons :",f_precision
3411    call wrtout(std_out,message,'COLL')
3412    write(message,'(2a,f13.6)') ch10, "  |---  Final  value for Fermi level",paw_dmft%fermie
3413    call wrtout(std_out,message,'COLL')
3414  end if
3415 
3416 !========================================================
3417 !Check convergence of fermi level during DMFT iterations
3418 !========================================================
3419  if(paw_dmft%idmftloop>=2) then
3420    if(abs(paw_dmft%fermie-fermi_old).le.paw_dmft%dmft_fepr) then
3421 !    write(message,'(a,8x,a,e9.2,a,8x,a,e12.5)') ch10,"|fermie(n)-fermie(n-1)|=<",paw_dmft%dmft_fepr,ch10,&
3422      write(message,'(a,8x,a,e9.2,a,e9.2,a,8x,a,e12.5)') ch10,"|fermie(n)-fermie(n-1)|=",&
3423 &     abs(paw_dmft%fermie-fermi_old),"<",paw_dmft%dmft_fepr,ch10,&
3424 &     "=> DMFT Loop: Fermi level is converged to:",paw_dmft%fermie
3425      call wrtout(std_out,message,'COLL')
3426      green%ifermie_cv=1
3427    else
3428      write(message,'(a,8x,a,2f12.5)') ch10,"DMFT Loop: Fermi level is not converged:",&
3429 &     paw_dmft%fermie
3430      call wrtout(std_out,message,'COLL')
3431      green%ifermie_cv=0
3432    end if
3433  end if
3434  write(message,'(2a)') ch10, "  |---------------------------------------------------|"
3435  call wrtout(std_out,message,'COLL')
3436 !
3437 
3438 !==========================================================
3439 !Recompute full green function including non diag elements
3440 !==========================================================
3441  call compute_green(cryst_struc,green,paw_dmft,pawang,0,self,opt_self=1,opt_nonxsum=1)
3442  call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=0,opt_ksloc=3) !,opt_nonxsum=1)
3443 
3444  return
3445 end subroutine fermi_green

m_green/fourier_fct [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 fourier_fct

FUNCTION

 Do fourier transformation from matsubara space to imaginary time
 (A spline is performed )

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

  ldiag    = option according to diagonal or non-diagonal elements
  opt_four = option for direct (1) or inverse (-1) fourier transform
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data

OUTPUT

SIDE EFFECTS

  fw= function is frequency space
  ft= function is time space

PARENTS

      local_ks_green,m_green

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2730 subroutine fourier_fct(fw,ft,ldiag,ltau,opt_four,paw_dmft)
2731 
2732  use m_paw_dmft, only : paw_dmft_type, construct_nwli_dmft
2733  use m_splines
2734 
2735 !This section has been created automatically by the script Abilint (TD).
2736 !Do not modify the following lines by hand.
2737 #undef ABI_FUNC
2738 #define ABI_FUNC 'fourier_fct'
2739 !End of the abilint section
2740 
2741  implicit none
2742 
2743 !Arguments ------------------------------------
2744 !type
2745  logical,intent(in) :: ldiag
2746  integer,intent(in) :: ltau,opt_four
2747  type(paw_dmft_type), intent(in) :: paw_dmft
2748  complex(dpc), intent(inout) :: fw(paw_dmft%dmft_nwlo)
2749  complex(dpc), intent(inout) :: ft(ltau)
2750 
2751 !local variables-------------------------------
2752  complex(dpc), allocatable ::  splined_li(:)
2753  complex(dpc), allocatable ::  tospline_li(:)
2754 ! complex(dpc), allocatable :: fw1(:)
2755  real(dp), allocatable :: ftr(:)
2756  real(dp) :: beta
2757  complex(dpc) :: xsto
2758  integer :: iflag,ifreq,itau,iwarn,log_direct
2759  character(len=500) :: message
2760  real(dp), allocatable :: omega_li(:)
2761 ! *********************************************************************
2762  beta=one/paw_dmft%temp
2763  iflag=0
2764  log_direct=1
2765 
2766  if(ldiag) iflag=1
2767 
2768 ! == inverse fourier transform
2769  if(opt_four==-1) then
2770 
2771    ABI_ALLOCATE(splined_li,(paw_dmft%dmft_nwli))
2772 !   allocate(fw1(0:paw_dmft%dmft_nwlo-1))
2773    if(paw_dmft%dmft_log_freq==1) then
2774      ABI_ALLOCATE(omega_li,(1:paw_dmft%dmft_nwli))
2775      call construct_nwli_dmft(paw_dmft,paw_dmft%dmft_nwli,omega_li)
2776      call spline_c(paw_dmft%dmft_nwlo,paw_dmft%dmft_nwli,paw_dmft%omega_lo,&
2777 &                   omega_li,splined_li,fw)
2778      ABI_DEALLOCATE(omega_li)
2779    else
2780      splined_li=fw
2781    endif
2782    call invfourier(splined_li,ft,paw_dmft%dmft_nwli,ltau,iflag,beta)
2783 !   deallocate(fw1)
2784    ABI_DEALLOCATE(splined_li)
2785 
2786 ! == direct fourier transform
2787  else if(opt_four==1) then
2788 
2789    ABI_ALLOCATE(ftr,(ltau))
2790 
2791    iwarn=0
2792    do itau=1,ltau
2793      if(abs(aimag(ft(itau)))>tol12) then
2794        if(ldiag) then
2795          write(message,'(a,a,2(e15.4))') ch10,&
2796 &          "green function is not real in imaginary time space",ft(itau)
2797          MSG_ERROR(message)
2798        else
2799          iwarn=iwarn+1
2800          ftr(itau)=real(ft(itau))
2801          xsto=ft(itau)
2802        endif
2803      else
2804        ftr(itau)=real(ft(itau))
2805      endif
2806 !       write(std_out,*) itau,ftr(itau)
2807    enddo
2808 
2809    ABI_ALLOCATE(tospline_li,(paw_dmft%dmft_nwli))
2810    if (log_direct==1) then
2811 !   do not have a physical meaning..because the log frequency is not
2812 !   one of the linear frequency.
2813 !   if log frequencies are also one of linear frequency, it should be good
2814      do ifreq=1, paw_dmft%dmft_nwlo
2815        call nfourier2(ftr,fw(ifreq),iflag,paw_dmft%omega_lo(ifreq),ltau,beta)
2816      enddo
2817    else
2818      call nfourier(ftr,tospline_li,iflag,paw_dmft%dmft_nwli-1,ltau,beta)
2819      do ifreq=1,paw_dmft%dmft_nwli
2820         write(1112,*) paw_dmft%temp*pi*real(2*ifreq-1,kind=dp),real(tospline_li(ifreq),kind=dp),aimag(tospline_li(ifreq))
2821         !write(1112,*) paw_dmft%omega_li(ifreq),real(tospline_li(ifreq)),aimag(tospline_li(ifreq))
2822      enddo
2823      if(paw_dmft%dmft_log_freq==1) then
2824        ABI_ALLOCATE(omega_li,(1:paw_dmft%dmft_nwli))
2825        call construct_nwli_dmft(paw_dmft,paw_dmft%dmft_nwli,omega_li)
2826        call spline_c(paw_dmft%dmft_nwli,paw_dmft%dmft_nwlo,omega_li,&
2827 &                 paw_dmft%omega_lo,fw,tospline_li)
2828        ABI_DEALLOCATE(omega_li)
2829      else
2830        fw=tospline_li
2831      endif
2832    endif
2833 
2834    ABI_DEALLOCATE(tospline_li)
2835 
2836    ABI_DEALLOCATE(ftr)
2837    if(iwarn>0) then
2838      write(message,'(a,a,2(e15.4))') ch10,&
2839 &     "WARNING: off-diag green function is not real in imaginary time space",xsto
2840      call wrtout(std_out,message,'COLL')
2841    endif
2842 
2843  endif
2844 
2845 end subroutine fourier_fct

m_green/fourier_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 fourier_green

FUNCTION

  integrate green function in the band index basis

INPUTS

  cryst_struc <type(crystal_t)>= crystal structure data.
  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawang <type(pawang)>=paw angular mesh and related data

OUTPUT

PARENTS

      impurity_solve,m_green,qmc_prep_ctqmc

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2013 subroutine fourier_green(cryst_struc,green,paw_dmft,pawang,opt_ksloc,opt_tw)
2014 
2015  use m_pawang, only : pawang_type
2016  use m_crystal, only : crystal_t
2017  use m_matlu, only : sym_matlu,init_matlu,destroy_matlu,print_matlu
2018  use m_paw_dmft, only : paw_dmft_type
2019  use m_oper, only : loc_oper
2020 
2021 !This section has been created automatically by the script Abilint (TD).
2022 !Do not modify the following lines by hand.
2023 #undef ABI_FUNC
2024 #define ABI_FUNC 'fourier_green'
2025 !End of the abilint section
2026 
2027  implicit none
2028 
2029 !Arguments ------------------------------------
2030 !type
2031  type(crystal_t),intent(in) :: cryst_struc
2032  type(green_type),intent(inout) :: green
2033  !type(MPI_type), intent(in) :: mpi_enreg
2034  type(pawang_type),intent(in) :: pawang
2035  type(paw_dmft_type), intent(in) :: paw_dmft
2036  integer,intent(in) :: opt_ksloc,opt_tw ! fourier on ks or local
2037 
2038 !local variables-------------------------------
2039  integer :: iatom,ib,ib1,ierr,ifreq,ikpt,im,im1,iparal,is,ispinor,ispinor1,itau
2040  integer :: mband,mbandc,myproc,natom,ndim,nkpt,nproc,nspinor,nsppol,spacecomm!,opt_four
2041  character(len=500) :: message
2042 ! complex(dpc):: ybcbeg,ybcend
2043 ! arrays
2044  complex(dpc), allocatable :: fw(:)
2045  complex(dpc), allocatable :: ft(:)
2046  type(green_type) :: green_temp
2047 ! *********************************************************************
2048 ! ybcbeg=czero
2049 ! ybcend=czero
2050 
2051 ! define spaceComm, myproc, and nproc from world communicator
2052 ! and mpi_enreg
2053  spacecomm=paw_dmft%spacecomm
2054  myproc=paw_dmft%myproc
2055  nproc=paw_dmft%nproc
2056 
2057 ! Initialise integers
2058  mband   = paw_dmft%mband
2059  mbandc  = paw_dmft%mbandc
2060  nkpt    = paw_dmft%nkpt
2061  nspinor = paw_dmft%nspinor
2062  nsppol  = paw_dmft%nsppol
2063  natom   = cryst_struc%natom
2064 
2065 ! Only imaginary frequencies here
2066  if(green%w_type=="real") then
2067    message = 'fourier_green not implemented for real frequency'
2068    MSG_BUG(message)
2069  endif
2070 
2071 ! Initialise temporary green function
2072  call init_green(green_temp,paw_dmft)
2073  call init_green_tau(green_temp,paw_dmft)
2074 
2075  !green%oper(:)%matlu(1)%mat(1,1,1,1,1)
2076  ABI_ALLOCATE(fw,(green%nw))
2077  ABI_ALLOCATE(ft,(green%dmftqmc_l))
2078 
2079 !==============================================
2080 ! == Inverse fourier transformation ===========
2081 !==============================================
2082 
2083  if(opt_tw==-1) then
2084 
2085 !  = For Kohn Sham green function
2086 !==============================================
2087    if(opt_ksloc==1.and.green%use_oper_tau_ks==1) then
2088 
2089      do is = 1 , nsppol
2090        do ikpt = 1, nkpt
2091          do ib = 1, mbandc
2092            do ib1 = 1, mbandc
2093              do ifreq=1,green%nw
2094                fw(ifreq)=green%oper(ifreq)%ks(is,ikpt,ib,ib1)
2095              enddo
2096              call fourier_fct(fw,ft,ib==ib1,green%dmftqmc_l,-1,paw_dmft) ! inverse fourier
2097              do itau=1,green%dmftqmc_l
2098                green_temp%oper_tau(itau)%ks(is,ikpt,ib,ib1)=ft(itau)
2099              enddo
2100              if(ib==ib1) then
2101                green%occup_tau%ks(is,ikpt,ib,ib1)=ft(1)+one
2102              else
2103                green%occup_tau%ks(is,ikpt,ib,ib1)=ft(1)
2104              endif
2105            enddo ! ib1
2106          enddo ! ib
2107        enddo ! ikpt
2108      enddo ! isppol
2109 
2110 !  = Post-treatment: necessary in the case of nspinor==2, but valid anywhere
2111 !    because G(tau)=G_{nu,nu'}(tau)+[G_{nu,nu'}(tau)]*
2112 
2113      do is = 1 , nsppol
2114        do ikpt = 1, nkpt
2115          do ib = 1, mbandc
2116            do ib1 = 1, mbandc
2117              do itau=1,green%dmftqmc_l
2118                green%oper_tau(itau)%ks(is,ikpt,ib,ib1)=&
2119 &                (green_temp%oper_tau(itau)%ks(is,ikpt,ib,ib1)+ &
2120 &                 conjg(green_temp%oper_tau(itau)%ks(is,ikpt,ib1,ib)))/two
2121                if(ib==ib1) then
2122                  green%occup_tau%ks(is,ikpt,ib,ib1)=green%oper_tau(1)%ks(is,ikpt,ib,ib1)+one
2123                else
2124                  green%occup_tau%ks(is,ikpt,ib,ib1)=green%oper_tau(1)%ks(is,ikpt,ib,ib1)
2125                endif
2126              enddo
2127            enddo ! ib1
2128          enddo ! ib
2129        enddo ! ikpt
2130      enddo ! isppol
2131      call loc_oper(green%occup_tau,paw_dmft,1)
2132      call sym_matlu(cryst_struc,green%occup_tau%matlu,pawang)
2133      write(message,'(a,a,i10,a)') ch10,"  green%occup_tau%matlu from green_occup_tau%ks"
2134      call wrtout(std_out,message,'COLL')
2135      call print_matlu(green%occup_tau%matlu,natom,prtopt=3)
2136    endif
2137 
2138 !  = For local green function
2139 !==============================================
2140    if(opt_ksloc ==2) then
2141 
2142      iparal=0
2143      do iatom=1, natom
2144        if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
2145          ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
2146          do is = 1 , nsppol
2147            do ispinor = 1, nspinor
2148              do ispinor1 = 1, nspinor
2149                do im=1,ndim
2150                  do im1=1,ndim
2151                    iparal=iparal+1
2152                    if(mod(iparal-1,nproc)==myproc) then
2153                      do ifreq=1,green%nw
2154                        fw(ifreq)=green%oper(ifreq)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2155                      enddo
2156                      ! inverse fourier
2157 !                     write(std_out,'(a)') "fourierbeforeposttreatement,ispinor,ispinor1,is,im,im1"
2158 !                     write(std_out,'(a,5i4,f12.5,f12.5)') "fourier",ispinor,ispinor1,is,im,im1
2159 !                     write(std_out,'(a,e12.5,e12.5)')&
2160 !                     &"green%oper(4)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)"&
2161 !&                     ,green%oper(4)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2162                      call fourier_fct(fw,ft,(im==im1).and.(ispinor==ispinor1),green%dmftqmc_l,-1,paw_dmft)
2163                      do itau=1,green%dmftqmc_l
2164                        green_temp%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(itau)
2165 !                     write(std_out,*) itau,green_temp%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2166                      enddo
2167 !                    if((im==im1).and.(ispinor==ispinor1)) then
2168 !                      green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(1)+one
2169 !                    else
2170 !                      green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(1)
2171 !                    endif
2172                    endif ! iparal
2173                  enddo
2174                enddo
2175              enddo ! ispinor1
2176            enddo ! ispinor
2177          enddo ! isppol
2178        endif ! lpawu.ne.-1
2179      enddo ! iatom
2180 
2181 !    Parallelisation must be finished here, because in the post
2182 !    treatment, value from different proc will be mixed.
2183      call xmpi_barrier(spacecomm)
2184      do iatom=1,natom
2185        do itau=1,green%dmftqmc_l
2186         call xmpi_sum(green_temp%oper_tau(itau)%matlu(iatom)%mat,spacecomm,ierr)
2187        enddo
2188      enddo
2189      call xmpi_barrier(spacecomm)
2190 
2191 !  = Post-treatment: necessary in the case of nspinor==2, but valid anywhere
2192 !    because G(tau)=G_{LL'}^{sigma,sigma'}(tau)+[G_{L'L}^{sigma',sigma}(tau)]*
2193 
2194      if(nspinor>0) Then
2195        do iatom=1, natom
2196          if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
2197            ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
2198            do is = 1 , nsppol
2199              do ispinor = 1, nspinor
2200                do ispinor1 = 1, nspinor
2201                  do im=1,ndim
2202                    do im1=1,ndim
2203 !                   write(std_out,'(a,5i4,f12.5,f12.5)') "fourier -1",ispinor,ispinor1,is,im,im1
2204                      do itau=1,green%dmftqmc_l
2205                        green%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=&
2206 &                      ((green_temp%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)+ &
2207 &                      conjg(green_temp%oper_tau(itau)%matlu(iatom)%mat(im1,im,is,ispinor1,ispinor))))/two
2208 !                       write(std_out,*) itau,green%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2209                        if((im==im1).and.(ispinor==ispinor1)) then
2210                          green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=&
2211 &                         green%oper_tau(1)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)+one
2212                        else
2213                          green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=&
2214 &                         green%oper_tau(1)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2215                        endif ! test diag
2216                      enddo  ! itau
2217                    enddo  ! im1
2218                  enddo ! im
2219                enddo ! ispinor1
2220              enddo ! ispinor
2221            enddo ! isppol
2222          endif
2223        enddo ! iatom
2224      endif ! nspinor>0
2225    endif ! opt_ksloc=2
2226  endif ! opt_tw==-1
2227 
2228 
2229 !==============================================
2230 ! == Direct fourier transformation
2231 !==============================================
2232 ! todo_ba ft useful only for diagonal elements ...
2233 
2234  if(opt_tw==1) then
2235 
2236 !  = For local green function
2237 !==============================================
2238    if(opt_ksloc ==2) then
2239 
2240      iparal=0
2241      do iatom=1, natom
2242 !  put to zero (usefull at least for parallelism)
2243        do ifreq=1,green%nw
2244          green%oper(ifreq)%matlu(iatom)%mat=czero
2245        enddo
2246        if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
2247          ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
2248          do is = 1 , nsppol
2249            do ispinor = 1, nspinor
2250              do ispinor1 = 1, nspinor
2251                do im=1,ndim
2252                  do im1=1,ndim
2253                    iparal=iparal+1
2254                    if(mod(iparal-1,nproc)==myproc) then
2255                      do itau=1,green%dmftqmc_l
2256                        ft(itau)=green%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
2257                      enddo
2258 !                     if(im==1.and.im1==1) write(std_out,*) "ft(itau=1)",is,ft(1) ! debug
2259                      call fourier_fct(fw,ft,(im==im1).and.(ispinor==ispinor1),green%dmftqmc_l,1,paw_dmft)
2260                      do ifreq=1,green%nw
2261                        green%oper(ifreq)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=fw(ifreq)
2262                      enddo
2263 !                     if(im==1.and.im1==1) write(std_out,*) "fw(ifreq=1)",is,fw(1) ! debug
2264                    endif
2265                  enddo
2266                enddo
2267              enddo ! ispinor1
2268            enddo ! ispinor
2269          enddo ! isppol
2270        endif ! lpawu=/-1
2271      enddo ! iatom
2272      call xmpi_barrier(spacecomm)
2273      do iatom=1,natom
2274        do ifreq=1,green%nw
2275        call xmpi_sum(green%oper(ifreq)%matlu(iatom)%mat,spacecomm,ierr)
2276        enddo
2277      enddo
2278    endif ! opt_ksloc=2
2279  endif ! opt_tw==-1
2280  ABI_DEALLOCATE(fw)
2281  ABI_DEALLOCATE(ft)
2282  call destroy_green_tau(green_temp)
2283  call destroy_green(green_temp)
2284 
2285 end subroutine fourier_green

m_green/green_type [ Types ]

[ Top ] [ m_green ] [ Types ]

NAME

  green_type

FUNCTION

  This structured datatype contains the necessary data

SOURCE

 80  type, public :: green_type ! for each atom
 81 
 82   integer :: dmft_nwlo
 83   ! dmft frequencies
 84 
 85   character(len=4) :: w_type
 86   ! type of frequencies used
 87 
 88   character(len=12) :: whichgreen
 89   ! describe the type of green function computed (LDA, DMFT, KS..)
 90 
 91   integer :: nw
 92   ! dmft frequencies
 93 
 94   integer :: fileprt_w
 95   ! 1 if file created
 96 
 97   integer :: fileprt_tau
 98   ! 1 if file created
 99 
100   integer :: dmft_nwli
101   ! dmft frequencies
102 
103   integer :: dmftqmc_l
104   ! number of time slices for QMC
105 
106   integer :: ifermie_cv
107 
108   integer :: ichargeloc_cv
109 
110   integer :: use_oper_tau_ks
111   ! 0 do not use oper_tau_ks
112   ! 1 use oper_tau_ks
113 
114   real(dp) :: charge_ks
115   ! Total charge computed from ks orbitals
116 
117   integer :: has_charge_matlu_solver
118   ! =0 charge_matlu_solver not allocated
119   ! =1 charge_matlu_solver is allocated
120   ! =2 charge_matlu_solver is calculated (ie calculation of LOCAL CORRELATED occupations is done  from
121   ! solver green function)
122 
123   integer :: has_greenmatlu_xsum
124   ! =1 green%oper%matlu xsum in compute_green
125   ! =0 green%oper%matlu non xsumed in compute_green
126   ! used in integrate_green to checked that green function was computed in
127   ! integrate_green.
128 
129   integer :: has_charge_matlu
130   ! =2 if calculation of LOCAL CORRELATED occupations is done  from
131 
132   integer :: has_charge_matlu_prev
133   ! =0 charge_matlu_prev not allocated
134   ! =1 charge_matlu_prev is allocated
135   ! =2 charge_matlu_prev is calculated (ie calculation of LOCAL CORRELATED occupations is done  from
136   ! solver green function)
137 
138   integer, allocatable :: procb(:,:)
139 
140   integer, allocatable :: proct(:,:)
141 
142   real(dp), allocatable :: charge_matlu_prev(:,:)
143   ! Total charge on correlated orbitals from previous iteration
144 
145   real(dp), allocatable :: charge_matlu(:,:)
146   ! Total charge on correlated orbitals
147 ! todo_ba name of charge_matlu is misleading: should be changed
148 
149   real(dp), allocatable :: charge_matlu_solver(:,:)
150   ! Total charge on correlated orbitals obtained from solver by
151   ! integration over frequencies.
152 
153   real(dp), allocatable :: tau(:)
154   ! value of time in imaginary space
155 
156   real(dp), pointer :: omega(:) => null()
157   ! value of frequencies
158 
159   real(dp), allocatable :: ecorr_qmc(:)
160   ! Correlation energy for a given atom in qmc
161 
162   type(oper_type), allocatable :: oper(:)
163   ! green function  in different basis
164 
165   type(oper_type), allocatable :: oper_tau(:)
166   ! green function  in different basis
167 
168   type(oper_type) :: occup
169   ! occupation in different basis
170 
171   type(oper_type) :: occup_tau
172   ! occupation in different basis
173 
174  end type green_type
175 
176 !----------------------------------------------------------------------
177 
178 
179 CONTAINS

m_green/icip_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

  icip_green

FUNCTION

  init, compute, integrate and print lda green function

INPUTS

  char1 = character which precises the type of green function computed.
  cryst_struc <type(crystal_t)>= crystal structure data.
  mpi_enreg=informations about MPI parallelization
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  green  <type(green_type)>= green function data
  pawang <type(pawang)>=paw angular mesh and related data
  pawprtvol  = option for printing
  self <type(self_type)>= variables related to self-energy
  opt_self = optional argument, if =1, upfold self-energy

OUTPUT

PARENTS

      m_dmft

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

1901 subroutine icip_green(char1,cryst_struc,green,paw_dmft,pawang,pawprtvol,self,opt_self)
1902 
1903  use m_pawang, only : pawang_type
1904  use m_crystal, only : crystal_t
1905  use m_matlu, only : sym_matlu
1906  use m_paw_dmft, only : paw_dmft_type
1907  use m_oper, only : loc_oper
1908  use m_self, only : self_type
1909 
1910 !This section has been created automatically by the script Abilint (TD).
1911 !Do not modify the following lines by hand.
1912 #undef ABI_FUNC
1913 #define ABI_FUNC 'icip_green'
1914 !End of the abilint section
1915 
1916  implicit none
1917 
1918 !Arguments ------------------------------------
1919 !type
1920  type(crystal_t),intent(in) :: cryst_struc
1921  !type(MPI_type), intent(in) :: mpi_enreg
1922  type(green_type),intent(inout) :: green
1923  type(pawang_type),intent(in) :: pawang
1924  type(paw_dmft_type), intent(inout) :: paw_dmft
1925  type(self_type), intent(inout) :: self
1926  integer, intent(in) :: pawprtvol
1927  character(len=*), intent(in) :: char1
1928  integer, optional, intent(in) :: opt_self
1929 !Local variables ------------------------------------
1930  integer :: option,optself,optlocks,prtopt_for_integrate_green,opt_nonxsum_icip
1931  character(len=500) :: message
1932 ! *********************************************************************
1933  green%whichgreen="LDA"
1934  prtopt_for_integrate_green=2
1935 
1936  if(present(opt_self)) then
1937    optself=opt_self
1938  else
1939    optself=0
1940  endif
1941  opt_nonxsum_icip=1
1942  if(paw_dmft%dmftcheck==1) then ! for fourier_green
1943  !  opt_nonxsum_icip=0
1944  endif
1945 
1946  call init_green(green,paw_dmft)
1947 
1948 ! useless test ok
1949 ! call printocc_green(green,1,paw_dmft,2)
1950 ! write(std_out,*)" printocc_green zero finished "
1951 
1952 !== Compute green%oper(:)%ks
1953 !== Deduce  green%oper(:)%matlu(:)%mat
1954  call compute_green(cryst_struc,green,paw_dmft,&
1955 & pawang,pawprtvol,self,optself,opt_nonxsum=opt_nonxsum_icip)
1956  if(paw_dmft%dmft_prgn>=1.and.paw_dmft%dmft_prgn<=2) then
1957    optlocks=paw_dmft%dmft_prgn*2+1 ! if dmft_prgn==2 => do not print
1958    if(paw_dmft%lpsichiortho==1.and.pawprtvol>-100)  then
1959      call print_green(char1,green,optlocks,paw_dmft,pawprtvol)
1960 !     call print_green('inicip',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
1961    endif
1962  endif
1963 
1964 !== Integrate green%oper(:)%ks
1965 !== Integrate green%oper(:)%matlu(:)%mat
1966  call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt_for_integrate_green,opt_ksloc=3)!,opt_nonxsum=opt_nonxsum_icip)
1967 !== Print green%oper(:)%ks
1968 !== Print green%oper(:)%matlu(:)%mat
1969  if(char1=="LDA") then
1970    option=1
1971    if(self%oper(1)%matlu(1)%lpawu/=-1) then
1972      if(abs(real(self%oper(1)%matlu(1)%mat(1,1,1,1,1)))>tol7) then
1973 ! todo_ab: generalise this
1974        write(message,'(a,a,2(e15.4))') ch10,&
1975 &        "Warning:  a LDA calculation is carried out and self is not zero"
1976        call wrtout(std_out,message,'COLL')
1977 !       call abi_abort('COLL')
1978      endif
1979    endif
1980  else
1981    option=5
1982  endif
1983  if(pawprtvol>-100) then
1984    call printocc_green(green,option,paw_dmft,3,chtype=char1)
1985  end if
1986 
1987 end subroutine icip_green

m_green/init_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 init_green

FUNCTION

  Allocate variables used in type green_type.

INPUTS

  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  opt_oper_ksloc (optional) = option for init_oper
  wtype = "real" Green function will be computed for real frequencies
        = "imag" Green function will be computed for imaginary frequencies

 OUTPUTS
 green  = variable of type green_type

PARENTS

      m_dmft,dyson,m_hubbard_one,m_green,qmc_prep_ctqmc,spectral_function

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

207 subroutine init_green(green,paw_dmft,opt_oper_ksloc,wtype)
208 
209  use m_crystal, only : crystal_t
210  use m_oper, only : init_oper
211  use m_paw_dmft, only: paw_dmft_type
212 
213 !This section has been created automatically by the script Abilint (TD).
214 !Do not modify the following lines by hand.
215 #undef ABI_FUNC
216 #define ABI_FUNC 'init_green'
217 !End of the abilint section
218 
219  implicit none
220 
221 !Arguments ------------------------------------
222 !scalars
223 !type
224  type(green_type), intent(out) :: green
225  type(paw_dmft_type), intent(in) :: paw_dmft
226  integer, optional, intent(in) :: opt_oper_ksloc
227  character(len=4), optional :: wtype
228 !local variables ------------------------------------
229  integer :: ifreq,nw,optoper_ksloc,nw_perproc
230 !************************************************************************
231  if(present(opt_oper_ksloc)) then
232    optoper_ksloc=opt_oper_ksloc
233  else
234    optoper_ksloc=2
235  endif
236  if(present(wtype)) then
237    green%w_type=wtype
238  else
239    green%w_type="imag"
240  endif
241 
242  if(green%w_type=="imag") then
243    nw=paw_dmft%dmft_nwlo
244    green%omega=>paw_dmft%omega_lo
245  else if(green%w_type=="real") then
246    nw=2*paw_dmft%dmft_nwr
247    green%omega=>paw_dmft%omega_r
248  endif
249 
250  green%dmft_nwlo=paw_dmft%dmft_nwlo
251  green%dmft_nwli=paw_dmft%dmft_nwli
252  green%charge_ks=zero
253  green%has_charge_matlu=0
254  ABI_ALLOCATE(green%charge_matlu,(paw_dmft%natom,paw_dmft%nsppol+1))
255  green%charge_matlu=zero
256  green%has_charge_matlu=1
257  green%has_greenmatlu_xsum=0
258 
259  green%has_charge_matlu_solver=0
260  ABI_ALLOCATE(green%charge_matlu_solver,(paw_dmft%natom,paw_dmft%nsppol+1))
261  green%charge_matlu_solver=zero
262  green%has_charge_matlu_solver=1
263 
264  green%has_charge_matlu_prev=0
265  ABI_ALLOCATE(green%charge_matlu_prev,(paw_dmft%natom,paw_dmft%nsppol+1))
266  green%charge_matlu_prev=zero
267  green%has_charge_matlu_prev=1
268 
269  call init_oper(paw_dmft,green%occup,opt_ksloc=3)
270 
271 !  built simple arrays to distribute the tasks in compute_green.
272  ABI_ALLOCATE(green%procb,(nw,paw_dmft%nkpt))
273  ABI_ALLOCATE(green%proct,(nw,0:paw_dmft%nproc-1))
274 
275  call distrib_paral(paw_dmft%nkpt,paw_dmft%nproc,nw,nw_perproc,green%procb,green%proct)
276  green%nw=nw
277 
278 !  need to distribute memory over frequencies
279 
280 !!  begin of temporary modificatios
281 ! ABI_DATATYPE_ALLOCATE(green%oper,(green%nw_perproc))
282 ! do ifreq=1,green%nw_perproc
283 !  call init_oper(paw_dmft,green%oper(ifreq),opt_ksloc=optoper_ksloc)
284 ! enddo
285 !
286 ! do ifreq=1,green%nw
287 !   if(green%proct(ifreq,myproc)==1) then
288 !     do ikpt = 1 , paw_dmft%nkpt
289 !       if (green%procb(ifreq,ikpt)==myproc) then
290 !       endif
291 !     enddo ! ikpt
292 !   endif ! parallelisation
293 ! enddo ! ifreq
294 !
295 !
296 !!  end of temporary modificatios
297  ABI_DATATYPE_ALLOCATE(green%oper,(green%nw))
298  do ifreq=1,green%nw
299   call init_oper(paw_dmft,green%oper(ifreq),opt_ksloc=optoper_ksloc)
300  enddo
301  green%ifermie_cv=0
302  green%ichargeloc_cv=0
303 
304  if(paw_dmft%dmft_solv>=4) then
305    ABI_ALLOCATE(green%ecorr_qmc,(paw_dmft%natom))
306  end if
307  if(paw_dmft%dmft_solv>=4) green%ecorr_qmc(:)=zero
308 
309 
310  green%fileprt_w=0
311  green%fileprt_tau=0
312 
313 
314 end subroutine init_green

m_green/init_green_tau [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 init_green_tau

FUNCTION

  Allocate variables used in type green_type.

INPUTS

  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data

 OUTPUTS
  green  <type(green_type)>= green function data

PARENTS

      impurity_solve,m_green

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

339 subroutine init_green_tau(green,paw_dmft,opt_ksloc)
340 
341  use m_oper, only : init_oper
342  use m_paw_dmft, only: paw_dmft_type
343 
344 !This section has been created automatically by the script Abilint (TD).
345 !Do not modify the following lines by hand.
346 #undef ABI_FUNC
347 #define ABI_FUNC 'init_green_tau'
348 !End of the abilint section
349 
350  implicit none
351 
352 !Arguments ------------------------------------
353 !scalars
354 !type
355  type(green_type), intent(inout) :: green
356  type(paw_dmft_type), intent(in) :: paw_dmft
357  integer, optional, intent(in) :: opt_ksloc
358 !local variables ------------------------------------
359  integer :: optksloc
360  integer :: itau
361 !************************************************************************
362  if(present(opt_ksloc)) then
363    optksloc=opt_ksloc
364  else
365    optksloc=3
366  endif
367  green%use_oper_tau_ks=0
368  if(green%use_oper_tau_ks==0) then
369    optksloc=2
370  endif
371 
372  green%dmftqmc_l=paw_dmft%dmftqmc_l
373  ABI_ALLOCATE(green%tau,(green%dmftqmc_l))
374  do itau=1,green%dmftqmc_l
375   green%tau(itau)=float(itau-1)/float(green%dmftqmc_l)/paw_dmft%temp
376  enddo
377 
378  call init_oper(paw_dmft,green%occup_tau,optksloc)
379 
380  ABI_DATATYPE_ALLOCATE(green%oper_tau,(paw_dmft%dmftqmc_l))
381  do itau=1,green%dmftqmc_l
382   call init_oper(paw_dmft,green%oper_tau(itau),optksloc)
383  enddo
384 
385 end subroutine init_green_tau

m_green/int_fct [ Modules ]

[ Top ] [ m_green ] [ Modules ]

NAME

 int_fct

FUNCTION

 Do integration in matsubara space

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

  ff= function is frequency space
  ldiag    = option according to diagonal or non-diagonal elements
  option = nspinor
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  proct= for parallelism

OUTPUT

  integral = integral of ff over matsubara frequencies

SIDE EFFECTS

  ft= function is time space

PARENTS

      m_green,qmc_prep_ctqmc

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2594 subroutine int_fct(ff,ldiag,option,paw_dmft,integral,procb,myproc)
2595 
2596  use m_paw_dmft, only : paw_dmft_type
2597 
2598 !This section has been created automatically by the script Abilint (TD).
2599 !Do not modify the following lines by hand.
2600 #undef ABI_FUNC
2601 #define ABI_FUNC 'int_fct'
2602 !End of the abilint section
2603 
2604  implicit none
2605 
2606 !Arguments ------------------------------------
2607 !type
2608  logical,intent(in) :: ldiag
2609  integer,intent(in) :: option
2610  complex(dpc),intent(out) :: integral
2611  type(paw_dmft_type), intent(in) :: paw_dmft
2612  complex(dpc), intent(in) :: ff(paw_dmft%dmft_nwlo)
2613  integer, optional, intent(in) :: procb(paw_dmft%dmft_nwlo)
2614  integer, optional, intent(in) :: myproc
2615 
2616 !local variables-------------------------------
2617  logical, allocatable :: procb2(:)
2618  character(len=500) :: message
2619  integer :: ifreq
2620 ! *********************************************************************
2621  ABI_ALLOCATE(procb2,(paw_dmft%dmft_nwlo))
2622  if(present(procb).and.present(myproc)) then
2623    do ifreq=1,paw_dmft%dmft_nwlo
2624     procb2(ifreq)=(procb(ifreq)==myproc)
2625    enddo
2626  else if(present(procb).and..not.present(myproc)) then
2627    write(message,'(a,a,2(e15.4))') ch10,&
2628 &    "BUG: procb is present and not myproc in int_fct"
2629    MSG_BUG(message)
2630  else if(.not.present(procb).and.present(myproc)) then
2631    write(message,'(a,a,2(e15.4))') ch10,&
2632 &    "BUG: procb is not present and myproc is in int_fct"
2633    MSG_BUG(message)
2634  else
2635    do ifreq=1,paw_dmft%dmft_nwlo
2636      procb2(ifreq)=(1==1)
2637    enddo
2638  endif
2639 
2640  integral=czero
2641  if(ldiag) then
2642 
2643   if(option==1) then ! nspinor==1
2644     do ifreq=1,paw_dmft%dmft_nwlo
2645       if(procb2(ifreq)) then
2646 !     write(500,*) paw_dmft%omega_lo(ifreq),real(ff(ifreq)),imag(ff(ifreq))
2647        integral=integral+2.d0*paw_dmft%temp *                         &
2648 &        real( ff(ifreq)-one / ( j_dpc*paw_dmft%omega_lo(ifreq) ) ) *  &
2649 &        paw_dmft%wgt_wlo(ifreq)
2650       endif
2651     enddo
2652     if(procb2(paw_dmft%dmft_nwlo)) integral=integral+half
2653 !    integral=integral+half
2654     ! the if is here, to count only one time this correction
2655   endif
2656 
2657   if(option==2) then ! nspinor==2
2658     do ifreq=1,paw_dmft%dmft_nwlo
2659       if(procb2(ifreq)) then
2660         integral=integral+2.d0*paw_dmft%temp *                         &
2661 &           ( ff(ifreq)-one / ( j_dpc*paw_dmft%omega_lo(ifreq) ) ) *  &
2662 &       paw_dmft%wgt_wlo(ifreq)
2663       endif
2664     enddo
2665     if(procb2(paw_dmft%dmft_nwlo)) integral=integral+half
2666 !    integral=integral+half
2667     ! the if is here, to count only one time this correction
2668   endif
2669 
2670 
2671  else   ! ldiag
2672 
2673 ! write(std_out,*) "nondiag"
2674   if(option==1) then
2675     do ifreq=1,paw_dmft%dmft_nwlo
2676       if(procb2(ifreq)) then
2677         integral=integral+2.d0*paw_dmft%temp *   &
2678 &        real( ff(ifreq) ) *                     &
2679 &       paw_dmft%wgt_wlo(ifreq)
2680       endif
2681     enddo
2682   endif
2683   if(option==2) then
2684     do ifreq=1,paw_dmft%dmft_nwlo
2685       if(procb2(ifreq)) then
2686         integral=integral+2.d0*paw_dmft%temp *   &
2687 &            ff(ifreq)   *                     &
2688 &        paw_dmft%wgt_wlo(ifreq)
2689       endif
2690     enddo
2691   endif
2692  endif  ! ldiag
2693  ABI_DEALLOCATE(procb2)
2694 
2695 end subroutine int_fct

m_green/integrate_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 integrate_green

FUNCTION

  integrate green function in the band index basis

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  green <type(green_type)>=green function  (green%oper(:))
  opt_ksloc=   1: do only integrate on the KS basis
               2: do the integration on the local basis
                 (This can work only if psichi are renormalized!!!)
               3: do both calculations and test the consistency of it
              -1: do the integration on the KS basis, but only
                      compute diagonal part of the band-band density matrix
                      in order to compute the total charge for fermi_green
  paw_dmft <type(m_paw_dmft)>= paw+dmft data
  pawang <type(pawang)>=paw angular mesh and related data
  opt_nonxsum 0 : green(ifreq)%ks was broadcast in compute_green.
              1 : green(ifreq)%ks was not broadcast in compute_green, so a xsum will be
                  necessary here to compute the total number of electron.

OUTPUT

   green%occup = occupations

PARENTS

      m_dmft,fermi_green,impurity_solve,m_green,newton

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

1423 subroutine integrate_green(cryst_struc,green,paw_dmft&
1424 &  ,pawang,prtopt,opt_ksloc,opt_after_solver,opt_diff,opt_nonxsum)
1425 
1426  use m_pawang, only : pawang_type
1427  use m_crystal,   only : crystal_t
1428  use m_matlu,     only : sym_matlu,print_matlu,init_matlu,&
1429 & destroy_matlu,diff_matlu,zero_matlu
1430  use m_paw_dmft,  only : paw_dmft_type
1431  use m_oper,      only : loc_oper,trace_oper,init_oper,destroy_oper
1432 ! use m_oper, only : loc_oper,trace_oper,upfold_oper,print_oper,identity_oper,init_oper,destroy_oper
1433  use m_errors
1434 
1435 !This section has been created automatically by the script Abilint (TD).
1436 !Do not modify the following lines by hand.
1437 #undef ABI_FUNC
1438 #define ABI_FUNC 'integrate_green'
1439 !End of the abilint section
1440 
1441  implicit none
1442 
1443 !Arguments ------------------------------------
1444 !type
1445  type(crystal_t),intent(in) :: cryst_struc
1446  type(green_type),intent(inout) :: green
1447  !type(MPI_type), intent(in) :: mpi_enreg
1448  type(pawang_type),intent(in) :: pawang
1449  type(paw_dmft_type), intent(inout) :: paw_dmft
1450  integer, intent(in) :: prtopt
1451  integer, optional, intent(in) :: opt_ksloc
1452  integer, optional, intent(in) :: opt_after_solver
1453  integer, optional, intent(in) :: opt_diff
1454  integer, optional, intent(in) :: opt_nonxsum
1455 
1456 !local variables-------------------------------
1457  real(dp) :: tsec(2)
1458  integer :: iatom,ib,ib1,icomp_chloc,ifreq,ikpt,im,im1,ispinor,ispinor1,is
1459  integer :: mband,mbandc,myproc,natom,ndim,nkpt,nproc,nspinor
1460  integer :: nsppol,option
1461  integer :: optksloc,spacecomm,optaftsolv,optnonxsum
1462  complex(dpc) :: integral
1463  character(len=500) :: message
1464  complex(dpc), allocatable :: ff(:)
1465  type(matlu_type), allocatable :: matlu_temp(:)
1466  type(oper_type) :: occup_temp
1467  real(dp) :: diff_chloc
1468 ! real(dp), allocatable :: charge_loc_old(:,:)
1469 ! type(oper_type)  :: oper_c
1470 ! *********************************************************************
1471 
1472  DBG_ENTER("COLL")
1473  call timab(625,1,tsec)
1474 
1475  if(prtopt>0) then
1476    write(message,'(2a,i3,13x,a)') ch10,'   ===  Integrate green function'
1477    call wrtout(std_out,message,'COLL')
1478  endif
1479  if(green%w_type=="real") then
1480    message = 'integrate_green not implemented for real frequency'
1481    MSG_BUG(message)
1482  endif
1483  if(present(opt_nonxsum)) then
1484    optnonxsum=opt_nonxsum
1485  else
1486    optnonxsum=0
1487  endif
1488 
1489 ! Initialise spaceComm, myproc, and master
1490  spacecomm=paw_dmft%spacecomm
1491  myproc=paw_dmft%myproc
1492  nproc=paw_dmft%nproc
1493 
1494 
1495 ! Initialise integers
1496  mband   = paw_dmft%mband
1497  mbandc  = paw_dmft%mbandc
1498  nkpt    = paw_dmft%nkpt
1499  nspinor = paw_dmft%nspinor
1500  nsppol  = paw_dmft%nsppol
1501  natom   = paw_dmft%natom
1502 
1503 ! Initialize green%oper before calculation (important for xmpi_sum)
1504 ! allocate(charge_loc_old(paw_dmft%natom,paw_dmft%nsppol+1))
1505 ! if(.not.present(opt_diff)) then  ! if integrate_green is called in m_dmft after calculation of self
1506 !   charge_loc_old=green%charge_matlu
1507 ! endif
1508  icomp_chloc=0
1509 
1510 ! Choose what to compute
1511  if(present(opt_ksloc)) then
1512    optksloc=opt_ksloc
1513  else
1514    optksloc=3
1515  endif
1516  if(present(opt_after_solver)) then
1517    optaftsolv=opt_after_solver
1518  else
1519    optaftsolv=0
1520  endif
1521  if(optaftsolv==1.and.abs(optksloc)/=2) then
1522     message = "integration of ks green function should not be done after call to solver : it has not been computed"
1523     MSG_BUG(message)
1524  endif
1525  if(abs(optksloc)>=2.and.green%has_greenmatlu_xsum==0) then
1526     write(message,'(4a)') ch10,&
1527 &     "BUG: integrate_green is asked to integrate local green function",ch10,&
1528 &    " and local green function was non broadcasted in compute_green"
1529     MSG_BUG(message)
1530  endif
1531 
1532 ! Allocations
1533  ABI_DATATYPE_ALLOCATE(matlu_temp,(natom))
1534  call init_matlu(natom,nspinor,nsppol,green%oper(1)%matlu(:)%lpawu,matlu_temp)
1535 
1536  ABI_ALLOCATE(ff,(green%nw))
1537 
1538 ! =================================================
1539 ! == Integrate Local Green function ===============
1540  if(abs(optksloc)/2==1) then ! optksloc=2 or 3
1541 ! =================================================
1542    call zero_matlu(green%occup%matlu,green%occup%natom)
1543 
1544 ! ==  Calculation of \int{G_{LL'}{\sigma\sigma',s}(R)(i\omega_n)}
1545    if(paw_dmft%lpsichiortho==1) then
1546 !  - Calculation of frequency sum over positive frequency
1547      if (nspinor==1) option=1
1548      if (nspinor==2) option=2
1549      do iatom=1, natom
1550        ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
1551        if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
1552              do ispinor1 = 1, nspinor
1553            do ispinor = 1, nspinor
1554          do is = 1 , nsppol
1555                do im=1,ndim
1556                  do im1=1,ndim
1557                    do ifreq=1,green%nw
1558                      ff(ifreq)= &
1559 &                     green%oper(ifreq)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
1560 !                     write(std_out,*) green%omega(ifreq),ff(ifreq)," integrate green fw_lo"
1561    !if(present(iii).and.im==1.and.im1==1) write(std_out_default,*) ch10,ifreq,ff(ifreq),"#ff"
1562                    enddo
1563 !                   call int_fct(ff,(im==im1).and.(ispinor==ispinor1),&
1564 !&                   option,paw_dmft,integral)
1565                    call int_fct(ff,(im==im1).and.(ispinor==ispinor1),&
1566 &                   2,paw_dmft,integral)  ! test_1
1567                    green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=integral
1568    !if(present(iii).and.im==1.and.im1==1) write(std_out_default,*) ch10,'integral',im,im1,ifreq,integral
1569 !                   if(im==2.and.im1==5.and.is==1.and.iatom==1) then
1570 !                     write(std_out,*) " occup        %matlu(1)%mat(2,5,1,1,1)",green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)
1571 !                   endif
1572                  enddo
1573                enddo
1574              enddo ! ispinor1
1575            enddo ! ispinor
1576          enddo ! is
1577          matlu_temp(iatom)%mat=green%occup%matlu(iatom)%mat
1578        endif ! lpawu=/-1
1579      enddo ! iatom
1580 
1581 !   Print density matrix if prtopt high
1582      if(abs(prtopt)>2) then
1583        write(message,'(a,a,i10,a)') ch10,"  = green%occup%matlu from int(gloc(w))"
1584        call wrtout(std_out,message,'COLL')
1585        call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1)
1586      endif
1587 
1588 !  - Symetrise: continue sum over k-point: Full BZ
1589      call sym_matlu(cryst_struc,green%occup%matlu,pawang)
1590      if(abs(prtopt)>2) then
1591        write(message,'(a,a,i10,a)') ch10,&
1592 &       "  = green%occup%matlu from int(gloc(w)) with symetrization"
1593        call wrtout(std_out,message,'COLL')
1594        call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1)
1595      endif
1596      call sym_matlu(cryst_struc,matlu_temp,pawang)
1597 
1598 !  - Post-treatment for summation over negative and positive frequencies:
1599 !    necessary in the case of nspinor==2 AND nspinor==1, but valid anywhere
1600 !    N(ll'sigmasigma')= (N(ll'sigmasigma')+ N*(l'lsigma'sigma))/2
1601 !    because [G_{LL'}^{sigma,sigma'}(iomega_n)]*= G_{L'L}^{sigma',sigma}(-iomega_n)
1602      if(nspinor>=1) Then
1603        do iatom=1, natom
1604          ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
1605          if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
1606                do ispinor1 = 1, nspinor
1607              do ispinor = 1, nspinor
1608            do is = 1 , nsppol
1609                  do im=1,ndim
1610                    do im1=1,ndim
1611                      green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)= &
1612 &                          (matlu_temp(iatom)%mat(im,im1,is,ispinor,ispinor1)+ &
1613 &                           conjg(matlu_temp(iatom)%mat(im1,im,is,ispinor1,ispinor)))/two
1614                    enddo
1615                  enddo
1616                enddo ! ispinor1
1617              enddo ! ispinor
1618            enddo ! isppol
1619            matlu_temp(iatom)%mat=green%occup%matlu(iatom)%mat
1620          endif
1621        enddo ! iatom
1622      endif
1623      if(abs(prtopt)>2) then
1624        write(message,'(a,a,i10,a)') ch10,&
1625 &       "  = green%occup%matlu from int(gloc(w)) symetrized with post-treatment"
1626        call wrtout(std_out,message,'COLL')
1627        call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1)
1628      endif
1629      if(optaftsolv==0) then
1630        call trace_oper(green%occup,green%charge_ks,green%charge_matlu,2)
1631        green%has_charge_matlu=2
1632        green%has_charge_matlu_solver=0
1633        icomp_chloc=1
1634      else if(optaftsolv==1.and.paw_dmft%dmft_solv/=4) then
1635 !     else if(optaftsolv==1) then
1636        if(paw_dmft%dmft_solv==4) then
1637          write(message,'(a,a,a)') ch10,&
1638 &        "  = WARNING: Double Counting will be computed with solver charge.." ,&
1639 &        "might be problematic with Hirsch Fye QMC"
1640          call wrtout(std_out,message,'COLL')
1641        endif
1642 !      This is done only when called from impurity_solver with solver
1643 !      green function. (And if QMC is NOT used).
1644        if(paw_dmft%dmft_solv>=4.and.green%has_charge_matlu_solver/=2) then
1645          write(message,'(a,a,i3)') ch10,&
1646 &        "  = BUG : has_charge_matlu_solver should be 2 and is",&
1647 &        green%has_charge_matlu_solver
1648          MSG_BUG(message)
1649        endif
1650        if(paw_dmft%dmft_solv<=4) then
1651          call trace_oper(green%occup,green%charge_ks,green%charge_matlu_solver,2)
1652          green%has_charge_matlu_solver=2
1653        endif
1654      endif
1655    else
1656      write(message,'(a,4x,a,a,a,4x,a)') ch10,&
1657 &     " Local basis is not (yet) orthonormal:",&
1658 &     " local green function is thus not integrated",ch10,&
1659 &     " Local occupations are computed from KS occupations"
1660      call wrtout(std_out,message,'COLL')
1661    endif
1662 
1663  endif ! optksloc
1664 ! =================================================
1665 
1666 
1667 
1668 ! =================================================
1669 ! == Integrate Kohn Sham Green function ===========
1670  if(mod(abs(optksloc),2)==1) then ! optksloc=1 or 3 or -1
1671 !   green%occup%ks=czero ! important for xmpi_sum
1672 !! =================================================
1673 !! ==  Calculation of \int{G_{\nu\nu'}{k,s}(i\omega_n)}
1674    call init_oper(paw_dmft,occup_temp)
1675 !   ff=czero
1676 !   do is = 1 , nsppol
1677 !     do ikpt = 1, nkpt
1678 !!020212       if(mod(ikpt-1,nproc)==myproc) then
1679 !!!         print'(A,3I4)', "P ",myproc,is,ikpt
1680 !         do ib = 1, mbandc
1681 !           do ib1 = 1, mbandc
1682 !             if(optksloc==-1.and.(ib/=ib1)) cycle
1683 !             ff(:)=czero
1684 !             do ifreq=1,green%nw
1685 !             !!  the following line is a quick and dirty tricks and should be removed when the data will
1686 !             !!  be correctly distributed.
1687 !!!               if((optnonxsum==1).and.(green%procb(ifreq,ikpt)==myproc)).or.(optnonxsum==0)) then
1688 !               if(green%procb(ifreq,ikpt)==myproc) then
1689 !                 ff(ifreq)=green%oper(ifreq)%ks(is,ikpt,ib,ib1)
1690 !!!                 print'(A,5I4,3(2E15.5,3x))', "P1",myproc,is,ikpt,ib,ib1,ff(:)
1691 !               endif
1692 !               !endif
1693 !             enddo
1694 !!             call int_fct(ff,ib==ib1,nspinor,paw_dmft,integral) ! here, option==1 even if nspinor==2
1695 !!             green%occup%ks(is,ikpt,ib,ib1)=integral
1696 !!!                 print'(A,5I4,3(2E15.5,3x))', "PP",myproc,is,ikpt,ib,ib1,ff(:)
1697 !             call int_fct(ff,ib==ib1,2,paw_dmft,integral,green%procb(:,ikpt),myproc) ! here, option==1 even if nspinor==2
1698 !!020212             call int_fct(ff,ib==ib1,2,paw_dmft,integral) ! here, option==1 even if nspinor==2
1699 !             green%occup%ks(is,ikpt,ib,ib1)=integral
1700 !!!                 print'(A,5I4,2E15.8)', "P2",myproc,is,ikpt,ib,ib1,integral
1701 !            !   write(std_out,*) "integral",ikpt,green%occup%ks(is,ikpt,ib,ib1)
1702 !            ! endif
1703 !!        write(std_out,'(a,4i6,e14.5,e14.5)') "ks",is,ikpt,ib,ib1,integral
1704 !           enddo ! ib1
1705 !         enddo ! ib
1706 !!020212       endif
1707 !     enddo ! ikpt
1708 !   enddo ! isppol
1709 
1710 
1711 !   call xmpi_barrier(spacecomm)
1712 !   call xmpi_sum(green%occup%ks,spacecomm,ierr)
1713 
1714 
1715 !!   do is = 1 , nsppol
1716 !!     do ikpt = 1, nkpt
1717 !!         do ib = 1, mbandc
1718 !!           do ib1 = 1, mbandc
1719 !!             write(6,'(A,5I4,2E15.8)') "AAAFTERXSUM",myproc,is,ikpt,ib,ib1,green%occup%ks(is,ikpt,ib,ib1)
1720 !!             print '(A,5I4,2E15.8)', "AFTERXSUM P",myproc,is,ikpt,ib,ib1,green%occup%ks(is,ikpt,ib,ib1)
1721 !!           enddo ! ib1
1722 !!         enddo ! ib
1723 !!     enddo ! ikpt
1724 !!   enddo ! isppol
1725    occup_temp%ks=green%occup%ks
1726                !write(std_out,*) "occup%ks ik1",green%occup%ks(1,1,1,3)
1727                !write(std_out,*) "occup%ks ik2",green%occup%ks(1,2,1,3)
1728 !  - Post-treatment for summation over negative and positive frequencies:
1729 !    necessary in the case of nspinor==2, but valid everywhere
1730 !    N(k,n_1,n_2)= (N(k,n_1,n_2)+ N*(k,n_2,n_1))/2
1731 !    because [G_{k}^{n_1,n_2}(iomega_n)]*= G_{k}^{n_2,n_1}(-iomega_n)
1732          do ib1 = 1, mbandc
1733        do ib = 1, mbandc
1734      do ikpt = 1, nkpt
1735    do is = 1 , nsppol
1736            green%occup%ks(is,ikpt,ib,ib1)=&
1737 &            (       occup_temp%ks(is,ikpt,ib,ib1)+ &
1738 &              conjg(occup_temp%ks(is,ikpt,ib1,ib))   )/two
1739          enddo ! ib1
1740        enddo ! ib
1741      enddo ! ikpt
1742    enddo ! isppol
1743                !write(std_out,*) "occup%ks ik1 BB",green%occup%ks(1,1,1,3)
1744                !write(std_out,*) "occup%ks ik2 AA",green%occup%ks(1,2,1,3)
1745    call destroy_oper(occup_temp)
1746          do ib1 = 1, mbandc
1747        do ib = 1, mbandc
1748      do ikpt = 1, nkpt
1749    do is = 1 , nsppol
1750            paw_dmft%occnd(1,paw_dmft%include_bands(ib),&
1751 &           paw_dmft%include_bands(ib1),ikpt,is)=real(green%occup%ks(is,ikpt,ib,ib1))
1752            if(nspinor==1) then
1753              paw_dmft%occnd(2,paw_dmft%include_bands(ib),&
1754 &             paw_dmft%include_bands(ib1),ikpt,is)=zero
1755              if(nsppol==1) then
1756                paw_dmft%occnd(1,paw_dmft%include_bands(ib),&
1757 &               paw_dmft%include_bands(ib1),ikpt,is)=two*real(green%occup%ks(is,ikpt,ib,ib1))
1758              endif
1759            else if (nspinor==2) then  ! and SOC
1760              paw_dmft%occnd(2,paw_dmft%include_bands(ib),&
1761 &             paw_dmft%include_bands(ib1),ikpt,is)=aimag(green%occup%ks(is,ikpt,ib,ib1))
1762            endif
1763          enddo
1764        enddo
1765      enddo
1766    enddo
1767 
1768    if(optksloc>0) then
1769 !  - Compute local occupations
1770      call loc_oper(green%occup,paw_dmft,1)
1771      if(abs(prtopt)>2) then
1772        write(message,'(a,a,i10,a)') ch10,&
1773 &        "  = green%occup%matlu from projection of int(gks(w)) without symetrization"
1774        call wrtout(std_out,message,'COLL')
1775        call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1)
1776      endif
1777 
1778 !  - Symetrise: continue sum over k-point: Full BZ
1779      call sym_matlu(cryst_struc,green%occup%matlu,pawang)
1780      if(abs(prtopt)>=2) then
1781 !       write(message,'(a,a,i10,a)') ch10,&
1782 !&        "  = green%occup%matlu from projection of int(gks(w)) with symetrization"
1783        write(message,'(a,a,i10,a)') ch10,&
1784 &        "  = Occupation matrix from KS occupations"
1785        call wrtout(std_out,message,'COLL')
1786        call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=0)
1787      endif
1788 
1789 !  - If the trace of occup matrix in the LOCAL basis was not done
1790 !  before because lpsichiortho/=1 , do it now
1791      if(paw_dmft%lpsichiortho/=1) then
1792        if(abs(prtopt)>0) then
1793          call trace_oper(green%occup,green%charge_ks,green%charge_matlu,2)
1794          green%has_charge_matlu=2
1795          icomp_chloc=1
1796        endif
1797      endif
1798    endif ! optksloc>0: only diagonal elements of G\nunu' are computed
1799 
1800 !  - Compute trace over ks density matrix
1801    call trace_oper(green%occup,green%charge_ks,green%charge_matlu,1)
1802    if(abs(prtopt)>0) then
1803      write(message,'(a,a,f12.6)') ch10,&
1804 &    "  ==  Total number of electrons from KS green function is :", green%charge_ks
1805      call wrtout(std_out,message,'COLL')
1806      write(message,'(8x,a,f12.6,a)') " (should be",paw_dmft%nelectval,")"
1807      call wrtout(std_out,message,'COLL')
1808    endif
1809  endif ! optksloc
1810 ! =================================================
1811 
1812 
1813 ! =================================================
1814 ! Tests and compute precision on local charge
1815 ! =================================================
1816 !  - Check that if, renormalized psichi are used, occupations matrices
1817 !    obtained directly from local green function or, through kohn sham
1818 !    occupations are the same.
1819  if(abs(optksloc)==3) then ! optksloc= 3
1820    if(paw_dmft%lpsichiortho==1) then
1821      call diff_matlu("Local_projection_of_kohnsham_occupations ",&
1822 &     "Integration_of_local_green_function ",&
1823 &       green%occup%matlu,matlu_temp,natom,1,tol4)
1824      write(message,'(2a)') ch10,&
1825 &     '  ***** => Calculations of Green function in KS and local spaces are coherent ****'
1826      call wrtout(std_out,message,'COLL')
1827    endif
1828  endif

m_green/local_ks_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 local_ks_green

FUNCTION

 Compute the sum over k-point of ks green function.
 do the fourier transformation and print it

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

  cryst_struc
  istep    =  step of iteration for LDA.
  lda_occup
  mpi_enreg=informations about MPI parallelization
  paw_dmft =  data for self-consistent LDA+DMFT calculations.

OUTPUT

  paw_dmft =  data for self-consistent LDA+DMFT calculations.

PARENTS

      m_dmft

CHILDREN

      fourier_fct,wrtout

SOURCE

3893 subroutine local_ks_green(green,paw_dmft,prtopt)
3894 
3895  use defs_basis
3896  use m_errors
3897  use m_abicore
3898 
3899  use m_crystal, only : crystal_t
3900  use m_paw_dmft, only : paw_dmft_type
3901 
3902 !This section has been created automatically by the script Abilint (TD).
3903 !Do not modify the following lines by hand.
3904 #undef ABI_FUNC
3905 #define ABI_FUNC 'local_ks_green'
3906 !End of the abilint section
3907 
3908  implicit none
3909 
3910 !Arguments ------------------------------------
3911 !scalars
3912  type(green_type), intent(in) :: green
3913  type(paw_dmft_type), intent(in)  :: paw_dmft
3914  integer, intent(in) :: prtopt
3915 
3916 !Local variables ------------------------------
3917  character(len=500) :: message
3918  integer :: iband,ifreq,ikpt,isppol,itau,lsub,ltau,mbandc,nkpt,nsppol
3919  character(len=1) :: tag_is
3920  character(len=fnlen) :: tmpfil
3921  integer,allocatable :: unitgreenlocks_arr(:)
3922  real(dp) :: beta
3923  real(dp), allocatable :: tau(:)
3924  complex(dpc), allocatable :: loc_ks(:,:,:)
3925  complex(dpc), allocatable :: loc_ks_tau(:,:,:),fw(:),ft(:)
3926 !scalars
3927 !************************************************************************
3928  mbandc=paw_dmft%mbandc
3929  nkpt=paw_dmft%nkpt
3930  nsppol=paw_dmft%nsppol
3931  ltau=128
3932  ABI_ALLOCATE(tau,(ltau))
3933  do itau=1,ltau
3934    tau(itau)=float(itau-1)/float(ltau)/paw_dmft%temp
3935  end do
3936  beta=one/paw_dmft%temp
3937 
3938 !Only imaginary frequencies here
3939  if(green%w_type=="real") then
3940    message = ' compute_energy not implemented for real frequency'
3941    MSG_BUG(message)
3942  end if
3943 
3944 !=========================================
3945 !Compute local band ks green function
3946 ! should be computed in compute_green: it would be less costly in memory.
3947 !=========================================
3948  ABI_ALLOCATE(loc_ks,(nsppol,mbandc,paw_dmft%dmft_nwlo))
3949  if(green%oper(1)%has_operks==1) then
3950    loc_ks(:,:,:)=czero
3951    do isppol=1,nsppol
3952      do iband=1,mbandc
3953        do ifreq=1,paw_dmft%dmft_nwlo
3954          do ikpt=1,nkpt
3955            loc_ks(isppol,iband,ifreq)=loc_ks(isppol,iband,ifreq)+  &
3956 &           green%oper(ifreq)%ks(isppol,ikpt,iband,iband)*paw_dmft%wtk(ikpt)
3957          end do
3958        end do
3959      end do
3960    end do
3961  else
3962    message = ' green fct is not computed in ks space'
3963    MSG_BUG(message)
3964  end if
3965 
3966 !=========================================
3967 !Compute fourier transformation
3968 !=========================================
3969 
3970  ABI_ALLOCATE(loc_ks_tau,(nsppol,mbandc,ltau))
3971  ABI_ALLOCATE(fw,(paw_dmft%dmft_nwlo))
3972  ABI_ALLOCATE(ft,(ltau))
3973  loc_ks_tau(:,:,:)=czero
3974  do isppol=1,nsppol
3975    do iband=1,mbandc
3976      do ifreq=1,paw_dmft%dmft_nwlo
3977        fw(ifreq)=loc_ks(isppol,iband,ifreq)
3978      end do
3979      call fourier_fct(fw,ft,.true.,ltau,-1,paw_dmft) ! inverse fourier
3980      do itau=1,ltau
3981        loc_ks_tau(isppol,iband,itau)=ft(itau)
3982      end do
3983    end do
3984  end do
3985  ABI_DEALLOCATE(fw)
3986  ABI_DEALLOCATE(ft)
3987  do isppol=1,nsppol
3988    do iband=1,mbandc
3989      do itau=1,ltau
3990        loc_ks_tau(isppol,iband,itau)=(loc_ks_tau(isppol,iband,itau)+conjg(loc_ks_tau(isppol,iband,itau)))/two
3991      end do
3992    end do
3993  end do
3994 
3995 !=========================================
3996 !Print out ksloc green function
3997 !=========================================
3998  if(abs(prtopt)==1) then
3999    ABI_ALLOCATE(unitgreenlocks_arr,(nsppol))
4000    do isppol=1,nsppol
4001      write(tag_is,'(i1)')isppol
4002      tmpfil = trim(paw_dmft%filapp)//'Gtau_locks_isppol'//tag_is
4003      write(message,'(3a)') ch10," == Print green function on file ",tmpfil
4004      call wrtout(std_out,message,'COLL')
4005      unitgreenlocks_arr(isppol)=500+isppol-1
4006      open (unit=unitgreenlocks_arr(isppol),file=trim(tmpfil),status='unknown',form='formatted')
4007      rewind(unitgreenlocks_arr(isppol))
4008      write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenlocks_arr(isppol)
4009      write(message,'(a,a)') ch10,"# New record : First 40 bands"
4010      call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4011      do lsub=1,mbandc/40+1
4012        do itau=1, ltau
4013          write(message,'(2x,50(e10.3,2x))') tau(itau), &
4014 &         (real(loc_ks_tau(isppol,iband,itau)),iband=40*(lsub-1)+1,min(40*lsub,mbandc))
4015          call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4016        end do
4017        write(message,'(2x,50(e10.3,2x))') beta, &
4018 &       ((-one-real(loc_ks_tau(isppol,iband,1))),iband=40*(lsub-1)+1,min(40*lsub,mbandc))
4019        call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4020        if(40*lsub<mbandc) then
4021          write(message,'(a,a,i5,a,i5)')    &
4022 &         ch10,"# Same record, Following bands : From ",    &
4023 &         40*(lsub),"  to ",min(40*(lsub+1),mbandc)
4024          call wrtout(unitgreenlocks_arr(isppol),message,'COLL')
4025        end if
4026      end do
4027 !    call flush(unitgreenlocks_arr(isppol))
4028    end do
4029    ABI_DEALLOCATE(unitgreenlocks_arr)
4030  end if
4031 
4032 !Deallocations
4033  ABI_DEALLOCATE(loc_ks)
4034  ABI_DEALLOCATE(loc_ks_tau)
4035  ABI_DEALLOCATE(tau)
4036 
4037 end subroutine local_ks_green

m_green/newton [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 newton

FUNCTION

  Compute root of a function with newton methods (newton/halley)

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

  x_input      : input of x
  x_precision  : required precision on x
  max_iter     : maximum number of iterations
  opt_noninter

OUTPUT

  f_precision  : output precision on function F
  ierr_hh      : different from zero if an error occurs

PARENTS

      fermi_green

CHILDREN

      compute_green,integrate_green

SOURCE

3479 subroutine newton(cryst_struc,green,paw_dmft,pawang,self,&
3480 & x_input,x_precision,max_iter,f_precision,ierr_hh,opt_noninter,opt_algo)
3481 
3482  use defs_basis
3483  use defs_abitypes
3484  use m_abicore
3485  use m_errors
3486 
3487  use m_pawang,    only : pawang_type
3488  use m_crystal,   only : crystal_t
3489  use m_paw_dmft,  only: paw_dmft_type
3490  use m_self,      only : self_type
3491 
3492 !This section has been created automatically by the script Abilint (TD).
3493 !Do not modify the following lines by hand.
3494 #undef ABI_FUNC
3495 #define ABI_FUNC 'newton'
3496 !End of the abilint section
3497 
3498  implicit none
3499 !Arguments ------------------------------------
3500 !scalars
3501  type(crystal_t),intent(in) :: cryst_struc
3502  type(green_type),intent(inout) :: green
3503  !type(MPI_type), intent(in) :: mpi_enreg
3504  type(paw_dmft_type), intent(inout) :: paw_dmft
3505  type(pawang_type),intent(in) :: pawang
3506  type(self_type), intent(inout) :: self
3507  integer,intent(in) :: opt_noninter,max_iter
3508  integer,intent(out) :: ierr_hh
3509  real(dp),intent(inout) :: x_input,x_precision
3510  real(dp),intent(inout) :: f_precision
3511  real(dp),intent(in), optional :: opt_algo
3512 !Local variables-------------------------------
3513  integer iter
3514  real(dp) Fx,Fxprime,Fxdouble,xold,x_before,Fxoptimum
3515  real(dp) :: nb_elec_x
3516  integer option,optalgo
3517  logical l_minus,l_plus
3518  real(dp) :: x_minus,x_plus,x_optimum
3519  character(len=500) :: message
3520 ! *********************************************************************
3521  x_minus=-10_dp
3522  x_plus=-11_dp
3523  xold=-12_dp
3524 
3525  if(present(opt_algo)) then
3526    optalgo=opt_algo
3527  else
3528    optalgo=1 ! newton
3529  end if
3530 
3531  x_input=paw_dmft%fermie
3532  ierr_hh=0
3533  option =2  ! Halley method
3534  option =1  ! Newton method
3535 
3536 !write(std_out,*) "ldaprint",opt_noninter
3537 
3538 !--- Start of iterations
3539  write(message,'(a,3a)') "     Fermi level   Charge    Difference"
3540  call wrtout(std_out,message,'COLL')
3541 !do iter=1, 40
3542 !x_input=float(iter)/100_dp
3543 !call function_and_deriv(cryst_struc,f_precision,green,iter,mpi_enreg,paw_dmft,pawang,self,&
3544 !&  x_input,x_before,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option)
3545 !write(std_out,*) x_input,Fx
3546 !enddo
3547 !call abi_abort('COLL')
3548 
3549  l_minus=.false.
3550  l_plus=.false.
3551  Fxoptimum=1_dp
3552 !========================================
3553 !start iteration to find fermi level
3554 !========================================
3555  do iter=1, max_iter
3556 
3557 !  ========================================
3558 !  If zero is located between two values: apply newton method or dichotomy
3559 !  ========================================
3560    if(l_minus.and.l_plus) then
3561 
3562 !    ==============================================
3563 !    Compute the function and derivatives for newton
3564 !    ==============================================
3565      call function_and_deriv(cryst_struc,f_precision,green,iter,paw_dmft,pawang,self,&
3566 &     x_input,x_before,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option)
3567 
3568 !    Apply stop criterion on Fx
3569      if(abs(Fx) < f_precision) then
3570 !      write(message,'(a,2f12.6)') "Fx,f_precision",Fx,f_precision
3571 !      call wrtout(std_out,message,'COLL')
3572        x_precision=x_input-xold
3573        return
3574      end if
3575      if(iter==max_iter) then
3576        write(message,'(a,2f12.6)') "   Fermi level could not be found"
3577        call wrtout(std_out,message,'COLL')
3578        x_input=x_optimum
3579        ierr_hh=-123
3580        return
3581      end if
3582 
3583 !    Cannot divide by Fxprime if too small
3584      if(abs(Fxprime) .le. 1.e-15)then
3585        ierr_hh=-314
3586        write(message,'(a,f12.7)') "Fxprime=",Fxprime
3587        call wrtout(std_out,message,'COLL')
3588        return
3589      end if
3590 
3591      x_precision=x_input-xold
3592 
3593 !    ==============================================
3594 !    Newton/Halley's  formula for next iteration
3595 !    ==============================================
3596      xold=x_input
3597      if(option==1) x_input=x_input-Fx/Fxprime
3598      if(option==2) x_input=x_input-2*Fx*Fxprime/(2*Fxprime**2-Fx*Fxdouble)
3599 
3600 !    ==============================================
3601 !    If newton does not work well, use dichotomy.
3602 !    ==============================================
3603      if(x_input<x_minus.or.x_input>x_plus) then
3604        call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
3605 &       Fx,opt_noninter,nb_elec_x,xold)
3606        write(message,'(a,3f12.6)') " ---",x_input,nb_elec_x,Fx
3607        call wrtout(std_out,message,'COLL')
3608        if(Fx>0) then
3609          x_plus=xold
3610        else if(Fx<0) then
3611          x_minus=xold
3612        end if
3613        x_input=(x_plus+x_minus)/2.d0
3614 
3615      end if
3616 !    write(std_out,'(a,2f12.6)') " Q(xold) and dQ/dx=",Fx,Fxprime
3617 !    write(std_out,'(a,f12.6)') " =>  new Fermi level",x_input
3618 !    ========================================
3619 !    Locate zero between two values
3620 !    ========================================
3621    else
3622      call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
3623 &     Fx,opt_noninter,nb_elec_x,x_input)
3624      write(message,'(a,3f12.6)') "  --",x_input,nb_elec_x,Fx
3625 !    Possible improvement for large systems, removed temporarely for
3626 !    automatic tests: more study is necessary: might worsen the convergency
3627 !    if(iter==1) then
3628 !    f_precision=max(abs(Fx/50),f_precision)
3629 !    write(message,'(a,4x,a,e12.6)') ch10," Precision required changed to:",f_precision
3630 !    call wrtout(std_out,message,'COLL')
3631 !    endif
3632      call wrtout(std_out,message,'COLL')
3633      if(Fx>0) then
3634        l_plus=.true.
3635        x_plus=x_input
3636        x_input=x_input-0.02
3637      else if(Fx<0) then
3638        l_minus=.true.
3639        x_minus=x_input
3640        x_input=x_input+0.02
3641      end if
3642 
3643    end if
3644 
3645    if(abs(Fx)<abs(Fxoptimum)) then
3646      Fxoptimum=Fx
3647      x_optimum=x_input
3648    end if
3649 
3650 
3651 
3652 !  if(myid==master) then
3653 !  write(std_out,'(a,i4,3f12.6)') "i,xnew,F,Fprime",i,x_input,Fx,Fxprime
3654 !  endif
3655 
3656 
3657 !  ! Apply stop criterion on x
3658 !  if(abs(x_input-xold) .le. x_input*x_precision) then
3659 !  !    write(std_out,'(a,4f12.6)') "x_input-xold, x_precision*x_input   "&
3660 !  !&    ,x_input-xold,x_precision*x_input,x_precision
3661 !  f_precision=Fx
3662 !  return
3663 !  endif
3664 
3665  end do
3666 !--- End of iterations
3667 
3668 
3669  ierr_hh=1
3670  return
3671 
3672  CONTAINS  !========================================================================================
3673 !-----------------------------------------------------------------------

m_green/occup_green_tau [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 occup_green_tau

FUNCTION

  Compute occup_tau from green%oper_tau

INPUTS

  green  <type(green_type)>= green function data

OUTPUT

PARENTS

      qmc_prep_ctqmc

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2963 subroutine occup_green_tau(green)
2964 
2965  use m_matlu, only : shift_matlu
2966 
2967 !This section has been created automatically by the script Abilint (TD).
2968 !Do not modify the following lines by hand.
2969 #undef ABI_FUNC
2970 #define ABI_FUNC 'occup_green_tau'
2971 !End of the abilint section
2972 
2973  implicit none
2974 
2975 !Arguments ------------------------------------
2976 !type
2977  type(green_type),intent(inout) :: green
2978 
2979 !local variables-------------------------------
2980  integer :: iatom,lpawu
2981  complex(dpc), allocatable :: shift(:)
2982 ! *********************************************************************
2983 
2984  ABI_ALLOCATE(shift,(green%oper_tau(1)%natom))
2985  do iatom=1,green%oper_tau(1)%natom
2986    shift(iatom)=-cone
2987    lpawu=green%oper_tau(1)%matlu(iatom)%lpawu
2988    if(lpawu/=-1) then
2989 !     do isppol=1,green%oper_tau(1)%nsppol
2990      green%occup_tau%matlu(iatom)%mat(:,:,:,:,:)= &
2991 &      green%oper_tau(1)%matlu(iatom)%mat(:,:,:,:,:)
2992    endif
2993  enddo
2994  call shift_matlu(green%occup_tau%matlu,green%occup_tau%natom,shift)
2995  ABI_DEALLOCATE(shift)
2996 
2997 
2998 end subroutine occup_green_tau

m_green/occupfd [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 occupfd

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

      wrtout

SOURCE

3019  function occupfd(eig,fermie,temp)
3020 
3021 
3022 !This section has been created automatically by the script Abilint (TD).
3023 !Do not modify the following lines by hand.
3024 #undef ABI_FUNC
3025 #define ABI_FUNC 'occupfd'
3026 !End of the abilint section
3027 
3028  implicit none
3029 
3030 !Arguments ------------------------------------
3031 !type
3032 ! Integrate analytic tail 1/(iw-mu)
3033  real(dp),intent(in) :: eig,fermie,temp
3034  real(dp) :: occupfd
3035 !Local variables-------------------------------
3036 ! *********************************************************************
3037 
3038  if((eig-fermie) > zero) then
3039    occupfd=exp(-(eig-fermie)/temp)/(one+exp(-(eig-fermie)/temp))
3040  else
3041    occupfd=one/(one+exp((eig-fermie)/temp))
3042  endif
3043 
3044  end function occupfd

m_green/print_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 print_green

FUNCTION

  print green function

INPUTS

  char1 = character which describes the type of green function
  green  <type(green_type)>= green function data
  option=1 print local green function
         2 print KS green function
         3 print both local and KS green function
         4 print spectral function is green%w_type="real"
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawprtvol = printing option
  opt_wt=1 print green function as a function of frequency
         2 print green function as a function of imaginary time

OUTPUT

PARENTS

      impurity_solve,m_green,qmc_prep_ctqmc,spectral_function

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

 740 subroutine print_green(char1,green,option,paw_dmft,pawprtvol,opt_wt,opt_decim)
 741 
 742  use m_crystal, only : crystal_t
 743  use m_oper, only : print_oper
 744  use m_paw_dmft, only : paw_dmft_type
 745 
 746  use m_fstrings, only : int2char4
 747 
 748 !This section has been created automatically by the script Abilint (TD).
 749 !Do not modify the following lines by hand.
 750 #undef ABI_FUNC
 751 #define ABI_FUNC 'print_green'
 752 !End of the abilint section
 753 
 754  implicit none
 755 
 756 !Arguments ------------------------------------
 757 !type
 758  type(paw_dmft_type), intent(in) :: paw_dmft
 759  type(green_type),intent(inout) :: green
 760  integer,intent(in) :: option,pawprtvol
 761  integer,intent(in),optional :: opt_wt,opt_decim
 762  character(len=*), intent(in) :: char1
 763 
 764 !local variables-------------------------------
 765  integer :: iall,iatom,ib,ifreq,ikpt,im,ispinor,isppol,itau
 766  integer :: lsub,mbandc,natom,ndim,nkpt,nspinor,nsppol,optwt,spf_unt,spcorb_unt
 767  character(len=2000) :: message
 768  integer,allocatable :: unitgreenfunc_arr(:)
 769  integer,allocatable :: unitgreenloc_arr(:)
 770  character(len=fnlen) :: tmpfil
 771  character(len=1) :: tag_is,tag_is2
 772  character(len=10) :: tag_at
 773  character(len=3) :: tag_ik
 774  real(dp) :: re, ima
 775  complex(dpc), allocatable :: sf(:),sf_corr(:)
 776 ! *********************************************************************
 777  if(present(opt_wt)) then
 778    optwt=opt_wt
 779  else
 780    optwt=1
 781  endif
 782 
 783 
 784  if(pawprtvol>200) then
 785  endif
 786  natom=green%oper(1)%natom
 787  nsppol=paw_dmft%nsppol
 788  nspinor=paw_dmft%nspinor
 789  mbandc=paw_dmft%mbandc
 790  nkpt=paw_dmft%nkpt
 791 
 792 !  == Print local Green Function
 793  if(option==1.or.option==3) then
 794    ABI_ALLOCATE(unitgreenfunc_arr,(natom*nsppol*nspinor))
 795    iall=0
 796    do iatom=1,natom
 797      if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
 798        ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
 799        call int2char4(iatom,tag_at)
 800        ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
 801        do isppol=1,nsppol
 802          write(tag_is,'(i1)')isppol
 803          do ispinor=1,nspinor
 804            iall=iall+1
 805            write(tag_is2,'(i1)')ispinor
 806 
 807 !         == Create names
 808            if(optwt==1) then
 809              tmpfil = &
 810 &             trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-omega_iatom'// &
 811 &             trim(tag_at)//'_isppol'//tag_is//'_ispinor'//tag_is2
 812            else
 813              tmpfil = &
 814 &             trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-tau_iatom'// &
 815 &             trim(tag_at)//'_isppol'//tag_is//'_ispinor'//tag_is2
 816            endif
 817            if(iall<=4) then
 818              write(message,'(3a)') ch10,"  == Print green function on file ",tmpfil
 819              call wrtout(std_out,message,'COLL')
 820            elseif(iall==5)  then
 821              write(message,'(3a)') ch10,"  == following values are printed in files"
 822              call wrtout(std_out,message,'COLL')
 823            endif
 824            unitgreenfunc_arr(iall)=300+iall-1
 825            if(optwt==1.or.green%fileprt_tau==0) then
 826              open(unit=unitgreenfunc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted',position='append')
 827            else if(optwt==2.and.green%fileprt_tau==1) then
 828              open(unit=unitgreenfunc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted',position='append')
 829            endif
 830 
 831 !         == Write in files
 832 !           rewind(unitgreenfunc_arr(iall))
 833 !           write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenfunc_arr(iall)
 834 !           call wrtout(std_out,message,'COLL')
 835            write(message,'(a,a)') ch10,"# New record :"
 836            call wrtout(unitgreenfunc_arr(iall),message,'COLL')
 837            if(optwt==1) then
 838              do ifreq=1,green%nw
 839                if(present(opt_decim)) then
 840                  write(message,'(2x,30(e23.16,2x))') &
 841 &                green%omega(ifreq),&
 842 &                (green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor),im=1,ndim)
 843                else
 844                  write(message,'(2x,30(e10.3,2x))') &
 845 &                green%omega(ifreq),&
 846 &                (green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor),im=1,ndim)
 847                endif
 848                call wrtout(unitgreenfunc_arr(iall),message,'COLL')
 849                re=real(green%oper(ifreq)%matlu(iatom)%mat(1,1,isppol,ispinor,ispinor))
 850                ima=aimag(green%oper(ifreq)%matlu(iatom)%mat(1,1,isppol,ispinor,ispinor))
 851 !               write(228,*) green%omega(ifreq),re/(re**2+ima**2),ima/(re**2+ima**2)+green%omega(ifreq)
 852              enddo
 853            else
 854              do itau=1,green%dmftqmc_l
 855                  write(message,'(2x,30(e10.3,2x))') &
 856 &                green%tau(itau),&
 857 &               (green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor),im=1,ndim)
 858                call wrtout(unitgreenfunc_arr(iall),message,'COLL')
 859              enddo
 860              write(message,'(2x,30(e10.3,2x))') &
 861 &            one/paw_dmft%temp,&
 862 &            (-green%oper_tau(1)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)-one,im=1,ndim)
 863              call wrtout(unitgreenfunc_arr(iall),message,'COLL')
 864            endif
 865          enddo
 866        enddo ! isppol
 867      endif ! lpawu=/-1
 868    enddo ! iatom
 869    ABI_DEALLOCATE(unitgreenfunc_arr)
 870  endif
 871 
 872 !  == Print ks green function
 873  if((option==2.or.option==3).and.green%oper(1)%has_operks==1) then
 874    ABI_ALLOCATE(unitgreenloc_arr,(nsppol*nkpt))
 875    iall=0
 876    do isppol = 1 , nsppol
 877      write(tag_is,'(i1)')isppol
 878      do ikpt = 1, nkpt
 879        write(tag_ik,'(i3)')ikpt
 880 !      do ib1 = 1, mbandc
 881        iall=iall+1
 882 
 883 !         == Create names
 884        if(optwt==1) then
 885          tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-omega_isppol'//tag_is//'_ikpt'//trim(adjustl(tag_ik))
 886        else
 887          tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-tau_isppol'//tag_is//'_ikpt'//trim(adjustl(tag_ik))
 888        endif
 889        if(iall<=4)  then
 890          write(message,'(3a)') ch10,"  == Print green function on file ",tmpfil
 891          call wrtout(std_out,message,'COLL')
 892        elseif(iall==5)  then
 893          write(message,'(3a)') ch10,"  == following values are printed in files"
 894          call wrtout(std_out,message,'COLL')
 895        endif
 896        unitgreenloc_arr(iall)=400+iall-1
 897        open (unit=unitgreenloc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted')
 898 !      rewind(unitgreenloc_arr(iall))
 899 !       write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenloc_arr(iall)
 900 !       call wrtout(std_out,message,'COLL')
 901 
 902 !         == Write in files
 903        write(message,'(a,a)') ch10,"# New record : First 20 bands"
 904        call wrtout(unitgreenloc_arr(iall),message,'COLL')
 905 !       call flush(std_out)
 906        do lsub=1,mbandc/20+1
 907          if(optwt==1) then
 908            do ifreq=1,green%nw
 909 !             call flush(std_out)
 910              write(message,'(2x,50(e10.3,2x))') &
 911 &            green%omega(ifreq), &
 912 &            (green%oper(ifreq)%ks(isppol,ikpt,ib,ib),ib=20*(lsub-1)+1,min(20*lsub,mbandc))
 913              call wrtout(unitgreenloc_arr(iall),message,'COLL')
 914                re=real(green%oper(ifreq)%ks(isppol,ikpt,1,1))
 915                ima=aimag(green%oper(ifreq)%ks(isppol,ikpt,1,1))
 916                if(ikpt==1) write(229,*) green%omega(ifreq),re/(re**2+ima**2),ima/(re**2+ima**2)+green%omega(ifreq)
 917            enddo
 918          else
 919            if(green%use_oper_tau_ks==1) then
 920              do itau=1,green%dmftqmc_l
 921 !               call flush(std_out)
 922                write(message,'(2x,50(e10.3,2x))') &
 923 &              green%tau(itau), &
 924 &              (green%oper_tau(itau)%ks(isppol,ikpt,ib,ib),ib=20*(lsub-1)+1,min(20*lsub,mbandc))
 925                call wrtout(unitgreenloc_arr(iall),message,'COLL')
 926              enddo
 927            endif
 928          endif
 929          if(20*lsub<mbandc) write(message,'(a,a,i5,a,i5)')    &
 930 &           ch10,"# Same record, Following bands : From ",    &
 931 &           20*(lsub),"  to ",min(20*(lsub+1),mbandc)
 932          call wrtout(unitgreenloc_arr(iall),message,'COLL')
 933        enddo
 934      enddo ! ikpt
 935    enddo ! isppol
 936    ABI_DEALLOCATE(unitgreenloc_arr)
 937  endif
 938 
 939  if((green%w_type=="real".and.option==4).and.green%oper(1)%has_operks==1) then
 940    write(message,'(a,a)') ch10,"  == About to print spectral function"
 941    call wrtout(std_out,message,'COLL')
 942    tmpfil = trim(paw_dmft%filapp)//'SpFunc-'//trim(char1)
 943    if (open_file(tmpfil, message, newunit=spf_unt, status='unknown', form='formatted') /= 0) then
 944      MSG_ERROR(message)
 945    end if
 946    ABI_ALLOCATE(sf,(green%nw))
 947    ABI_ALLOCATE(sf_corr,(green%nw))
 948    iall=0
 949    sf=czero
 950    do isppol = 1 , nsppol
 951      do ikpt = 1, nkpt
 952        do ib=1,mbandc
 953          iall=iall+1
 954          do ifreq=1,green%nw
 955            sf(ifreq)=sf(ifreq)+green%oper(ifreq)%ks(isppol,ikpt,ib,ib)*green%oper(1)%wtk(ikpt)
 956          enddo
 957        enddo
 958      enddo
 959    enddo
 960    sf_corr=czero
 961    do iatom=1,natom
 962      call int2char4(iatom,tag_at)
 963      ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
 964      tmpfil = trim(paw_dmft%filapp)//'SpFunc-cor_orb-'//trim(char1)//'_iatom'//trim(tag_at)
 965      if (open_file(tmpfil, message, newunit=spcorb_unt, status='unknown', form='formatted') /= 0) then
 966        MSG_ERROR(message)
 967      end if
 968      write(message,*) "#", nspinor,nsppol,ndim,green%nw
 969      call wrtout(spcorb_unt,message,'COLL')
 970      if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then
 971        write(message,*) "#", green%oper(1)%matlu(iatom)%lpawu
 972        call wrtout(spcorb_unt,message,'COLL')
 973        ndim=2*green%oper(1)%matlu(iatom)%lpawu+1
 974        do isppol = 1 , nsppol
 975          do ispinor=1, nspinor
 976            do im=1,ndim
 977              do ifreq=1,green%nw
 978                sf_corr(ifreq)=sf_corr(ifreq)+ green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)
 979              enddo
 980            enddo
 981          enddo
 982        enddo
 983      endif
 984      do ifreq=1,green%nw
 985        write(message,*) green%omega(ifreq),(-aimag(sf_corr(ifreq)))/3.141592653589793238_dp
 986        call wrtout(spcorb_unt,message,'COLL')
 987      enddo
 988      close(spcorb_unt)
 989    enddo
 990    do ifreq=1,green%nw
 991      write(message,*) green%omega(ifreq),(-aimag(sf(ifreq)))/3.141592653589793238_dp
 992      call wrtout(spf_unt,message,'COLL')
 993    enddo
 994    close(spf_unt)
 995    ABI_DEALLOCATE(sf)
 996    ABI_DEALLOCATE(sf_corr)
 997  endif
 998 
 999  if(optwt==2.and.(option==1.or.option==3)) then
1000     green%fileprt_tau=1  ! file for G(tau) has been created here
1001  endif
1002 
1003 end subroutine print_green

m_green/printocc_green [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 printocc_green

FUNCTION

  print occupations

INPUTS

  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  green  <type(green_type)>= green function data
 option= 1 :for G(w)
         2 :for G(tau)
         3 :for G(tau) and check % G(w)
         4
         <5: write diagonal part of KS occupation matrix
         5: for G(w)
         6: for G(tau)
         7 :for G(tau) and check % G(w)
         >8: write all elements of KS occup. matrix.
         9: for G(w)

OUTPUT

PARENTS

      m_dmft,impurity_solve,m_green,qmc_prep_ctqmc

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

622 subroutine printocc_green(green,option,paw_dmft,pawprtvol,opt_weissgreen,chtype)
623 
624  use m_crystal, only : crystal_t
625  use m_oper, only : print_oper
626  use m_matlu, only : diff_matlu,print_matlu,trace_matlu
627  use m_paw_dmft, only : paw_dmft_type
628 
629 !This section has been created automatically by the script Abilint (TD).
630 !Do not modify the following lines by hand.
631 #undef ABI_FUNC
632 #define ABI_FUNC 'printocc_green'
633 !End of the abilint section
634 
635  implicit none
636 
637 !Arguments ------------------------------------
638 !type
639  type(paw_dmft_type), intent(in) :: paw_dmft
640  type(green_type),intent(in) :: green
641  integer,intent(in) :: option,pawprtvol
642  integer,optional,intent(in) :: opt_weissgreen
643  character(len=*), optional, intent(in) :: chtype
644 
645 !local variables-------------------------------
646  character(len=500) :: message
647  integer :: optweissgreen,i_tau
648 ! *********************************************************************
649  if(present(opt_weissgreen)) then
650    optweissgreen=opt_weissgreen
651  else
652    optweissgreen=2
653  endif
654 
655 
656  if(mod(option,4)==1) then
657    if(optweissgreen==2) then
658      if(present(chtype)) then
659        write(message,'(4a)') ch10,"  == The ",trim(chtype)," occupations are  == "
660      else
661        write(message,'(2a)') ch10,"  == The occupations (integral of the Green function) are  == "
662      endif
663    else if(optweissgreen==1) then
664      write(message,'(2a)') ch10,"  == The integrals of the Weiss function are  == "
665    endif
666    call wrtout(std_out,message,'COLL')
667    call print_oper(green%occup,option,paw_dmft,pawprtvol)
668  endif
669  if(mod(option,4)>=2) then
670    if(optweissgreen==2) then
671      write(message,'(2a)') ch10,"  == The occupations (value of G(tau) for tau=0-) are  == "
672    else if(optweissgreen==1) then
673      write(message,'(2a)') ch10,"  == Values of G_0(tau) for tau=0- are  == "
674    endif
675    call wrtout(std_out,message,'COLL')
676    call print_oper(green%occup_tau,option,paw_dmft,pawprtvol)
677 !   write(message,'(2a)') ch10," == check: occupations from Green functions are  == "
678 !   call wrtout(std_out,message,'COLL')
679 !   call print_oper(green%occup,1,paw_dmft,pawprtvol)
680    if(mod(option,4)>=3) then
681      call diff_matlu("Local occup from integral of G(w) ",&
682 &       "Local occup from G(tau=0-) ",&
683 &       green%occup%matlu,green%occup_tau%matlu,paw_dmft%natom,1,tol4)
684      write(message,'(2a)') ch10,&
685 &     '  *****  => Calculations of occupations in omega and tau spaces are coherent ****'
686      call wrtout(std_out,message,'COLL')
687    endif
688  endif
689 
690  if(present(chtype)) then
691    if(paw_dmft%prtvol>=4.and.&
692 &      (chtype=="DMFT (end of DMFT loop)".or.chtype=="converged DMFT")&
693 &      .and.green%occup%has_opermatlu==1) then
694      write(message,'(4a)') ch10,"  == The DFT+DMFT occupation matrix for correlated electrons is == "
695      call wrtout(ab_out,message,'COLL')
696      call print_matlu(green%occup%matlu,paw_dmft%natom,pawprtvol,opt_ab_out=1)
697      write(message,'(a)') "  "
698      call wrtout(ab_out,message,'COLL')
699    endif
700  endif
701 
702  if(mod(option,4)>=2) then
703    i_tau=1
704    if(optweissgreen==1) i_tau=-1
705    call trace_matlu(green%occup_tau%matlu,paw_dmft%natom,itau=i_tau)
706  endif
707 
708 end subroutine printocc_green

m_green/spline_fct [ Functions ]

[ Top ] [ m_green ] [ Functions ]

NAME

 spline_fct

FUNCTION

 Do fourier transformation from matsubara space to imaginary time
 (A spline is performed )

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

  opt_spline = option for the spline
              -1 log to linear.
               1 linear to log.
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data

OUTPUT

SIDE EFFECTS

  fw= function is frequency space
  ft= function is time space

PARENTS

CHILDREN

      loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

2881 subroutine spline_fct(fw1,fw2,opt_spline,paw_dmft)
2882 
2883  use defs_basis
2884  use m_paw_dmft, only : paw_dmft_type, construct_nwli_dmft
2885  use m_splines
2886 
2887 !This section has been created automatically by the script Abilint (TD).
2888 !Do not modify the following lines by hand.
2889 #undef ABI_FUNC
2890 #define ABI_FUNC 'spline_fct'
2891 !End of the abilint section
2892 
2893  implicit none
2894 
2895 !Arguments ------------------------------------
2896 !type
2897  integer,intent(in) :: opt_spline
2898  type(paw_dmft_type), intent(in) :: paw_dmft
2899  complex(dpc), intent(inout) :: fw1(:)
2900  complex(dpc), intent(inout) :: fw2(:)
2901  integer :: size_fw1
2902  integer :: size_fw2
2903  real(dp), allocatable :: omega_li(:)
2904 
2905 ! *********************************************************************
2906 
2907  size_fw1 = size(fw1)
2908  size_fw2 = size(fw2)
2909 
2910 ! == inverse fourier transform
2911  if(opt_spline==-1) then
2912 
2913 !   allocate(fw1(0:paw_dmft%dmft_nwlo-1))
2914    if(paw_dmft%dmft_log_freq==1) then
2915      ABI_ALLOCATE(omega_li,(1:size_fw2))
2916      call construct_nwli_dmft(paw_dmft,size_fw2,omega_li)
2917      call spline_c(size_fw1,size_fw2,paw_dmft%omega_lo(1:size_fw1),&
2918 &                   omega_li(1:size_fw2),fw2,fw1)
2919      ABI_DEALLOCATE(omega_li)
2920    else
2921      fw2=fw1
2922    endif
2923 
2924 ! == direct fourier transform
2925  else if(opt_spline==1) then
2926 
2927 
2928    if(paw_dmft%dmft_log_freq==1) then
2929      ABI_ALLOCATE(omega_li,(1:size_fw2))
2930      call construct_nwli_dmft(paw_dmft,size_fw2,omega_li)
2931      call spline_c(size_fw2,size_fw1,omega_li(1:size_fw2),&
2932 &                 paw_dmft%omega_lo(1:size_fw1),fw1,fw2)
2933      ABI_DEALLOCATE(omega_li)
2934    else
2935      fw1=fw2
2936    endif
2937 
2938  endif
2939 
2940 end subroutine spline_fct

newton/compute_nb_elec [ Functions ]

[ Top ] [ newton ] [ Functions ]

NAME

 compute_nb_elec

FUNCTION

  Compute nb of electrons as a function of Fermi level

INPUTS

 fermie       : input of energy
 opt_noninter

 OUTPUTS
 Fx           : Value of F(x).
 nb_elec_x    : Number of electrons for the value of x

PARENTS

      newton

CHILDREN

      compute_green,integrate_green

SOURCE

3812 subroutine compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
3813 &  Fx,opt_noninter,nb_elec_x,fermie)
3814 
3815  use m_abicore
3816 
3817  use defs_basis
3818  use defs_abitypes
3819  use m_errors
3820  use m_crystal, only : crystal_t
3821  use m_paw_dmft, only: paw_dmft_type
3822  use m_self, only : self_type
3823 
3824 !This section has been created automatically by the script Abilint (TD).
3825 !Do not modify the following lines by hand.
3826 #undef ABI_FUNC
3827 #define ABI_FUNC 'compute_nb_elec'
3828 !End of the abilint section
3829 
3830  implicit none
3831 !Arguments ------------------------------------
3832 !scalars
3833  type(crystal_t),intent(in) :: cryst_struc
3834  type(green_type),intent(inout) :: green
3835  !type(MPI_type), intent(in) :: mpi_enreg
3836  type(paw_dmft_type), intent(inout) :: paw_dmft
3837  type(pawang_type),intent(in) :: pawang
3838  type(self_type), intent(inout) :: self
3839  integer,intent(in) :: opt_noninter
3840  real(dp),intent(in)  :: fermie
3841  real(dp),intent(out) :: Fx,nb_elec_x
3842 ! *********************************************************************
3843    paw_dmft%fermie=fermie
3844    call compute_green(cryst_struc,green,paw_dmft,pawang,0,self,opt_self=1,&
3845 &   opt_nonxsum=1,opt_nonxsum2=1)
3846    call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=0,&
3847 &   opt_ksloc=-1) !,opt_nonxsum=1)
3848 !  opt_ksloc=-1, compute total charge
3849    nb_elec_x=green%charge_ks
3850    Fx=nb_elec_x-paw_dmft%nelectval
3851 
3852    if(opt_noninter==1) then
3853    end if
3854  end subroutine compute_nb_elec

newton/function_and_deriv [ Functions ]

[ Top ] [ newton ] [ Functions ]

NAME

  function_and_deriv

FUNCTION

  Compute value of a function and its numerical derivatives

INPUTS

  x_input      : input of x
  option       : if 1 compute only first derivative
                 if 2 compute the two first derivatives.
  opt_noninter

 OUTPUTS
  Fx           : Value of F(x)
  Fxprime      : Value of F'(x)
  Fxdouble     : Value of F''(x)

PARENTS

      newton

CHILDREN

      compute_green,integrate_green

SOURCE

3702 subroutine function_and_deriv(cryst_struc,f_precision,green,iter,paw_dmft,pawang,self&
3703 & ,x_input,x_old,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option)
3704 
3705  use m_abicore
3706 
3707  use defs_basis
3708  use defs_abitypes
3709  use m_errors
3710  use m_crystal, only : crystal_t
3711  use m_paw_dmft, only: paw_dmft_type
3712  use m_self, only : self_type
3713 
3714 !This section has been created automatically by the script Abilint (TD).
3715 !Do not modify the following lines by hand.
3716 #undef ABI_FUNC
3717 #define ABI_FUNC 'function_and_deriv'
3718 !End of the abilint section
3719 
3720  implicit none
3721 !Arguments ------------------------------------
3722 !scalars
3723  type(crystal_t),intent(in) :: cryst_struc
3724  type(green_type),intent(inout) :: green
3725  !type(MPI_type), intent(in) :: mpi_enreg
3726  type(paw_dmft_type), intent(inout) :: paw_dmft
3727  type(pawang_type),intent(in) :: pawang
3728  type(self_type), intent(inout) :: self
3729  integer,intent(in) :: iter,opt_noninter,option
3730  real(dp),intent(inout)  :: f_precision,x_input,x_precision
3731  real(dp),intent(out) :: Fx,Fxprime,Fxdouble
3732  real(dp),intent(inout) :: x_old
3733 !Local variables-------------------------------
3734  real(dp) :: deltax,nb_elec_x,Fxmoins,Fxplus,xmoins,xplus,x0
3735  character(len=500) :: message
3736 ! *********************************************************************
3737 
3738 !  choose deltax: for numeric evaluation of derivative
3739    if(iter==1) then
3740 !    deltax=0.02
3741    end if
3742 !  deltax=max((x_input-x_old)/10.d0,min(0.00001_dp,x_precision/100_dp))
3743    deltax=min(0.00001_dp,x_precision/100_dp)  ! small but efficient
3744 !  endif
3745 !  write(std_out,*) "iter,x_input,deltax",iter,x_input,deltax
3746    x0=x_input
3747    xmoins=x0-deltax
3748    xplus=x0+deltax
3749 
3750    call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
3751 &   Fx,opt_noninter,nb_elec_x,x0)
3752 
3753    write(message,'(a,3f12.6)') "  - ",x0,nb_elec_x,Fx
3754    call wrtout(std_out,message,'COLL')
3755 !  write(std_out,*) "Fx", Fx
3756    if(abs(Fx)<f_precision) return
3757 
3758    call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
3759 &   Fxplus,opt_noninter,nb_elec_x,xplus)
3760 
3761    write(message,'(a,3f12.6)') "  - ",xplus,nb_elec_x,Fxplus
3762    call wrtout(std_out,message,'COLL')
3763 
3764    if(option==2) then
3765      call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
3766 &     Fxmoins,opt_noninter,nb_elec_x,xmoins)
3767 
3768      write(message,'(a,3f12.6)') "  - ",xmoins,nb_elec_x,Fxmoins
3769      call wrtout(std_out,message,'COLL')
3770    end if
3771 
3772    if(option==1) then
3773      Fxprime=(Fxplus-Fx)/deltax
3774    else if (option==2) then
3775      Fxprime=(Fxplus-Fxmoins)/(2*deltax)
3776      Fxdouble=(Fxplus+Fxmoins-2*Fx)/(deltax**2)
3777    end if
3778 !  write(std_out,*) "after computation of Fxprime",myid
3779    if(Fxprime<zero) then
3780      write(message,'(a,f12.6)') "  Warning: slope of charge versus fermi level is negative !", Fxprime
3781      call wrtout(std_out,message,'COLL')
3782    end if
3783    x_old=x_input
3784 
3785    return
3786  end subroutine function_and_deriv